how to reset/change the initial data?

Hi there, I am working on a non linear system of 4 ordinary differential equations with time varying coefficients. I solved them numerically using ode 15s. It worked well. But now i want to stop my solver every 2 or 3 years and change/reset the initial data and run it again. I dont know how to do it. Any guidance/help is greatly appreciated.

[John Petersen]

How did you start it in the first place? Can't you just do that?

[Rose]

i gave certain initial values to the dependent variables in he beginning and ran the code for almost 20 years. but now i want to change/ rest the initial values every 2 or 3 years. i am not sure how to do this.

Answer 1

[Star Strider]

You can give the ODE solvers a vector of time points at which to evaluate your ODE function. I've done this often, especially because it constrains the integration time without significantly affecting the accuracy of the integration.

From the documentation for the ODE solvers http://www.mathworks.com/help/techdoc/ref/ode23.html :

tspan — A vector specifying the interval of integration, [t0,tf]. The solver imposes the initial conditions at tspan(1), and integrates from tspan(1) to tspan(end). To obtain solutions at specific times (all increasing or all decreasing), use tspan = [t0,t1,...,tf].

For tspan vectors with two elements [t0 tf], the solver returns the solution evaluated at every integration step. For tspan vectors with more than two elements, the solver returns solutions evaluated at the given time points. The time values must be in order, either increasing or decreasing.

I also refer you to http://www.mathworks.com/help/techdoc/ref/odeset.html.

[Rose]

Thank you for your guidance. Let me try this and i'll let you know asap. Thanks once again.

[Rose]

I tried something like this:

finaltime = 6650;
deltat=730;
tspan = 0:deltat:finaltime;

y0=[10;13221;0;10];

options = odeset('RelTol',1e-12,'AbsTol',1e-12);

N=9;

for i=1:N

      tspan=[(i-1)*deltat, i*deltat];
      sol = ode15s(@para,tspan,y0,options);
      sol(2)=sol(2)-sol(2).*0.6;
      y0=sol;

end

X = sol.x

Y = (sol.y)'

plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3), X, Y(:,4))

[Rose]

but i got the error:

??? Index exceeds matrix dimensions.

Error in ==> mai at 30 sol(2)=sol(2)-sol(2).*0.6;

[Star Strider]

I do not know what the dimensions of ‘sol’ are, but you later refer to it as a structure ‘sol.x’ and ‘sol.y’, and then in your plot command you refer to Y being an [N x 4] matrix. I suspect you need to refer to it differently in the line that gives you the error. (I do not know what you want to do in that line so I cannot suggest a way to correct it.)