# fsolve function give poor results for multiple equations

3 views (last 30 days)
Peng Liu on 20 Oct 2020
Commented: Peng Liu on 21 Oct 2020
I want to get 7 parameters used in three different equations (let's say equation CP, S, and H), however, I am only happy for the fitted parameters for H equation. The fitted results for CP and S is really poor. How to improve it? What I wish is that the fitted R square is minimun for the three equations. The code is listed below. Thanks in advance.
function y = nonlinmodel(x,Temp,CP,S,H)
y = [CP - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
S - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
H - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];
end
Temp=[1200 1300 1500 1600 1700 1800 2000 2200 2500];
S=[520.042428000000 567.730080000000 582.090804000000 595.781640000000 608.844456000000 633.337236000000 655.862220000000 686.593332000000 731.182752000000];
CP=[206.367372000000 220.686228000000 224.245008000000 227.343240000000 230.064660000000 234.586404000000 238.145184000000 242.164512000000 246.644388000000];
H=[-121082.256000000 -56898.6120000000 -34624.8360000000 -12057.9840000000 10801.9440000000 57275.4240000000 104586.264000000 176641.092000000 298937.520000000];
f = @(x)nonlinmodel(x,Temp,CP,S,H);
optimal_x = fsolve(f,[1;0.2;0.03;0.004;0.00005;6;72]);

Matt J on 20 Oct 2020
Edited: Matt J on 20 Oct 2020
Your equations are linear, so there is no reason to be using fsolve. I was able to obtain the coefficient matrix A for your linear system using func2mat from the File Exchange,
A=full(func2mat( @(x)CoeffFun(x,Temp,CP,S,H) ,zeros(7,1))),
function y = CoeffFun(x,Temp,CP,S,H)
Temp=Temp(:);
y = -[ - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
- 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
- 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];
end
The condition number for you A matrix is very poor, so it's not too surprising that you don't get a good result.
>> cond(A)
ans =
3.2818e+19
I wonder if Temp is expressed in the correct units? The current units result in an insane order of magnitude for the coefficients,
>> max(abs(A(:)))
ans =
1.6238e+17

Show 6 older comments
Peng Liu on 21 Oct 2020
I got the func2mat function from your shared link. However, new error comes "Undefined function or variable 'vecnorm'".
Matt J on 21 Oct 2020
If you can, upgrade to R2017b
Otherwise,
a=sqrt(sum(A.^2,1)) %vecnorm for older Matlab versions
Peng Liu on 21 Oct 2020
Yeah, It works now. Thank you!
Having a nice day!

### More Answers (1)

Alan Weiss on 20 Oct 2020
I wasn't able to find a very good answer either. I'm not sure that one exists. I think that lsqnonlin is a more appropriate solver. Here is what I did using MultiStart to search:
Temp=[1200 1300 1500 1600 1700 1800 2000 2200 2500];
S=[520.042428000000 567.730080000000 582.090804000000 595.781640000000 608.844456000000 633.337236000000 655.862220000000 686.593332000000 731.182752000000];
CP=[206.367372000000 220.686228000000 224.245008000000 227.343240000000 230.064660000000 234.586404000000 238.145184000000 242.164512000000 246.644388000000];
H=[-121082.256000000 -56898.6120000000 -34624.8360000000 -12057.9840000000 10801.9440000000 57275.4240000000 104586.264000000 176641.092000000 298937.520000000];
f = @(x)nonlinmodel(x,Temp,CP,S,H);
[optimal_x, resnorm] = lsqnonlin(f,[1;0.2;0.03;0.004;0.00005;6;6]);
lb = [-5000;-2;-.01;-.01;-.01;-10;0];
ub = [2;2;.01;.01;.01;10;20000];
ms = MultiStart;
x0 = lb + rand(size(lb)).*(ub - lb);
prob = createOptimProblem('lsqnonlin','x0',x0,'objective',f);
prob.lb = lb;
prob.ub = ub;
[x,fval,exitflag,output,solutions] = run(ms,prob,1000);
disp(x(2:6))
function y = nonlinmodel(x,Temp,CP,S,H)
y = [CP - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
S - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
H - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];
end
The smallest residual I found was near 1e9. Not so good.
Alan Weiss
MATLAB mathematical toolbox documentation

#### 1 Comment

Peng Liu on 21 Oct 2020
Dear Alan,
I am very appreciate for your feedback. Is it possible to constrain the optimazation by minimizing the sum of squared difference between the estimated values and the actual values? I find that this constrain condition is very useful in Excel for my case.