Difficulty applying fsolve in this context

Hello. I have been working on trying to solve a relatively simple equation but it has been a real doozy getting MATLAB to work with it. I am definitely doing something wrong, but after a few days of messing around with solvers, I have made no headway whatsoever.

I am plotting a function whereby every point is determined by an equation whose parameters must be numerically solved for. But, I am stuck on the solving part. I have tried using both fzero and fsolve and I just end up with a jagged mess, unity, or no solution at all.

Here is my function file that is called by fsolve:

function F = maugis(m,i,lambda,A)
F = 0.5*lambda*A(i).^2*(sqrt(m.^2-1)+(m.^2-2).*atan(sqrt(m.^2-1)))+(4/3)*lambda^2*A(i).*(sqrt(m.^2-1).*atan(sqrt(m.^2-1))-m+1)-1;

Here are the values of the constants, for clarity:

w   =   .07; % Estimate for surface energy
E1  =   400e9;  % E of IFM tip
E2  =   170e9;  % E of sample
n1  =   .28;    % Poisson of IFM tip
n2  =   .22;     % Poisson of sample
R   =   7e-6;   % Radius of tip
Es   =   ((1-n1^2)/E1+(1-n2^2)/E2)^(-1);
A=linspace(.02,2,100);
s0=13.7e3; % max of Lennard-Jones-derived stress
lambda=(2*s0)/(pi*w*Es^2/R)^(1/3); % interaction parameter
a=A*((3*pi*w*R^2)/(4*Es))^(1/3);

The solver itself is embedded in a "for" loop on the main code as follows:

for i=1:length(A)
    m0=2;
    m_md(i)=fsolve(@(m) maugis(m,i,lambda,A),m0);
end

Once m_md is acquired, I run it through this code:

Pmd=A.^3-lambda*A.^2.*(sqrt(m_md.^2-1)+m_md.^2.*(tan(sqrt(m_md.^2-1))).^2);
fmd=Pmd.*(pi.*w.*R);
figure(2)
axis([-1e-5,1e-5,0,8e-8])
plot(fmd,a)

Is there a better way to set up this solver? I have plotted this function in a separate plot and I get explicit zero crossings everywhere. I am not sure what is going on. If I set the initial search point m0 as 1, I get a completely different solution than if I set it as 2, for instance. m0=1 shows a smooth line after plotting, while m0=2 generates an almost random scattering of points.

[Yu Jiang]

Hi Charles

I wonder what are the expected solutions for m_md and why would you think "zero crossing everywhere" is not correct? Should there be bounds on m_md?

- Yu Jiang

[Charles]

m_md should be unbounded, but should never be less than unity. If A is 0, then m_md would theoretically be infinity. As A increases, it should decrease monotonically towards 1.

Right now I am baffled as to the erratic behavior of the solver in how I set it up. It seems to produce different results depending on my starting value. But, having plotted the equation in a separate program, it should only cross 0 at one point.

[Charles]

Oh, I should also never have complex numbers in my plot. I think Pmd is outputting complex numbers for some reason, given these initial values. I have no idea. It's quite frustrating since I am dealing with lots of real-world numbers, but also attempting to mimic some results from some journal articles that have pretty much brushed over the explicit solution steps for the model.

Answer 1

[Matt J]

You are doing nothing to enforce the constraint that m>=1, which leads to undefined or complex values for your objective F as the solver iterates. Use lsqnonlin to add bounds on m, or maybe even just fminbnd as applied to F^2.

Also, if the m_md(i) are supposed to be varying continuously with i as the loop increments, solve for m_md(i) using an initial value of m0 = m_md(i-1).

[Matt J]

Still better might be to make the change of variables

z^2=sqrt(m.^2-1)

and solve in terms of z instead of m. This leads to an objective,

F = 0.5*lambda*A(i).^2*(z.^2+(z.^4-1).*atan(z.^2)+...
    (4/3)*lambda^2*A(i).*(z.^2.*atan(z.^2)-sqrt(1+z.^4)+1)-1;

Notice that, in contrast to F(m), the function F(z) is defined and also differentiable at all z, which is what fsolve requires.

[Charles]

Thanks for the prompt reply and the suggestions. I tried all three, and decided on combining the change of variables approach you suggested with fminbnd for F^2, along with an interval over which to iterate.

In a separate file, I plotted the F with respect to a range of z values, each curve representing a different value for A. (I ran it through a loop), just to get an idea of where the minima should show up. I get a plot that looks like this:

A ranges from 0 to 2 (Red -> Blue) in this plot. So, I should expect z values from about 10 to past 30.

Here is my modified function file for an fminbnd approach:

function F = maugis(z,i,lambda,A)
F = (0.5*lambda*A(i).^2*(z.^2+(z.^4-1).*atan(z.^2))+...
    (4/3)*lambda^2*A(i).*(z.^2.*atan(z.^2)-sqrt(1+z.^4)+1)-1)^2;

Here is my code for the altered "for" loop using fminbnd.

for i=1:length(A)
    z(i)=fminbnd(@(z) maugis(z,i,lambda,A),0,40);
end

Running the plot again, I get (as the black curve) fmd vs. a as:

Which is the same problem I had previously when using fsolve! Frustrating. However, looking at the equation for Pmd more closely... it looks like the paper I was reading had a typo. I cross-referenced with some other papers and it seems that

Pmd=A.^3-lambda*A.^2.*(sqrt(m_md.^2-1)+m_md.^2.*(tan(sqrt(m_md.^2-1))).^2);

should ACTUALLY be

Pmd=A.^3-lambda*A.^2.*(sqrt(m_md.^2-1)+m_md.^2.*atan(sqrt(m_md.^2-1)));

By doing so, I get the following plot (along with other parts of the code I activated):

It seems the culprit was a misprint in the journal article! Go figure (per se).

I still am getting some imaginary numbers for the blue data, but I figure that is another problem that I will use the skills I got sleuthing around for this bug to solve. Thanks for all the help, everyone.