Using lsqcurvefit with fsolve
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Linear and Nonlinear Regression finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!