MATLAB Answers

Need a better guess y0 for consistent initial conditions.

Dursman Mchabe 님이 질문을 제출함. 14 Jan 2019
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

표시 이전 댓글 수: 33
I am using version '9.5.0.944444 (R2018b)'
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
It works!!!! Thank you! Thank You! Thank you!

로그인 to comment.

태그

답변 수: 0


Translated by