matlab for loop within ode solver

Question

i am using ODE solver for solving my system of ODE. Now i also want to calculate sensitivity of various parameters by varying them. My basic ODE solver program is something like: Function:

function F = ODE_site(t,y)
%%Parameter block
k1f=1e8;
k1b=1e5; %%and so on

%%Elementary rate block
r1=k1f*y(1)-k1b*P*y(5);
and so on

%%Mass balance block
F=[  r1-r5-r78+r86 ;
    and so on ];
%% End of function

The main ODE solver looks like:

optode = odeset('NonNegative',1:41,'Abstol',1E-20,'RelTol',1E-20);
y0=zeros(41,1);
y0(41) = 1;
[t,y]=ode15s(@ODE_site,[0,50000000000],y0,optode);
%% Plus some other lines to plot the solution
%% End of ODE solver

Now, I want to introduce a for loop in the main solver but the problem is if i use for loop in the main solver, it cannot pass the command to the function. I tried to use global to assign the variable to all over the environment but it didn't work. How can i assign the for loop for both the ODE solver code and function? Any help will be highly appreciated. Thanks, Mamun

BerndGit BerndGit · Accepted Answer · 2015-03-23T20:02:38

Define the parameters, which you want to modify as addional input parameter to the function:

function F = ODE_site(t,y,param)
%%Parameter block
k1f=param(1);
k1b=param(2); %%and so on

%%Elementary rate block
r1=k1f*y(1)-k1b*P*y(5);
and so on

%%Mass balance block
F=[  r1-r5-r78+r86 ;
    and so on ];
%% End of function

You can hand over the parameters to your function in ode15() after the optimisation Parameters:

optode = odeset('NonNegative',1:41,'Abstol',1E-20,'RelTol',1E-20);
y0=zeros(41,1);
y0(41) = 1;
para=[1e5, 1e8]
[t,y]=ode15s(@ODE_site,[0,50000000000],y0,optode,para);

Two notes:

*) Code is untested. Since your example was not executable I didn't care to make a running example too.

*) It looks like you are solving chemical reaction kinetics. Most likely you could further optimize your code by using matrix and vector operations. So instead of usning separate variables r1, r2, r3, .. you could store it in one vector r(1), r(2), r(3), ... Your last line, than could be written by the product of the reaction vector with the stoechiometric matrix.

matlab for loop within ode solver

1 Answers