Need a better guess y0 for consistent initial conditions.

6 次查看（过去 30 天）
Dursman Mchabe
Dursman Mchabe 2019-1-14
评论： Dursman Mchabe 2019-1-16
Hi all,
I get the the error message:
Error using daeic12 (line 166)
Need a better guess y0 for consistent initial conditions.
Error in ode15s (line 310)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
Error in IgorWaterODE15s/kinetics (line 58)
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in IgorWaterODE15s (line 166)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
When I run the code below:
function IgorWaterODE15s
function C =kinetics(theta,t)
% Initial conditions
%c0 = [1; 1; 1; 1];
c0 = [0;0;0;0];
%c0 = [3.16E-02;0.33;1.62E-04;1.35E-04];
% Time
t=[0
960
1920
2880
3840
4800
5760
6720
7680
8640
9600
10560
11520
12480
13440
14400
15360
16320
17280
18240
19200
20160
21120
22080
23040
24000
24960
25920
26880
27840
28800
29760
30720];
tspan = linspace(min(t), max(t),33);
tspan = tspan';
% A constant, singular mass matrix
M = [1 0 0 0
0 1 0 0
0 0 0 0
0 0 0 0];
% Options
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-6 1e-6 1e-6],'Vectorized','on');
% Solving
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
%
function dC = DifEq(t,c)
% theta(1) = 5.97E-04;
% theta(2) =1.95E-03;
% theta(3) =8.314;
% theta(4) =323.15;
% theta(5) =5.3E-8;
% theta(6) =143.40;
% theta(7) =6.24;
% theta(8) =5.68E-05;
% theta(9) =2.04E-09;
% theta(10) =2.85E-09;
% theta(11) =1.70E-09;
% theta(12) =5.00E-06;
% theta(13) =4.72E-03;
dC = [(theta(1)./theta(2)).*((1930.65./1000000)- (c(1,:) .*theta(3) .*theta(4) ./ 875000 )) - (c(1,:).* theta(3).* theta(4) - theta(6).* ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) - ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2) - ((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))).* theta(13)))
(c(1,:).* theta(3).* theta(4) - theta(6).* ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) - ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2) - ((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))).* theta(13)))
- c(3,:) + (theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) + 2.* (theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2)+ (theta(5)./c(3,:))
- c(4,:) + ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))) + 2.*((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))) + (theta(5)./c(4,:))];
end
C=Cv
end
t=[0
960
1920
2880
3840
4800
5760
6720
7680
8640
9600
10560
11520
12480
13440
14400
15360
16320
17280
18240
19200
20160
21120
22080
23040
24000
24960
25920
26880
27840
28800
29760
30720];
c= [3.16E-02 0.33 1.62E-04 1.35E-04
3.56E-02 0.64 2.69E-04 2.24E-04
3.83E-02 0.93 4.46E-04 3.72E-04
4.12E-02 1.18 2.51E-03 2.09E-03
4.34E-02 1.40 0.43 0.35
4.58E-02 1.60 1.20 1.00
4.74E-02 1.78 1.70 1.41
4.92E-02 1.94 1.99 1.66
5.08E-02 2.08 2.81 2.34
5.21E-02 2.20 2.95 2.45
5.34E-02 2.31 3.23 2.69
5.49E-02 2.41 3.31 2.75
5.56E-02 2.49 3.38 2.82
5.63E-02 2.57 3.16 2.63
5.72E-02 2.63 3.08 2.57
5.80E-02 2.69 3.38 2.82
5.88E-02 2.74 3.23 2.69
5.92E-02 2.78 3.23 2.69
5.98E-02 2.82 3.62 3.02
6.03E-02 2.85 3.71 3.09
6.06E-02 2.87 3.71 3.09
6.10E-02 2.89 3.54 2.95
6.14E-02 2.91 3.54 2.95
6.15E-02 2.93 3.71 3.09
6.18E-02 2.94 3.79 3.16
6.20E-02 2.95 3.71 3.09
6.22E-02 2.96 3.88 3.24
6.24E-02 2.96 3.88 3.24
6.25E-02 2.97 3.88 3.24
6.26E-02 2.97 4.07 3.39
6.27E-02 2.97 4.07 3.39
6.28E-02 2.98 4.16 3.47
6.29E-02 2.98 4.16 3.47];
%theta0=[1;1;1;1;1;1;1;1;1;1;1;1;1];
%theta0 = [0;0;0;0;0;0;0;0;0;0;0;0;0];
theta0 = [0;1;0;1;0;1;0;1;0;1;0;1;0];
%theta0 = [5.97E-04;1.95E-03;8.314;323.15;5.3E-8;143.40;6.24;5.68E-05;2.04E-09;2.85E-09;1.70E-09;5.00E-06;4.72E-03];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
%[theta]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t),33);
tv = tv';
Cv = kinetics(theta, tv)
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cv);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend('Exp-SO_{2,gas}','Exp-S_{total}','Exp-H^+_{liquid}','Exp-H^+_{Interface}','Mod-SO_{2,gas}','Mod-S_{total}','Mod-H^+_{liquid}','Mod-H^+_{Interface}')
end
What better guess can I try?
  36 个评论
显示 34更早的评论
Torsten
Torsten 2019-1-16
This is the code I tested, and it worked. Since I don't have the optimization toolbox, I can't use "fsolve" or "lsqcurvefit":
function main
t=[0;960;1920;2880;3840;4800;5760;6720;7680;8640;9600;10560;11520;12480;13440;14400;15360;16320;17280;18240;19200;20160;21120;22080;23040;24000;24960;25920;26880;27840;28800;29760;30720];
c= [3.16E-02 0.33 1.62E-04 1.35E-04; ...
3.56E-02 0.64 2.69E-04 2.24E-04; ...
3.83E-02 0.93 4.46E-04 3.72E-04; ...
4.12E-02 1.18 2.51E-03 2.09E-03; ...
4.34E-02 1.40 0.43 0.35; ...
4.58E-02 1.60 1.20 1.00; ...
4.74E-02 1.78 1.70 1.41; ...
4.92E-02 1.94 1.99 1.66; ...
5.08E-02 2.08 2.81 2.34; ...
5.21E-02 2.20 2.95 2.45; ...
5.34E-02 2.31 3.23 2.69; ...
5.49E-02 2.41 3.31 2.75; ...
5.56E-02 2.49 3.38 2.82; ...
5.63E-02 2.57 3.16 2.63; ...
5.72E-02 2.63 3.08 2.57; ...
5.80E-02 2.69 3.38 2.82; ...
5.88E-02 2.74 3.23 2.69; ...
5.92E-02 2.78 3.23 2.69; ...
5.98E-02 2.82 3.62 3.02; ...
6.03E-02 2.85 3.71 3.09; ...
6.06E-02 2.87 3.71 3.09; ...
6.10E-02 2.89 3.54 2.95; ...
6.14E-02 2.91 3.54 2.95; ...
6.15E-02 2.93 3.71 3.09; ...
6.18E-02 2.94 3.79 3.16; ...
6.20E-02 2.95 3.71 3.09; ...
6.22E-02 2.96 3.88 3.24; ...
6.24E-02 2.96 3.88 3.24; ...
6.25E-02 2.97 3.88 3.24; ...
6.26E-02 2.97 4.07 3.39; ...
6.27E-02 2.97 4.07 3.39; ...
6.28E-02 2.98 4.16 3.47; ...
6.29E-02 2.98 4.16 3.47];
theta0 = [5.97E-04;1.95E-03;8.314;323.15;5.3E-8;143.40;6.24;5.68E-05;2.04E-09;2.85E-09;1.70E-09;5.00E-06;4.72E-03];
C = kinetics(theta0,t)
end
function C = kinetics(theta,t)
c0 = [3.16E-02;0.33];
fun = @(x) sum([-x(1)+(theta(7)*theta(6)*c0(1)*theta(3)*theta(4)/x(1))+2.0*(theta(8)*theta(7)*theta(6)*c0(1)*theta(3)*theta(4)/(x(1))^2)+(theta(5)/x(1)); -x(2)+((c0(2)*theta(7)*x(2))/(x(2)^2+theta(7)*x(2)+theta(7)*theta(8)))+2.0*((c0(2)*theta(7)*theta(8))/(x(2)^2+theta(7)*x(2)+theta(7)*theta(8)))+(theta(5)/x(2))].^2);
sol = fminsearch(fun,[1;1]);
fun(sol)
c0 = [c0;sol]
tspan = t;
% A constant, singular mass matrix
M = [1 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0];
% Options
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
% Solving
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
end
function dC = DifEq(t,c,theta)
dC = zeros(4,1);
dC(1) = (theta(1)/theta(2))*((1930.65/1000000)- (c(1)*theta(3)*theta(4) / 875000 )) - (c(1)* theta(3)* theta(4) - theta(6)* ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8))))/((1.0/theta(12)) + theta(6)/((1 + theta(9)* ((theta(7)* theta(6)* c(1)* theta(3)* theta(4) / c(3)) - ((c(2)* theta(7)* c(4))/(c(4)^2 + theta(7)* c(4) + theta(7) *theta(8))))/2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))) + theta(11)* ((theta(8)* theta(7)* theta(6)* c(1)* theta(3)* theta(4) / (c(3))^2) - ((c(2)* theta(7)* theta(8))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))))/ 2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))))* theta(13)));
dC(2) = (c(1)* theta(3)* theta(4) - theta(6)* ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8))))/((1.0/theta(12)) + theta(6)/((1 + theta(9)* ((theta(7)* theta(6)* c(1)* theta(3)* theta(4) / c(3)) - ((c(2)* theta(7)* c(4))/(c(4).^2 + theta(7) * c(4) + theta(7) *theta(8))))/2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))) + theta(11)* ((theta(8)* theta(7)* theta(6)* c(1)* theta(3)* theta(4) / (c(3))^2) - ((c(2)* theta(7)* theta(8))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))))/ 2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))))* theta(13)));
dC(3) = - c(3) + (theta(7)* theta(6)* c(1)* theta(3)* theta(4) / c(3)) + 2.0* (theta(8)* theta(7)* theta(6)* c(1)* theta(3)* theta(4) / (c(3))^2)+ (theta(5)/c(3));
dC(4) = - c(4) + ((c(2)* theta(7)* c(4))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))) + 2.0*((c(2)* theta(7)* theta(8))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))) + (theta(5)/c(4));
end

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Nonlinear Least Squares (Curve Fitting) 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by