Need a better guess y0 for consistent initial conditions.

Question

Dursman Mchabe 2019-1-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/439651-need-a-better-guess-y0-for-consistent-initial-conditions

评论： 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更早的评论隐藏 34更早的评论

Dursman Mchabe 2019-1-14

在 MATLAB Online 中打开

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);

Dursman Mchabe 2019-1-14

在 MATLAB Online 中打开

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.

Dursman Mchabe 2019-1-15

在 MATLAB Online 中打开

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);

Dursman Mchabe 2019-1-15

在 MATLAB Online 中打开

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 2019-1-15

编辑：Torsten 2019-1-15

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 2019-1-16

在 MATLAB Online 中打开

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 2019-1-16

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 2019-1-16

在 MATLAB Online 中打开

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 2019-1-16

在 MATLAB Online 中打开

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);

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