Nonlinear equation solver issue

49 views (last 30 days)
Rick
Rick on 26 Nov 2014
Commented: John D'Errico on 26 Nov 2014
Hello, I have an equation
N = (DP/RT * ln(P/(P-P*)))/(L/A) + (l/a))
P,R,T,P*,(L+l), and A are known. L+l is always .05
There is some experimental data
l = .004, N = 6.6e-11
l = .008, N = 4.9e-11
l = .012, N = 4.0e-11
I am trying to find D and a. So I make a function and I try to use fsolve, but when I do I just get an erroneous answer. I am saying D(1) = D, and D(2) = a in my function
function F = stefanflow(D)
P = 1e5;
T = 294;
R = 8.314e3;
A = 1e-4;
P_sat = 1.12e4;
O = log(P/(P-P_sat));
I = P/(R*T);
F = zeros(3,1);
F(1) = 6.6e-11 - (D(1)*I*O)/(((.05-.004)/A)+(.004/D(2)));
F(2) = 4.9e-11 - (D(1)*I*O)/(((.05-.008)/A)+(.008/D(2)));
F(3) = 4.0e-11 - (D(1)*I*O)/(((.05-.012)/A)+(.012/D(2)));
I run this in my script file
options = optimoptions('fsolve','Display','iter');
[D,fval] = fsolve(@stefanflow,[1e-5;1e-5],options);
And this is my output
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle
non-square systems; using Levenberg-Marquardt algorithm instead.
> In fsolve at 286
First-Order Norm of
Iteration Func-count Residual optimality Lambda step
0 3 2.59598e-22 1.19e-16 0.01
1 6 2.59598e-22 1.19e-16 0.001 1.37935e-14
Equation solved, fsolve stalled.
fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance and the vector of function values
is near zero as measured by the default value of the function tolerance.
<stopping criteria details>
EDU>> D
D =
1.0e-04 *
0.1000
0.1000
What is going on, and how can I solve this system?
  1 Comment
John D'Errico
John D'Errico on 26 Nov 2014
Fsolve is happy with the starting values. It told you that fact. They meet the convergence criteria, so why look beyond that point?

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Matt J
Matt J on 26 Nov 2014
Edited: Matt J on 26 Nov 2014
The solution fsolve gave you appears to satisfy the equations very well,
>> stefanflow(D)
ans =
1.0e-11 *
0.9493
0.9167
0.9243
So, either you have incorrectly coded your equations, or they have multiple solutions that you don't expect.
You could also try making the tolerances more strict,
options = optimoptions('fsolve','Display','iter','TolFun',1e-30,'TolX',1e-30);
When I try this, it produces a small change in the solution,
D =
1.0e-04 *
0.1013
0.1391
but maybe that's significant to you.
  3 Comments
Show 1 older comment
Alan Weiss
Alan Weiss on 26 Nov 2014
Edited: Alan Weiss on 26 Nov 2014
I only see two elements in D. The first is 1e-4*0.1013, the second is 1e-4*0.1391.
Alan Weiss
MATLAB mathematical toolbox documentation
John D'Errico
John D'Errico on 26 Nov 2014
Rick - I think you miss that the first "number" shown is not an element of D, but an exponent applied to those numbers. Note the little * after that number.

Sign in to comment.

More Answers (0)

Categories

Find more on Choose a Solver in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!