Hi everyone,
I am following:
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; endC=Cv;endt=[1200240010200];c=[0.01721262399 0.00000008759 67.022108400 0.00000487921 0.00000487921 49.17075376105 0.000000000000.01700907268 0.00000010281 138.644926345 0.00000511161 0.00000511161 49.60948155721 0.000000000000.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 onhlp = plot(tv, Cfit);hold offxlabel('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
Best Answer