fsolve not solving system of nonlinear equations

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"

[Star Strider]

What solutions did Mathematica provide?

Answer 1

[John D'Errico]

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.

[bloodtalon]

Your answer is good. Thanks. However,

1) I usually avoid syms since it gives long answers but I guess using a double(res.A) works fine.

2) Since you clearly didn't realize, (1e17,1.282e+05) was an answer MATLAB gave me when I used my initial guess as (1e17,1e5). Mathematica gave me your answer, the right answer, which is clearly correct to 16 digits too.

3) Hand calculation? That's so nonsensical. Why would I do all these huge computations by hand when

a) it's time consuming

b) I have to correct it each time I have a tiny mistake, which it turns out I did

c) have to type in the answers each time into the computer for further processing

d) is just not the point of this question?

[John D'Errico]

I showed you how to solve it TRIVIALLY using MATLAB with the symbolic toolbox. Did you bother to look?

In fact, I showed you how to solve it in two different ways. That the second way was using pencil and paper is my choice. It was a simple way to prove that the answer I got was consistent, since you did not bother to provide what you felt was the true solution.

As for the solution, you stated, clearly and in your own words:

"A and E_a are 1e+17 and 1.282e+05 respectively"

Should I have been able to read your mind when you made that statement that this is what you got from fsolve, instead of what you got from Mathematica? You stated that as if it was the true values of those variables. Sorry, but it is so difficult to read your mind. I can only read what you wrote.

As far as using hand computations, feel free to forget completely how to do anything by hand. That will be your loss one day. If you don't use it, you lose it. Worse, if you rely on a computer for everything, and never know how to check the results, then expect to be surprised when one day garbage comes from your computations. (I hope you are not involved with a Mars expedition when that happens.) Clearly, even here you thought you could just throw a computation at fsolve without really understanding what was involved. I.e. without understanding floating point arithmetic, and the limitations thereof.