Trouble finding more than one solution with solve

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

Answer 1

[Paul]

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')

[Timothy Le]

Hi @Paul

Thanks for this solution! From the plot you gave, it seemed really promising. However, when I ran the code for myself and had a look at Sd, I could see that there was a point where it failed to find all three solutions when there are meant to be three and at others there were even 4 solutions.

Also, how would you go about separating the three solutions into three lines as I have in the plot in this comment?

[Paul]

Hi Timothy,

If any value of T_star yielded four solutions, then S (not Sd) would have four columns. But S only has three columns. Can you tell me the value of T_star where you see four solutions?

Similarly, which value of T_star did not find all three solutions? This code looks for sign changes in y1, so it will miss a case where an element of y1 is exactly equal to zero. So you'd have to look for that special case.

Also, with a brute force approach such as this, it's possible to miss a zero if OPval isn't dense enough. I was original using a spacing of 1e-6, but that was taking forever for erfi() to run.

Did you link to the correct comment? In any case, for any value of T_star you'll have to figure out how to sort the values of OPzeros so they are consistent with results from the previous value of T_star. Given the apparent structure of the of the solution, that might not be too difficult, beause one root seem so to be always zero, and the others are monotic if they exist. I vaguely recall a post here not too long ago that considered a very similar, if not identical, problem and had a very clever solution. Don't recall much more about it than that.

[Timothy Le]

Hi Paul,

Oh okay, my apologies! I thought that I should be looking at Sd. I had at look at S and everything seems fine apart from getting NaN at the point where the blue dotted line on the plot from my comment on the previous answer crosses the black line. I'd assume that the solution at this point is exactly 0 so I'll take your explanation for it and manually insert the solution I guess.

In that case, looking at S, I don't think it will be too difficult to figure out a way to sort the values.

Thanks again for your solution!

[Walter Roberson]

At OP = 0 exactly, the Ms(i) is 0 / 0 which give NaN. But if you take left and right limit, then from both directions the result is 0, so for that case in theory you can substitute 0 for that NaN.

[Paul]

Nothing to apologize for. I'm still confused about references to Sd here, as my code doesn't have an Sd defined.

Anyway, please keep in mind that the solution doesn't return NaN, at least not in the way that I think of that. I was just using the NaN values to initalize S so that a) I woulnd't have to use a cell array in case a different number of OPzeros were returned for different values of T_star, and b) so it would be easy to see how many solutions were actually returned for different values of T_star, and c) so it would be easy to plot the result because plot() ignores NaN values.

As for the situation at OP = zero, I'm not convinced that OP = 0 is a root of MS for all values of T_star and a double root at T_star = 1. As @Walter Roberson pointed out, the limit of MS may be 0 as OP approaches 0, but I'm not sure that means OP is a root, unless you were to explicity define MS that way, e.g.,

MS(OP) = { 0 for OP = 0; etc, for OP ~=-0}

Regardless, if we assert the root at OP = 0, there's no need to solve for it numerically and the entries in S above that correspond to the root at OP = 0 can be fixed accordingly. If you look at the row of S that corresponds to T_star = 1, you should see that my code didn't even find the root at OP = 0. I think that's because T_star = 1 is a special case where MS doesn't change sign on either side of OP = 0.

The prototype code above uses linear interpolation to approximate the solution. No reason a few more points couldn't be used with cubic or spline interpolation if you see a need to do so.

[Paul]

@Walter Roberson

Did you mean "in theory" in the mathemetical sense, i.e., that there is a theory that applies in this situation? Or did you mean it in the sense that because the limit is zero, we can just define MS(0) = 0 and get on with things? I guess it depends on whether or not @Timothy Le cares about what happens at OP = 0.

[Walter Roberson]

Every closed compact set near OP = 0 gives 0, for every finite value of T_star , and there is no specific definition of MS at OP = 0 that would disagree with 0. Complex Analysis say that means that the limit in that situation is 0.

[Timothy Le]

@Paul

Sorry for confusing you! It was my mistake! I had some leftover variables from when I was running the codes for previous solutions. Forgive me!

Also, I have tried running the code for

V0 = -0.001;
T_star_min = 0.001;
T_star_max = 2.000;
interval = 0.001;

And for some reason, I'm not sure why, I get NaN values in the first 10 rows of S. Do you have an idea of what could be causing this?

[Paul]

I think that func isn't capable of accurate computations for some values of T_star with V0 = -0.001. Didn't @Walter Roberson mention this risk somewhere in this thread?

Regardless, consider

V0 = -0.001;
% 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);

OPval = -1:1e-2:1;

% plots with T_star = 0.001

%this looks ok
figure
plot(OPval,double(vpa(subs(MS,[OP T_star],{OPval , 0.001})))),grid

% this doesn't
figure
plot(OPval,(func(OPval,0.001))),grid

Warning: Imaginary parts of complex X and/or Y arguments ignored.

% plots with T_star = 0.1

% both look ok
figure
plot(OPval,double(vpa(subs(MS,[OP T_star],{OPval , 0.1})))),grid

figure
plot(OPval,(func(OPval,0.1))),grid

Note that OP = 0 is no longer a singularity

subs(MS,[OP T_star],[0 0.1])

ans = 

vpa(ans)

ans = 

0.0020028514026002125900561128929534

Did the singularity move to some other value of OP, or is it gone altogether with V0 ~= 0?

Anyway, these results suggest an alternative approach that doesn't use matlabFunction

T_star_min = 0.001;
T_star_max = 2.0;

I had to change the interval to meet run time constraint on Answers

interval = 0.001;
%T_star_val = T_star_min:interval:T_star_max;
T_star_val = [T_star_min:interval:0.01 .02:.01:T_star_max];
S = nan(numel(T_star_val),3);
for ii = 1:numel(T_star_val)
  % convert to double here, but could do all of this in VPA math if desired?  
  y1 = double(vpa(subs(MS,[OP T_star],{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 = 209×3

-0.4998 -0.0031 0.9998 -0.4996 -0.0018 0.9996 -0.4994 -0.0013 0.9994 -0.4992 -0.0010 0.9992 -0.4990 -0.0009 0.9990 -0.4988 -0.0007 0.9988 -0.4986 -0.0007 0.9986 -0.4984 -0.0006 0.9984 -0.4982 -0.0005 0.9982 -0.4980 -0.0005 0.9980

plot(T_star_val,S,'bo')

I did not explore what's happening at T_star = 1.

[Walter Roberson]

Those erfi() lead to quite inaccurate computations when you use matlabFunction()

[Timothy Le]

@Paul

I don't recall him mentioning this and had a quick look and didn't find anything. Perhaps I missed it.

I don't know much about the singularity moving, part of my dissertation is finding out how V0 affects things so I'm going into this with no idea what to expect haha. That's the reason why I wanted a better code for solving the equation because the code that I wrote took way too long (more than 24 hours to run!) when V0 ~= 0.

I don't think what's happening at T_star = 1 is important. It's probably beyond the scope of my dissertation. I just need to see how it changes and discuss it really and then compare it to simulations.

I just ran your code for an interval of 0.001 and it seems to have worked fine, thank you!

[Timothy Le]

@Walter Roberson

Would it be better to vpa the function to get rid of erfi then? Sorry if this is a silly question, I'm not much of an expert when it comes to MATLAB :)

[Walter Roberson]

See the "Note: several paragraphs down in https://www.mathworks.com/matlabcentral/answers/1705785-trouble-finding-more-than-one-solution-with-solve#comment_2129175

[Walter Roberson]

When you have a symbolic expression, the symbolic toolbox can rewrite it to any equivalent form, and the order of evaluation is not locked in. For example

   f(x)/(A*x+b)

with constant A and b might be internally represented as (1/A) * f(x) / (x + (b/A)) but the actual order of the division might change as you know more about f(x)... and the whole internal approach sometimes gets rewritten if you have a negative!

If you end up substituting a rational value for x, the expression might get reordered. simplify() might reorder further. The question of which order the expression will "really" be evaluated is pretty difficult to answer.

When you vpa() you do not change which functions are called, but you do trigger the symbolic engine to choose some particular order of operations, and convert rational constants to software floating point (except for exponents), and try to do numeric approximations for all functions that do not have unbound variables.

vpa() would not rewrite erfi with a unbound variable into something like an infinite sum: if there is an unbound variable it would leave it as a symbolic call to erfi.

matlabFunction goes further in binding to a very specific and inflexible order of operations, to be evaluated at 53 bits. If the erfi is close to 1 or close to 0, the loss of precision by working with double precision can be crucial.

[Paul]

Does vpa() apply any simplifications before evaluation? For example:

syms f(x)
f(x) = cos(x)^2 + sin(x)^2;
vpa(f(.01))

ans = 

1.0

Does this result obtain because of the identity, or because that's just what cos^2(0.01) + sin^2(0.01) evaluate to in VPA computations?

[Walter Roberson]

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()

[Timothy Le]

@Paul

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!