Need a better guess y0 for consistent initial conditions.

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?

[Walter Roberson]

See https://groups.google.com/forum/#!topic/comp.soft-sys.matlab/F558nu1baPE

[Dursman Mchabe]

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]

ODE15s complains about a better guess for c0, not for theta.

[Dursman Mchabe]

Ooops. Sorry.

[Dursman Mchabe]

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]

Why do the first two values from c0 change ? They are your initial concentrations and have to be 0, I guess.

[Dursman Mchabe]

I am not sure why did fsolve give different values compared to experimental values.

[Torsten]

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]

Thankyou. I try that.

[Dursman Mchabe]

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]

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]

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]

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.

[Dursman Mchabe]

First call.

[Torsten]

Please print "dC" in "DifEq". What do you get ?

[Dursman Mchabe]

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]

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]

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]

Remove the ";" after the definition

dC = [(theta(1)./...];

[Dursman Mchabe]

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]

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]

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]

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]

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]

Thanks a lot, I try that.

[Dursman Mchabe]

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]

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]

Sorry, I am a bit confused. Does "old definition" in this case mean, I should replace "fun" by "DifEq"?

[Torsten]

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]

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.

[Torsten]

For me, it works.

[Dursman Mchabe]

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]

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]

I am using version '9.5.0.944444 (R2018b)'

[Torsten]

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]

It works!!!! Thank you! Thank You! Thank you!