ode23s, ode15s, what to use?

5 views (last 30 days)
Laura
Laura on 4 Oct 2012
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
  1 Comment
Babak
Babak on 4 Oct 2012
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?

Sign in to comment.

Answers (2)

Laura
Laura on 4 Oct 2012
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.

Babak
Babak on 4 Oct 2012
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!