ode23s, ode15s, what to use?

Hi, I have a set of 5 ODE's with 12 unknown parameters (yes, I know this is a LOT of parameters!) that I am trying to fit to data. I know my problem is stiff but it is not clear to me which ode solver I should use, ode23s or ode15s.

I had some initial guess for the values of the parameters so I ran the ode with both ode23s and ode15s.

I was expecting ode23s and ode15s to behave similarly, but the differences are huge! Please see attached results.

I read the documentation but cannot figure out which one is better. So, my question is: How to decide between ode15s and ode23s?

this is my code: function xout = testStiffODE clear all clc %clf

%initialConditions:
y0c1 = [297*10^3,0.01,0,0];
y0c2 = [999*10^3,0,0,0];
y0exp1 = [870*10^3,185*10^3,0,0];
y0exp2 = [710*10^3,953*10^3,0,0];

y0 = [y0c1;y0c2;y0exp1;y0exp2];

tvalues = 7*[0.0  0.4  1.0  1.4  2.0  3.0  4.0  6.0  8.0 10.0 12.0 16.0 20.0 24.0 28.0 32.0];

%parameters
a = 1.0;
aprime = 2*10^(-1);
beta = 1*10^(-8);   
c = 4;
d = 10^(-3);
delta = 0.88;
k = 5*10^3;
kappa = 1*10^(-2);
lambdaS = 1*10^(-2);
lambdaR = 4*10^(-2);
m = 10^3;
n = 1*10^(4);
p = 10^(-5);
rho = 1    ;
theta = 0.0001;
V0 =  1.e-4;
p0eq32 = [aprime, beta, c, d, ...
          k, kappa, lambdaS,...
          lambdaR,m, n, theta,V0];

%extra parameters:
extraparams = [[297.0*10^3,0];[999.0*10^3,0];[870.0*10^3,185*10^3];[710.0*10^3,444*10^3];[775.72*10^3, 84.28*10^3]];
extraparams1 = [[297.0*10^3,0, delta];[999.0*10^3,0, delta];[870.0*10^3,185*10^3, delta];[710.0*10^3,444*10^3, delta];[775.72*10^3, 84.28*10^3,delta]];

%create test data:
myfunName = @myeq;
ival = 3;
initCond = [y0(ival,:), p0eq32(1,end)];
[tout, yout] = ode23s(myfunName,tvalues,initCond,[],p0eq32(1:end-1),extraparams1(ival,:));
[tout1, yout1] = ode15s(myfunName,tvalues,initCond,[],p0eq32(1:end-1),extraparams1(ival,:));

abs(yout1 - yout)

function dydt = myeq(t,y,params,extraparams)
myparams = num2cell(params);
[aprime, beta, c, d, k, kappa, lambdaS, lambdaR,m, n, theta] = myparams{:};

S0 = extraparams(1,1);
R0 = extraparams(1,2);
delta = extraparams(1,3);

S= y(1,:);  R = y(2,:); I = y(3,:); L = y(4,:); V = y(5,:);
%definition of g:
g = kappa*(I+L)/(I+L+m);

dydt0 = lambdaS*S0 -(d + aprime*(V./(V+n)))*S- beta*S*V;
dydt1 = lambdaR*R0 -(d + aprime*(V./(V+n)))*R -beta*theta*R*V;
dydt2 = beta*S*V - (delta + aprime*(V./(V+n)) + g)*I;
dydt3 = theta*beta*R*V - (delta + aprime*(V./(V+n)) + g)*L;
dydt4 = k*(I+L)  - c*V; 
dydt = [dydt0; dydt1; dydt2; dydt3; dydt4];

This is the output (I just printed the last column for the sake of length):

Column 5

                         0
          0.19043550432143
          433296.853832826
          105131.802523792
          10655.6375275331
          298.415963380943
           22.032777486922
          21.1542593133763
          15.3459447373316
          45.0865582953029
          14.4720186743434
          12.6140635681131
          23.0151598653374
          4.96503961268536
          13.0580526547419
          13.9938160755173

[Babak]

Seems you want the integrators integrate the equations in myfunName = @shiveq33;

but

Where did you define shiveq33?

I don't find it in your code above.

You are also showing myeq() which is not used here at all...

By the way, what are the 12 unknowns you are talking about?

Answer 1

[Laura]

Hi Babak,

Yes, sorry, I meant to say @myeq which is defined above. The 12 unknowns I am talking about are the parameters of the ODE, defined in the vector p0eq32. In my real problem, I don't know what these parameters are, and this is just an initial guess. Once I figure out what is the correct solver to use, I can ask the following question, which relates to fitting the ode to data that I have.

Answer 2

[Babak]

As for the discrepansy of ODE15s and ODE23s or ODE23tb 's responses, I would suggest you to integrate the equations for a very short period of time with a fine timespan. Sometimes the ODE solvers show different behaviour when they reacha singularity. This might have happenned in your case because I guess you probably just randomly set the coefficients in p0eq32.