Using lsqcurvefit with fsolve
2 views (last 30 days)
Show older comments
Hi everyone,
I am trying to solve a nonlinear equation using fsolve. The equation is:
for i = 1:length(t)
sol(i) = fsolve(@(x) (Param + x)-(K5/x + 2.*(ydata3(i).*K3.*K4./(x^2 + x.*K3 + K3.*K4)) +...
(ydata3(i).*x.*K3./(x^2 + x.*K3 + K3.*K4))+ ...
2.*(Param.*K1.*K2./(x^2 + x.*K1 + K1.*K2)) +...
(ydata3(i).*x.*K1./(x^2 + x.*K1 + K1.*K2))+ ...
(ydata2(i).*K6)./(K6 + x)), x0);
end
However, the parameter, Param, is not known, I have to regress it from lsqcurve fitting using:
Param0 = 0.035;
Param = lsqcurvefit(@chargeBal,Param0,t,ydata1);
Where the YDATA will have the size of 4 x 1, whereas the Param, has a size of 1 x 1.
How can I make this to work?
The complete code is:
function pH_Example
global t pH K1 K2 K3 K4 K5 K6 ydata1 ydata2 ydata3 Param sol
pH = [11.48
6.86
6.72
6.25];
t = [0
30
60
100];
ydata3 = [2.58E-01
6.28E-02
5.76E-02
4.83E-02];
ydata2 = [ 1.25E-02
9.75E-03
7.93E-03
5.83E-03];
ydata1 = [3.31131E-12
1.38038E-07
1.38038E-07
1.38038E-07];
K1 = 1.09e-6;
K2 = 3.69e-10;
K3 = 6.2e-4;
K4 = 3.11e-8;
K5 = 3.02e-14;
K6 = (1.91e-2 + 5.38e-4)/2;
function H = chargeBal(Param,t)
x0 = 10^(-11.48);
for i = 1:length(t)
sol(i) = fsolve(@(x) (Param + x)-(K5/x + 2.*(ydata3(i).*K3.*K4./(x^2 + x.*K3 + K3.*K4)) +...
(ydata3(i).*x.*K3./(x^2 + x.*K3 + K3.*K4))+ ...
2.*(Param.*K1.*K2./(x^2 + x.*K1 + K1.*K2)) +...
(ydata3(i).*x.*K1./(x^2 + x.*K1 + K1.*K2))+ ...
(ydata2(i).*K6)./(K6 + x)), x0);
end
function sol = pHfun(t,x)
K1 = 1.09e-6;
K2 = 3.69e-10;
K3 = 6.2e-4;
K4 = 3.11e-8;
K5 = 3.02e-14;
K6 = (1.91e-2 + 5.38e-4)/2;
x0 = 10^(-11.48);
t = [0
30
60
100];
pH = [11.48
6.86
6.72
6.25];
sol = (Param + x)-(K5./x + 2.*(ydata3.*K3.*K4./(x^2 + x.*K3 + K3.*K4)) +...
(ydata3.*x.*K3./(x^2 + x.*K3 + K3.*K4))+ ...
2.*(Param.*K1.*K2./(x^2 + x.*K1 + K1.*K2)) +...
(ydata3.*x.*K1./(x^2 + x.*K1 + K1.*K2))+ ...
(ydata2.*K6)./(K6 + x));
end
H = sol(i);
end
Param0 = 0.035;
Param = lsqcurvefit(@chargeBal,Param0,t,ydata1);
fprintf(1,'\tParameters:\n')
for k1 = 1:length(Param)
fprintf(1, '\t\tParam(%d) = %8.5f\n', k1, Param(k1))
end
ph = 3 - log10(sol(i));
plot(t,ph,'r--')
hold on
plot (t,pH,'kd')
xlabel('Time (sec)')
ylabel('pH')
hold off
legend('pH-Mod','pH-Exp')
end
3 Comments
Answers (0)
See Also
Categories
Find more on Linear and Nonlinear Regression in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!