Using lsqcurvefit with fsolve
3 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!