Need a better guess y0 for consistent initial conditions.

14 Jan 2019
0 Answers
Updated 16 Jan 2019
7 Views (30 days)

0 votes

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 Comments
Show 34 older comments Hide 34 older comments

Dursman Mchabe
Dursman Mchabe on 14 Jan 2019
Hi Walter, following the discussions on https://groups.google.com/forum/#!topic/comp.soft-sys.matlab/F558nu1baPE
I decided to use fsolve to evaluate the 'better guesses', which I did as follows:
%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];
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];
options = optimoptions('fsolve','Display','iter');
[theta,fval] = fsolve(@DifEq,theta0,options)
function dC = DifEq(theta)
dC = [(theta(1)./theta(2)).*((1930.65./1000000)- (3.16E-02 .*theta(3) .*theta(4) ./ 875000 )) - (3.16E-02.* theta(3).* theta(4) - theta(6).* ((0.33.* 1.35E-04.^2)./(1.35E-04.^2 + theta(7).* 1.35E-04 + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* 3.16E-02.* theta(3).* theta(4) ./ 1.62E-04) - ((0.33.* theta(7).* 1.35E-04)./(1.35E-04.^2 + theta(7) .* 1.35E-04 + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* 3.16E-02.* theta(3).* theta(4)) - ((0.33.* 1.35E-04.^2)./(1.35E-04.^2 + theta(7).* 1.35E-04 + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* 3.16E-02.* theta(3).* theta(4) ./ (1.62E-04).^2) - ((0.33.* theta(7).* theta(8))./(1.35E-04.^2 + theta(7) .* 1.35E-04 + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* 3.16E-02.* theta(3).* theta(4)) - ((0.33.* 1.35E-04.^2)./(1.35E-04.^2 + theta(7).* 1.35E-04 + theta(7).* theta(8))))).* theta(13)))
(3.16E-02.* theta(3).* theta(4) - theta(6).* ((0.33.* 1.35E-04.^2)./(1.35E-04.^2 + theta(7).* 1.35E-04 + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* 3.16E-02.* theta(3).* theta(4) ./ 1.62E-04) - ((0.33.* theta(7).* 1.35E-04)./(1.35E-04.^2 + theta(7) .* 1.35E-04 + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* 3.16E-02.* theta(3).* theta(4)) - ((0.33.* 1.35E-04.^2)./(1.35E-04.^2 + theta(7).* 1.35E-04 + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* 3.16E-02.* theta(3).* theta(4) ./ (1.62E-04).^2) - ((0.33.* theta(7).* theta(8))./(1.35E-04.^2 + theta(7) .* 1.35E-04 + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* 3.16E-02.* theta(3).* theta(4)) - ((0.33.* 1.35E-04.^2)./(1.35E-04.^2 + theta(7).* 1.35E-04 + theta(7).* theta(8))))).* theta(13)))
- 1.62E-04 + (theta(7).* theta(6).* 3.16E-02.* theta(3).* theta(4) ./ 1.62E-04) + 2.* (theta(8).* theta(7).* theta(6).* 3.16E-02.* theta(3).* theta(4) ./ (1.62E-04).^2)+ (theta(5)./1.62E-04)
- 1.35E-04 + ((0.33.* theta(7).* 1.35E-04)./(1.35E-04.^2 + theta(7) .* 1.35E-04 + theta(7) .*theta(8))) + 2.*((0.33.* theta(7).* theta(8))./(1.35E-04.^2 + theta(7) .* 1.35E-04 + theta(7) .*theta(8))) + (theta(5)./1.35E-04)];
end
The above code gives the following inital guesses for theta:
theta0 = [2.1782e-07;2.1153e-03;8.3140e+00;3.2315e+02;2.2292e-05;1.4340e+02;6.2400e+00;-8.1000e-05;2.0400e-09;2.8500e-09;1.7000e-09;8.4011e-10; 4.7200e-03];
However, when I run with theta0, I again get the same error message of better guess:
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/objective (line 278)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});
Error in snls (line 338)
newfvec = feval(funfcn{3},xcurr,varargin{:});
Error in lsqncommon (line 167)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
Error in lsqcurvefit (line 271)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...
Error in IgorWaterODE15s (line 169)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Torsten
Torsten on 14 Jan 2019
ODE15s complains about a better guess for c0, not for theta.
Dursman Mchabe
Dursman Mchabe on 14 Jan 2019
Thanks a lot, Torsten. I found
c0 = [4.3061e-03;1.2720e+01;1.0175e+02;5.6965e+00];
which is working, and leads to the new error:
Matrix dimensions must agree.
Error in lsqcurvefit/objective (line 279)
F = F - YDATA;
Error in snls (line 338)
newfvec = feval(funfcn{3},xcurr,varargin{:});
Error in lsqncommon (line 167)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
Error in lsqcurvefit (line 271)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...
Error in IgorWaterODE15s (line 171)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Torsten
Torsten on 14 Jan 2019
Why do the first two values from c0 change ? They are your initial concentrations and have to be 0, I guess.
Dursman Mchabe
Dursman Mchabe on 14 Jan 2019
I am not sure why did fsolve give different values compared to experimental values.
Torsten
Torsten on 14 Jan 2019
You must fix the first two values c(1) and c(2) to the experimental values and only adjust c(3) and c(4) such that the last two algebraic equations are satisfied.
Dursman Mchabe
Dursman Mchabe on 14 Jan 2019
I have tried that using the fsolve code below:
c0 = [1.62E-04;1.35E-04];
options = optimoptions('fsolve','Display','iter');
[c,fval] = fsolve(@DifEq,c0,options)
function dC = DifEq(c)
dC = [(5.97E-04./1.95E-03).*((1930.65./1000000)- (3.16E-02 .*8.314 .*323.15 ./ 875000 )) - (3.16E-02.* 8.314.* 323.15 - 143.40.* ((0.33.* c(2).^2)./(c(2).^2 + 6.24.* c(2) + 6.24.* 5.68E-05)))./((1./5.00E-06) + 143.40./((1 + 2.04E-09.* ((6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ c(1)) - ((0.33.* 6.24.* c(2))./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)))./2.85E-09.* ((143.40.* 3.16E-02.* 8.314.* 323.15) - ((0.33.* c(2).^2)./(c(2).^2 + 6.24.* c(2) + 6.24.* 5.68E-05))) + 1.70E-09.* ((5.68E-05.* 6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ (c(1)).^2) - ((0.33.* 6.24.* 5.68E-05)./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)))./ 2.85E-09.* ((143.40.* 3.16E-02.* 8.314.* 323.15) - ((0.33.* c(2).^2)./(c(2).^2 + 6.24.* c(2) + 6.24.* 5.68E-05)))).* 4.72E-03))
(3.16E-02.* 8.314.* 323.15 - 143.40.* ((0.33.* c(2).^2)./(c(2).^2 + 6.24.* c(2) + 6.24.* 5.68E-05)))./((1./5.00E-06) + 143.40./((1 + 2.04E-09.* ((6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ c(1)) - ((0.33.* 6.24.* c(2))./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)))./2.85E-09.* ((143.40.* 3.16E-02.* 8.314.* 323.15) - ((0.33.* c(2).^2)./(c(2).^2 + 6.24.* c(2) + 6.24.* 5.68E-05))) + 1.70E-09.* ((5.68E-05.* 6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ (c(1)).^2) - ((0.33.* 6.24.* 5.68E-05)./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)))./ 2.85E-09.* ((143.40.* 3.16E-02.* 8.314.* 323.15) - ((0.33.* c(2).^2)./(c(2).^2 + 6.24.* c(2) + 6.24.* 5.68E-05)))).* 4.72E-03))
- c(1) + (6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ c(1)) + 2.* (5.68E-05.* 6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ (c(1)).^2)+ (5.3E-8./c(1))
- c(2) + ((0.33.* 6.24.* c(2))./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)) + 2.*((0.33.* 6.24.* 5.68E-05)./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)) + (5.3E-8./c(2))];
end
Which gave me:
c =
4.1324e-01
3.1424e-01
And on using that in the ODE15s code, I get this 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 60)
[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 171)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
Torsten
Torsten on 15 Jan 2019
According to your data, you have to work with c0(1)=3.16e-2 and c0(2)=0.33. Since the first two equations are differential equations, these values have to remain fixed. Now insert these values in the third and the fourth algebraic equations and solve them for c0(3) and c0(4). Then put [c0(1),c0(2),c0(3),c0(4)] together in one vector and use it as initial guess for ODE15S. This procedure should be done each time "kinetics" is called. Of course, the guesses for c0(3) and c0(4)will vary because "theta" varies with each iteration.
Dursman Mchabe
Dursman Mchabe on 15 Jan 2019
Thanks a lot Torsten, I used
c0 = [1;1];
options = optimoptions('fsolve','Display','iter');
[c,fval] = fsolve(@DifEq,c0,options)
function dC = DifEq(c)
dC = [ - c(1) + (6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ c(1)) + 2.* (5.68E-05.* 6.24.* 143.40.* 3.16E-02.* 8.314.* 323.15 ./ (c(1)).^2)+ (5.3E-8./c(1))
- c(2) + ((0.33.* 6.24.* c(2))./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)) + 2.*((0.33.* 6.24.* 5.68E-05)./(c(2).^2 + 6.24 .* c(2) + 6.24 .*5.68E-05)) + (5.3E-8./c(2))];
end
to evaluate c0(3) and c0(4) for the ode15s, then I got:
c =
2.7562e+02
3.1424e-01
Then I ran the ode15s as
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;2.7562e+02;3.1424e-01];
%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];
%theta0 = [2.1782e-07;2.1153e-03;8.3140e+00;3.2315e+02;2.2292e-05;1.4340e+02;6.2400e+00;-8.1000e-05;2.0400e-09;2.8500e-09;1.7000e-09;8.4011e-10; 4.7200e-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
Strangely, I still get 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 60)
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
Error in lsqcurvefit/objective (line 278)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});
Error in snls (line 338)
newfvec = feval(funfcn{3},xcurr,varargin{:});
Error in lsqncommon (line 167)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
Error in lsqcurvefit (line 271)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...
Error in IgorWaterODE15s (line 171)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Torsten
Torsten on 15 Jan 2019
Does ODE15S fail in the first call to "kinetics" or later ? As I told you, the initial guess for C0 has to be modified each time "kinetics" is called because the theta vector changes.
Torsten
Torsten on 15 Jan 2019
Please print "dC" in "DifEq". What do you get ?
Dursman Mchabe
Dursman Mchabe on 15 Jan 2019
When I put
fprintf('t= %f dC= %f\n',mat2str(t),mat2str(dC))
, in
DifEq
I get :
t= 50.000000 dC= 52.000000
t= 53.000000 dC= 53.000000
t= 57.000000 dC= 54.000000
t= 54.000000 dC= 54.000000
t= 55.000000 dC= 52.000000
t= 91.000000 dC= 49.000000
t= 46.000000 dC= 48.000000
t= 57.000000 dC= 54.000000
t= 53.000000 dC= 54.000000
t= 57.000000 dC= 48.000000
t= 57.000000 dC= 54.000000
t= 50.000000 dC= 49.000000
t= 52.000000 dC= 57.000000
t= 55.000000 dC= 101.000000
t= 45.000000 dC= 48.000000
t= 53.000000 dC= 59.000000
t= 48.000000 dC= 46.000000
t= 48.000000 dC= 48.000000
t= 48.000000 dC= 52.000000
t= 57.000000 dC= 53.000000
.
.
.
This goes on for a longer period
Torsten
Torsten on 15 Jan 2019
Why do you get one value for dC ? dC should be a vector with four elements.
Just add the line
dC
at the end of "DifEq"
What is written to console ?
Dursman Mchabe
Dursman Mchabe on 15 Jan 2019
dC will be undefined outside the function (after end), I will get:
Undefined function or variable 'dC'.
Error in IgorWaterODE15s/kinetics (line 85)
fprintf('t= %f dC= %f\n',mat2str(t),mat2str(dC))
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in IgorWaterODE15s (line 173)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
I expected dC to give a vector of 4 elements.
Torsten
Torsten on 15 Jan 2019
Remove the ";" after the definition
dC = [(theta(1)./...];
Dursman Mchabe
Dursman Mchabe on 15 Jan 2019
I now get:
t= 52.000000 dC= 55.000000
t= 51.000000 dC= 48.000000
t= 57.000000 dC= 53.000000
t= 56.000000 dC= 51.000000
t= 55.000000 dC= 56.000000
t= 52.000000 dC= 52.000000
t= 55.000000 dC= 56.000000
t= 52.000000 dC= 59.000000
t= 56.000000 dC= 46.000000
t= 50.000000 dC= 50.000000
t= 56.000000 dC= 55.000000
t= 50.000000 dC= 50.000000
t= 55.000000 dC= 51.000000
t= 57.000000 dC= 57.000000
t= 52.000000 dC= 53.000000
t= 57.000000 dC= 52.000000
t= 101.000000 dC= 45.000000
t= 48.000000 dC= 53.000000
t= 93.000000 dC=
dC =
1.1273e-05
4.3613e-04
-2.7260e-01
8.3383e-04
t= 49.000000 dC= 49.000000
t= 52.000000 dC= 52.000000
t= 50.000000 dC= 46.000000
t= 54.000000 dC= 55.000000
t= 48.000000 dC= 51.000000
t= 49.000000 dC= 57.000000
t= 57.000000 dC= 56.000000
t= 57.000000 dC= 54.000000
t= 91.000000 dC= 49.000000
t= 46.000000 dC= 49.000000
t= 50.000000 dC= 55.000000
t= 50.000000 dC= 53.000000
t= 57.000000 dC= 53.000000
t= 53.000000 dC= 52.000000
t= 51.000000 dC= 49.000000
t= 52.000000 dC= 49.000000
t= 55.000000 dC= 101.000000
t= 45.000000 dC= 48.000000
t= 53.000000 dC= 59.000000
t= 48.000000 dC= 46.000000
t= 48.000000 dC= 48.000000
t= 48.000000 dC= 52.000000
t= 51.000000 dC= 54.000000
t= 49.000000 dC= 50.000000
t= 57.000000 dC= 49.000000
t= 56.000000 dC= 53.000000
t= 54.000000 dC= 52.000000
t= 48.000000 dC= 49.000000
t= 49.000000 dC= 52.000000
t= 59.000000 dC= 45.000000
t= 48.000000 dC= 46.000000
t= 50.000000 dC= 55.000000
t= 50.000000 dC= 53.000000
t= 57.000000 dC= 56.000000
t= 51.000000 dC= 54.000000
t= 53.000000 dC= 48.000000
t= 48.000000 dC= 57.000000
t= 54.000000 dC= 53.000000
t= 49.000000 dC= 59.000000
t= 48.000000 dC= 46.000000
t= 48.000000 dC= 48.000000
t= 48.000000 dC= 56.000000
t= 51.000000 dC= 51.000000
t= 56.000000 dC= 51.000000
t= 48.000000 dC= 52.000000
t= 57.000000 dC= 57.000000
t= 56.000000 dC= 53.000000
t= 51.000000 dC= 55.000000
t= 50.000000 dC= 49.000000
t= 93.000000 dC=
dC =
1.2018e-05
4.3542e-04
1.8401e-01
-2.3003e-04
.
.
.
repeating many times
Torsten
Torsten on 15 Jan 2019
Edited: Torsten on 15 Jan 2019
Try
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];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
end
function C = kinetics(theta,t)
c0 = [3.16E-02;0.33];
fun=@(x) [-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))];
sol = fsolve(fun,[1,1]);
c0 = [c0;sol]
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,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
Dursman Mchabe
Dursman Mchabe on 15 Jan 2019
Thanks a lot Torsten, the code gives the following error:
Not enough input arguments.
Error in Torsten (line 3)
fun=@(x) [-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))];
>>
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
I truly appreciate that you went an extra mile in helping. I have also tried the last suggestion. It gave this error message:
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Error in main1>kinetics (line 43)
c0 = [c0;sol]
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in main1 (line 37)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
Torsten
Torsten on 16 Jan 2019
With this code for "kinetics", the first iteration step worked for me.
I don't have the optimization toolbox licenced - so you may want to replace the call to "fminsearch" by a call to "fsolve".
Note that I deleted the "Vectorized"-option in the call to "ode15s".
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]);
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
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
When I try fsolve:
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 = fsolve(fun,[1;1]);
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
The old error comes back.
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 main1>kinetics (line 50)
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
Error in lsqcurvefit/objective (line 278)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});
Error in snls (line 338)
newfvec = feval(funfcn{3},xcurr,varargin{:});
Error in lsqncommon (line 167)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
Error in lsqcurvefit (line 271)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...
Error in main1 (line 37)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Can I please ask for the code that worked for you, so that I can try it as well.
Torsten
Torsten on 16 Jan 2019
Edited: Torsten on 16 Jan 2019
You have to replace "fun" with the old definition.
Try
function C = kinetics(theta,t)
c0 = [3.16E-02;0.33];
fun=@(x) [-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))];
sol = fsolve(fun,[1;1]);
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
or take the code from above using "fminsearch".
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
Sorry, I am a bit confused. Does "old definition" in this case mean, I should replace "fun" by "DifEq"?
Torsten
Torsten on 16 Jan 2019
Try
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];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
end
function C = kinetics(theta,t)
c0 = [3.16E-02;0.33];
fun=@(x) [-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))];
sol = fsolve(fun,[1;1]);
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
If it doesn't work, replace function "kinetics" by
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]);
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
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
That lead to too many input arguments, resulting in an error message:
Error using
main1>@(x)[-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))]
Too many input arguments.
Error in main1>@(t,y)DifEq(t,y,theta) (line 50)
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 150)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in main1>kinetics (line 50)
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in main1 (line 37)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
It seems there is a very small mistake hidden. When I the first suggestion I get an 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 main2>kinetics (line 50)
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
Error in lsqcurvefit/objective (line 278)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});
Error in snls (line 338)
newfvec = feval(funfcn{3},xcurr,varargin{:});
Error in lsqncommon (line 167)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
Error in lsqcurvefit (line 271)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...
Error in main2 (line 37)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
And when I replace "kinetics" with:
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]);
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
I get an 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 main2>kinetics (line 50)
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
Error in lsqcurvefit/objective (line 278)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});
Error in snls (line 338)
newfvec = feval(funfcn{3},xcurr,varargin{:});
Error in lsqncommon (line 167)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
Error in lsqcurvefit (line 271)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...
Error in main2 (line 37)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
Thanks a lot for the help. I truly appreciate. I will try to figure out my mistake. I will comment when I win.
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
I am using version '9.5.0.944444 (R2018b)'
Torsten
Torsten on 16 Jan 2019
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
Dursman Mchabe
Dursman Mchabe on 16 Jan 2019
It works!!!! Thank you! Thank You! Thank you!

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Tags

Asked:

Dursman Mchabe
Dursman Mchabe
on 14 Jan 2019

Commented:

Dursman Mchabe
Dursman Mchabe
on 16 Jan 2019

Community Treasure Hunt

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

Start Hunting!