Trouble finding more than one solution with solve
3 views (last 30 days)
Show older comments
I have a code that solves an equation using the solve function in MATLAB. However, the problem is it only returns 1 solution to the equation. For example, I know that the solutions to MSvpa(1) should be a value close to -0.5, one close to 0.9 and one at 0. When I run the code I only get the solution close to -0.5.
%% Variables
% User-defined
V0 = 0;
T_star_min = 0.1;
T_star_max = 2.0;
interval = 0.1;
% Calculated
T_star = T_star_min:interval:T_star_max;
n = size(T_star);
%
%% Solving the Maier-Saupe (MS) equation for the order parameter S
% Setting up the equations
syms x OP
numerator = int((3*x^2/2-1/2)*exp((5*OP-V0)*(3*x^2/2-1/2)./T_star), x, 0 ,1);
denominator = int(exp((5*OP-V0)*(3*x^2/2-1/2)./T_star), x, 0, 1);
MS = (numerator./denominator) - OP;
% Approximating imaginary error function erfi
MSvpa = vpa(MS);
% Solving for each value of T*
S = cell([1 n(2)]);
for i = 1:n(2)
S{i} = solve(MSvpa(i), OP);
end
0 Comments
Accepted Answer
Paul
on 28 Apr 2022
Here's a brute force approach that at least seems to give a reasonable result for V0 = 0. This code runs very slowly because erfi() with a double input still is evaluated via the Symbolic version of erfi().
V0 = 0;
% Calculated
%
%% Solving the Maier-Saupe (MS) equation for the order parameter S
% Setting up the equations
syms x OP T_star
numerator = int((3*x^2/2-1/2)*exp((5*OP-V0)*(3*x^2/2-1/2)./T_star), x, 0 ,1);
denominator = int(exp((5*OP-V0)*(3*x^2/2-1/2)./T_star), x, 0, 1);
MS = (numerator./denominator) - OP;
func = matlabFunction(MS);
% Solving OP for each value of T*
OPval = -1:1e-2:1;
MS as defined above has a singularity at OP = 0, so remove OP = 0 from the values to be evaluated. Note that this will cause a root to be found at OP = 0, don't know if that's correct or not.
OPval(OPval == 0) = [];
T_star_min = 0.1;
T_star_max = 2.0;
interval = 0.01;
T_star_val = T_star_min:interval:T_star_max;
S = nan(numel(T_star_val),3);
for ii = 1:numel(T_star_val)
y1 = func(OPval,T_star_val(ii));
y = [0; y1(:)];
ys= [y1(:); 0];
k = find(sign(y) .* sign(ys) == -1);
OPzeros = -y1(k-1).*(OPval(k)-OPval(k-1))./(y1(k)-y1(k-1))+OPval(k-1);
S(ii,1:numel(OPzeros)) = sort(OPzeros);
end
S
plot(T_star_val,S,'bo')
17 Comments
Walter Roberson
on 1 May 2022
Edited: Walter Roberson
on 1 May 2022
vpa() does not apply simplifications. simplifications can be applied:
- if you use rewrite()
- when using int()
- internally when you use isAlways()
- when using solve() or dsolve()
More Answers (0)
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!