Thread Subject: error message on ode15s?

Hi all,

I'm trying to use ode15s to solve 6 ODE's, (since ode45 failed) however I'm getting the error message:

Warning: Failure at t=0.000000e+000. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (0.000000e+000) at time t.

I have checked for divisions by 0, and not sure why else I would be getting this error message, can anyone please explain?

Thanks in advance!



% Model parameters -----------------------------------------------------
      m = 1;
      chi_L = 1.28;
      chi_0 = 0;
      psi_L = 0.279;
      omega_L = 1.56;
      sigma = 2.74;
      gamma = 0.000294;
      gamma_rr = 30;
      gamma_s = 0.54;
      lambda = 1.56;
      Rchol_L = 3400;
      K = 2;
      f = 0.9;
      W = 1500;

% Initial values -------------------------------------------------------
        y0(1)=1; %lE(0) = l_0 = 1
        y0(2)=0; %lB(0) = 0
        y0(3)=0; %lI(0) = 0
        y0(4)=1; %pF(0) = p0 = 1
        y0(5)=0 %pI(0) = 0
        y0(6)=0.7 % c(0) = 0.7

% ODE solver ------------------------------------------------------------
% [t, y]=ode15s('name of equation', [t0:step:tf], y0);
        [t y] = ode15s('ndLDL', [0:0.01:10000], y0);

% Plot results ----------------------------------------------------------
        lE=y(:,1);
        lB=y(:,2);
        lI=y(:,3);
        pF=y(:,4);
        pI=y(:,5);
         c=y(:,6);

plot(t,lE,'r',t,lB,'b',t,lI,'g',t,pF,'c',t,pI,'y',t,c,'k');
legend ('[lE]', '[lB]', '[lI]', '[pF]', '[pI]', '[c]');
xlabel('Time')
ylabel('Concentration')
title('Plot of Endocysotis Process over Time')

function dydt = ndLDL(t, y)

% Rate constants -------------------------------------------------------
      m = 1;
      chi_L = 1.28;
      chi_0 = 0;
      psi_L = 0.279;
      omega_L = 1.56;
      sigma = 2.74;
      gamma = 0.000294;
      gamma_rr = 30;
      gamma_s = 0.54;
      lambda = 1.56;
      Rchol_L = 3400;
      K = 2;
      f = 0.9;
      W = 1500;

% Work in a more comprehensible set of variables -----------------------
      
      lE = y(1);
      lB = y(2);
      lI = y(3);
      pF = y(4);
      pI = y(5);
      c = y(6);

% Dynamic equations-----------------------------------------------------

  dydt(1) = (-lE*pF + psi_L*lB)/W % dlE/dt
  dydt(2) = pF*lE - psi_L*lB - chi_L*lB % dlB/dt
  dydt(3) = chi_L*lB - omega_L*lI % dlI/dt
  dydt(4) = (gamma_rr*sigma*pI - chi_0*pF - (m*lE*pF - m*psi_L*lB + (m*chi_L*lB*pF)/(1 - pF)))/sigma % dpF/dt
  dydt(5) = ((gamma_s*sigma)/(K + c) + chi_0*f*pF + f*(1 +(pF/1 - pF))*(m*chi_L*lB) - sigma*gamma_rr*pI)/sigma % dpI/dt
  dydt(6) = gamma*(omega_L*Rchol_L*lI) - lambda*(c - 1) % dc/dt
  
  dydt = dydt'; %MATLAB v6 wants a column vector

> Hi all,
>
> I'm trying to use ode15s to solve 6 ODE's, (since
> ode45 failed) however I'm getting the error message:
>
> Warning: Failure at t=0.000000e+000. Unable to meet
> integration tolerances without reducing the step size
> below the smallest value allowed (0.000000e+000) at
> time t.
>
> I have checked for divisions by 0, and not sure why
> else I would be getting this error message, can
> anyone please explain?
>

Yes, division by zero is most probably the reason
for aborting the integration.

> Thanks in advance!
>
>
>
> % Model parameters
> -----------------------------------------------------
> m = 1;
> chi_L = 1.28;
> chi_0 = 0;
> psi_L = 0.279;
> omega_L = 1.56;
> sigma = 2.74;
> gamma = 0.000294;
> gamma_rr = 30;
> gamma_s = 0.54;
> lambda = 1.56;
> Rchol_L = 3400;
> K = 2;
> f = 0.9;
> W = 1500;
>
> % Initial values
> ------------------------------------------------------
> -
> y0(1)=1; %lE(0) = l_0 = 1
> y0(2)=0; %lB(0) = 0
> y0(3)=0; %lI(0) = 0
> y0(4)=1; %pF(0) = p0 = 1
> y0(5)=0 %pI(0) = 0
> y0(6)=0.7 % c(0) = 0.7
>
> % ODE solver
> ------------------------------------------------------
> ------
> % [t, y]=ode15s('name of equation', [t0:step:tf],
> y0);
> [t y] = ode15s('ndLDL', [0:0.01:10000], y0);
>
> % Plot results
> ------------------------------------------------------
> ----
> lE=y(:,1);
> lB=y(:,2);
> lI=y(:,3);
> pF=y(:,4);
> pI=y(:,5);
> c=y(:,6);
>
> plot(t,lE,'r',t,lB,'b',t,lI,'g',t,pF,'c',t,pI,'y',t,c,
> 'k');
> legend ('[lE]', '[lB]', '[lI]', '[pF]', '[pI]',
> '[c]');
> xlabel('Time')
> ylabel('Concentration')
> title('Plot of Endocysotis Process over Time')
>
> function dydt = ndLDL(t, y)
>
> % Rate constants
> ------------------------------------------------------
> -
> m = 1;
> chi_L = 1.28;
> chi_0 = 0;
> psi_L = 0.279;
> omega_L = 1.56;
> sigma = 2.74;
> gamma = 0.000294;
> gamma_rr = 30;
> gamma_s = 0.54;
> lambda = 1.56;
> Rchol_L = 3400;
> K = 2;
> f = 0.9;
> W = 1500;
>
> % Work in a more comprehensible set of variables
> -----------------------
>
> lE = y(1);
> lB = y(2);
> lI = y(3);
> pF = y(4);
> pI = y(5);
> c = y(6);
>
> % Dynamic
> equations---------------------------------------------
> --------
>
> dydt(1) = (-lE*pF + psi_L*lB)/W %
> % dlE/dt
> dydt(2) = pF*lE - psi_L*lB - chi_L*lB %
> % dlB/dt
> dydt(3) = chi_L*lB - omega_L*lI %
> % dlI/dt
> dydt(4) = (gamma_rr*sigma*pI - chi_0*pF - (m*lE*pF
> pF - m*psi_L*lB + (m*chi_L*lB*pF)/(1 - pF)))/sigma
> % dpF/dt

You divide by 1 - pF here which is zero at t=0.

> dydt(5) = ((gamma_s*sigma)/(K + c) + chi_0*f*pF +
> + f*(1 +(pF/1 - pF))*(m*chi_L*lB) -

Maybe you mean
+ f*(1 +(pF/(1 - pF)))*(m*chi_L*lB) -
here which would also produce a division by zero at t=0.

> sigma*gamma_rr*pI)/sigma % dpI/dt
> dydt(6) = gamma*(omega_L*Rchol_L*lI) - lambda*(c -
> - 1) % dc/dt
>
> dydt = dydt'; %MATLAB v6 wants a
> a column vector

Best wishes
Torsten.