The entries in tspan must strictly increase or decrease.

Hi everyone,

I am following:

https://www.mathworks.com/matlabcentral/answers/314965-parameter-estimation-for-a-system-of-differential-equations#comment_411192

in fitting experimental data to a system of coupled differential equations. Yet I get an error message:

"
Error using odearguments (line 87)
The entries in tspan must strictly increase or decrease.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in ODE4520Aug2018a/IntegratedModel (line 8)
[T,Cv]=ode45(@ODEloop,t,c0);

Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in ODE4520Aug2018a (line 89)
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,c,t);

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue."

Please see the m-file on the attachment.

function ODE4520Aug2018a
%--------------------------------------------------------------------------
 function C=IntegratedModel(k,t)

c0 = [0.01721262399 0.00000008759   67.022108400    0.00000487921   0.00000487921   49.17075376105  0.00000000000];

[T,Cv]=ode45(@ODEloop,t,c0);

function dC=ODEloop(t,c)

%%PARAMETERS  

A  = 1.5e-6;        
B   = 1.66667e-5;    
CC  = 6.51332e-2;    
D   = 0;             
E   = 8.314;         
F   = 323.15;       
G   = 149;           
H   = 4.14e-6;       
I   = 1.39e-9;       
J   = 2.89e-9;       
K   = 8.4e-4;        
L   = 9.598e-4;      
M   = 5.15e+3;       
N   = 3.53e-9;      
O   = 1.07e-7;      
P   = 10;           
Q   = 8.825e-3;     
R   = 12.54;        
S   = 100.0869;     
%k(1)= 0.84;         
TT   = 2703;         
U    = 1.7e-3;       
V    = 6.55e-8;      
W    = 6.24;         
X    =5.68e-5;       
Y    = 258.30;       
Z    = 2540;         
AA   = 0.00000933254;

dcdt=zeros(7,1);

dcdt(1) = (1/A) * (B * CC - B * c(1)) - ((c(1) * E * F - G * (((c(3) *  AA.^2) / (AA.^2  +  W * AA +  W * X))))/((1/H) + (G/((1 + (I* c(5))/(J* c(3))) * K))));

dcdt(2) = (1/A) * (B * D - B * c(2)) - (L * (1 + (I* c(5))/(N * c(4)))) * (c(2)* E * F/M - (((c(3) *  AA.^2) / (AA.^2  +  U * AA +  U * V))));

dcdt(3) = ((c(1) * E * F - G * ((c(3) *  AA.^2) / (AA.^2  +  W * AA +  W * X)))/((1/H) + (G/((1 + (I* c(5))/(J* c(3)))*K)) - (0.162 * exp(-5153 / F) * ((( c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA +  W * X) )) / O) -  1).^3 * (P / (( c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA +  W * X) )) / O)))));

dcdt(4) = (L * (1 + (I* c(5)) / (N * c(4))) * (((c(2)* E * F) / M)) - ((c(3) *  AA.^2) / (AA.^2  +  U * AA +  U * V))) - (-Q * R * S * c(6) * c(7) * (1 - (k(1) * c(7)) / (1 + k(1) * c(7))));

dcdt(5) = (-Q * R * S * c(6) * c(7) *(1 -(k(1) * c(7)) / (1 + k(1) * c(7)))) - (0.162 * exp(-5153/F) * (((c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA +  W * X) )) / O) - 1).^3 * (P / ((c(5) * ( (c(3) * W * X) / (AA.^2 + W * AA +  W * X) )) / O)));

dcdt(6) = -c(6) * (-Q * R * S * c(6) * c(7) * (1 - (k(1) * c(7)) / ( 1 + k(1) * c(7)))) * S / TT;

dcdt(7) = (-Q * R * S * c(6) * c(7) * (1- (k(1) * c(7))/(1 + k(1) * c(7))))* (Y/ Z);

dC=dcdt;
    end
C=Cv;
end

t=[1200
2400
10200
];

c=[
0.01721262399   0.00000008759   67.022108400    0.00000487921   0.00000487921   49.17075376105  0.00000000000
0.01700907268   0.00000010281   138.644926345   0.00000511161   0.00000511161   49.60948155721  0.00000000000
0.01741617529   0.00000036495   1002.322509102  0.00000535393   0.00000535393   30.91118867616  9.33482826985
]; 

k0 = [0.84];

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,c,t);

tv = linspace(min(t), max(t));
Cfit = IntegratedModel(k, tv);

figure(1)
plot(t, c, 'v')
hold on
hlp = plot(tv, Cfit);
hold off
xlabel('Time (sec)')
ylabel('Concentration (mol/m^3)')
legend(hlp, 'CSO_2,Headspace', 'C_{CO_2,Headspace}','C_(S_total}','C_{C_total}','C_{Ca^{2+}_total}','C_{CaCO_3}', 'Ca^2+', 'C_{H^{+}}', 'Location','E')

end

Answer 1

[Walter Roberson]

You have

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,c,t);

k0 is being passed in the x0 position -- an initial guess at the parameters to be found.

c is being passed in the xdata position.

t is being passed in the ydata position.

lsqcurvefit will then call IntegratedModel passing in "a vector or matrix x, and a vector or matrix xdata". Here, x would refer to the trial parameter whose effect is to be calculated at the positions given by xdata.

You are receiving the values in IntegratedModel under the names k and t. k is consistent with your naming use of passing in k0 and expecting k to be output. However, the second position is having xdata passed in, and your xdata is your c matrix, a 3 x 7 matrix which is not in numeric ascending or descending order.

In IntegratedModel you attempt to pass this "t" in to ode45 as if it were tspan values... but it isn't.

Perhaps your code should have been

    [k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,t,c);

[Dursman Mchabe]

Thanks a lot Walter, changing the positions of c and t, solves the problem. However, it creates an error of function value and YDATA sizes not equal.

" Error using lsqcurvefit (line 262) Function value and YDATA sizes are not equal.

Error in ODE4520Aug2018a (line 86) [k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,t,c);"

[Walter Roberson]

Well, you get

>> ODE4520Aug2018a
Warning: Failure at t=1.200008e+03.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t. 
> In ode45 (line 360)
  In ODE4520Aug2018a/IntegratedModel (line 8)
  In lsqcurvefit (line 213)
  In ODE4520Aug2018a (line 86)

so you are encountering a singularity in the ODE function that is causing ode45() to not return the size of array it is expected to.

[Dursman Mchabe]

Thanks a lot Walter. I guess I have to read more, currently I have no idea what should be my next move.

[Walter Roberson]

I suggest adding an appropriate options structure to the ode45() call, that has plotting configured into it, so that you can get an idea of which output is "exploding".

I notice that some of your dcdt entries divide by linear factors that are based upon the input c, which suggests to me that you might be getting a division by 0.

[Dursman Mchabe]

Thanks once more for pointing me in the right direction. I have included options

options = odeset('RelTol',1e-3,'AbsTol',1e-3,'Stats','on','OutputFcn',@odeplot);

The figure I get is on the attachment. The current error that I am working on now is

"

Error using snls (line 47) Objective function is returning undefined values at initial point. lsqcurvefit cannot continue.

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 ODE4520Aug2018a (line 87) [k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@IntegratedModel,k0,t,c);"

[Walter Roberson]

The figure did not make it.

[Walter Roberson]

Are you still getting the integration failures?

[Dursman Mchabe]

Objective function is returning undefined values at initial point. lsqcurvefit cannot continue. I'll try to attach the figure and a table of those 'undefined' values.

[Walter Roberson]

I would have to test to be sure but it appears to me that your code would have a divide by 0 if c(5) were 0, even only because the ode routine used 0 as a probe value while estimating gradients.

[Dursman Mchabe]

Thanks a lot Walter. In the mean time I will google 'how to check if a matlab code have a divide by 0'.

[Walter Roberson]

dbstop if naninf

[Dursman Mchabe]

Thank you. I try that.

[Dursman Mchabe]

It does not pick it. Maybe there is something else.

[Walter Roberson]

Some of your values get large. Then in dcdt(5), the expression being cubed becomes about -5E130, and when that is cubed that gives infinity. That infinity poisons the rest of your calculations.

[Dursman Mchabe]

Thanks a lot. I will start by squaring it. If I still get an error, then, I will remove the power.

[Dursman Mchabe]

Wow! Wow! Wow! Thank you so so so so much. I truly appreciate your help. It worked wonders. I thank you.

[Walter Roberson]

If you have the symbolic toolbox, I would recommend reading the first example for odeFunction to see how a system of ode can be written symbolically and converted for use with a numeric ode routine. That reduces the possibilities of making typing mistakes in the implementation.

[Dursman Mchabe]

I have a student access to Symbolic Math Toolbox Version 8.1 (R2018a). I will study the example straight-away.