Difficulty applying fsolve in this context
1 view (last 30 days)
Show older comments
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.
3 Comments
Accepted Answer
Matt J
on 13 Aug 2014
Edited: Matt J
on 13 Aug 2014
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).
2 Comments
Matt J
on 13 Aug 2014
Edited: Matt J
on 13 Aug 2014
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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!