MATLAB: Fsolve not solving system of nonlinear equations

equationfsolvefunctionsMATLABnonlinearoptimizationsolvesystem

Have a system of 2 nonlinear equations which Mathematica solved in 5 seconds but MATLAB is having trouble with for some reason. Probably due to the large numbers, which perhaps result in very small gradients. Tried changing the tolerance but it didn't seem to work. And different initial guesses gave very different answers for one of the variables (the second), while the other could never be found.
function Q1
clear;clc;
fun = @root2d;
x0 = [1e13,1e5];
%options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt);
options = optimoptions('fsolve','Display','iter','FunctionTolerance',1e-20);
x = fsolve(fun,x0,options);
fprintf('A and E_a are %.4g and %.4g respectively\n',x);
function F = root2d(x)
% x(1) is A and x(2) is E_a
F(1)=x(1).*exp(-x(2)./(8.31447*(50+273.15)))-0.005;
F(2)=x(1).*exp(-x(2)./(8.31447*(100+273.15)))-0.1;
The results, with some initial guesses, are "A and E_a are 1e+15 and 1.141e+05 respectively", "A and E_a are 1e+17 and 1.282e+05 respectively"

Best Answer

That mathematica was able to solve it is just a reflection that mathematica uses high precision. MATLAB works in double precision, which is much faster for many problems, but is limited.
So had you made the effort to properly scale things so the variables are more alike, fsolve would have been successful too. Or would it have been? As I will show, the problem lies more deeply, in that you have not solved the same problem as claim to have solved in Mathematica.
Personally, I'd just use pencil and paper. We have
x(1).*exp(-x(2)./(8.31447*(50+273.15))) = 0.005
Take the log of each equation.
log(x(1)) - x(2)./(8.31447*(50+273.15)) = log(0.005)
log(x(1)) - x(2)./(8.31447*(100+273.15))) = log(0.1);
Subtract, solve for x(2).
x(2)*(1/(8.31447*(50+273.15)) - 1/(8.31447*(100+273.15))) = log(0.1) - log(0.005)
Therefore, x(2) is:
x2 = (log(0.1) - log(0.005))/(1/(8.31447*(50+273.15)) - 1/(8.31447*(100+273.15)))
x2 =
        60069.6595700855
Now, solve for x1, given x2.
x1 = exp(x2./(8.31447*(100+273.15))) * 0.1
x1 =
        25618686.6641389
Could MATLAB have solved it directly? Of course. I'd not use fsolve though, because of the dynamic range of the variables. Instead, use something comparable to what Mathematica would have done.
syms x1 x2
F1=x1.*exp(-x2./(8.31447*(50+273.15)))-0.005 == 0;
F2=x1.*exp(-x2./(8.31447*(100+273.15)))-0.1 == 0;
res = solve(F1,F2)
res.x1
ans =
20^(2954190909812263/457092822189736)/10
res.x2
ans =
(10077583391870757373634230713737*log(20))/502578872970562362079707136
vpa(res.x1)
ans =
25618686.664139188215079770087505
vpa(res.x2)
ans =
60069.65957008549587848872243846
Which as you can see, are essentially the same to the full precision that I got from my hand computations.
Are these the same results that you claim Mathematica yielded? OF COURSE NOT! But what you claim to be the solutions simply do not solve the equations you have given us. So you apparently solved a different problem in Mathematica. Do I need to convince you of that fact? Try substituting the supposedly correct solutions you have into those equations.
vpa(subs(F1,{x1,x2},{1e17,1.282e+05}))
ans =
-0.0048103697511556315807357424524785 == 0.0
As you can see, this is clearly incorrect. Whereas...
vpa(subs(F1,{x1,x2},{25618686.6641389,60069.6595700855}))
ans =
-0.000000000000000070252340008864116230491806750275 == 0.0
Which agrees to within the 16 digits or so that I put in.
The moral is, if you will make the effort to compare software, at least try to solve the same problem in each tool.
Related Question