Trouble finding more than one solution with solve

4 views (last 30 days)
Timothy Le
Timothy Le on 26 Apr 2022
Commented: Timothy Le on 16 May 2022
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

Accepted Answer

Paul
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
S = 191×3
-0.4791 -0.0008 0.9794 -0.4769 -0.0007 0.9773 -0.4747 -0.0007 0.9752 -0.4725 -0.0006 0.9730 -0.4702 -0.0006 0.9709 -0.4679 -0.0006 0.9687 -0.4656 -0.0005 0.9665 -0.4633 -0.0005 0.9643 -0.4609 -0.0005 0.9621 -0.4586 -0.0005 0.9598
plot(T_star_val,S,'bo')
  17 Comments
Timothy Le
Timothy Le on 16 May 2022
I ended up using your method of solving the equation in my code. My sincere thanks for your help with this!
Also thanks to @Walter Roberson and @Star Strider for their inputs, I really appreciated it!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!