Can you use ODE45 in a for loop?

Hi there, I am having to solve a pair of coupled 1st ODEs. I also need to change the integration time like this:

ti = 0
tf = 10
timestep = 0.1
num    = (tf-ti)/timestep;
timerange = linspace(ts,tf,num);

for i = 1:length(timespan)

        tt = timerange(i); %defining the integration time start.
        [t,x]=ode45('func',[tt:0.01:tf],[0,0]);  %solving the coupled odes.

end

So what is happening is that I am changing the start integration in the [T0 TFINAL] vector. Is there a way that I can save all of the t's and the x's as a matrix/for each iteration so that I can plot all of the solutions to the coupled odes?

Thanks

Answer 1

[Matt Tearle]

Given that ode45 is a variable step solver, you don't know how many t and x values you'll get each time, so the simplest solution would be to save them all in cell arrays:

n = length(timerange);
allx = cell(n,1);
allt = cell(n,1);
for i = 1:n
   ...
   allx{i} = x;
   allt{i} = t;
end

But if all you're trying to do is plot them all, why not just plot each one as you go? Use a hold on or hold all and then plot(t,x) in the loop.

(Also, [t,x] = ode45(@func,[tt... -- function handles are cooler than strings :) )

[Martin]

Thanks, I was thinking cell arrays might be useful. Basically once I have all of the t's and x' I need to do a few things on them before plotting, I'll give it a go and let you know

Thanks

[Matt Tearle]

Each integration is independent, so I suspect you can do all your analysis and post-processing each time in the loop.... unless, of course, you actually need to compare between runs as part of your post-processing. Saving all the results, then processing them in a loop is slightly less efficient than doing everything in a single loop and throwing away the results from each run once you're done with them. But... a cell array of 100 elements is no great strain on your system, so it's not worth stressing over.

[Martin]

Good point, basically I only want to include the solutions to the differential equations up to a point where the sign of the solution changes for the first time(if that makes sense) so I was thinking of adding a continue clause in it to stop unnecessary integration because the solution switches sign a lot over the total timespan.

[Martin]

actually that can't be done unless I change the ode45.m file! Whoops!

[Matt Tearle]

Actually, it can! Define an event function and use odeset to tell ode45 to look for the event(s) defined by that function.

function [val,isterm,dir] = signchange(t,x)
val = x;    % look for x to go through 0
isterm = 1; % stop when it happens
dir = 0;    % no matter which direction (+ -> -  or - -> +)

opts = odeset('Events',@signchange);
[t,x] = ode45(...,opts);

(Here I'm assuming x is a scalar, and you don't care which direction the sign change occurs in.)

[Martin]

This is really helpful. I'm going to look this up and see if I can make the for loop more efficient. Thanks, I'll let you know how it goes =]

[Martin]

so is the event function defined separately as a .m file. and then called in to the script where the ode45 solver is (with the fact there is already a .m function for the ODEs..)

[Matt Tearle]

Right. So you'd have func.m that defines the ODEs and signchange.m that defines the event (as above). Then use odeset to "attach" signchange.m to the ODEs as an event (the second to last line of code in my previous comment); this creates a structure of ODE options (that I called opts). Pass this to ode45 as the last input.

[Martin]

Ahh I see, thanks Matt. That is really really helpful and I appreciate your time taken to explain this!

[Martin]

I'm getting a problem:

"SWITCH expression must be a scalar or string constant."

I read what the switch function does, and I'm assuming in the lines you wrote for the signchange.m function requires x to be a number, which I'm pretty sure it is if I add the line:

function [val,isterm,dir] = signchange(t,x)
xx = x(1); %making element 1 of the vector x
val = xx;    % look for x to go through 0
isterm = 1; % stop when it happens
dir = 0;    % no matter which direction (+ -> -  or - -> +)

I did this because x will actually be a two element vector due to being coupled odes with position in the x(1,1) not x(1,2) element? This still doesn't seem to work though...

[Matt Tearle]

Doesn't work how? Do you get the same error message (about switch) or something else? You might have to give the whole error message traceback. What you've done is correct -- if x is a vector, then you need to pull out a specific scalar value that you're looking to go through 0. You could also just do val = x(1); of course.

[Martin]

Basically I press run on the program and it immediately brings the error message up, so I am unsure what it is doing. This is the full error message I get.

 SWITCH expression must be a scalar or string constant.

Error in odeevents (line 32)

    switch lower(eventFcn)

Error in ode45 (line 148)
[haveEventFcn,eventFcn,eventArgs,valt,teout,yeout,ieout] = ...

Error in ode45testfunc2 (line 39)
          [t,x]=ode45('cwfield',tt:0.01:tf,[0,0],opts);  %solving the coupled odes.

Maybe the problem is is that some solutions never pass through zero? other than that I was thinking of having to do it the old way I was but using the cellfun function.

[Matt Tearle]

I'm guessing you're using the code I provided to define the event (using a function handle in odeset). If so, the problem might be that you're using a string to define the ODE rate equations and a function handle to define the event function. Try changing your call to ode45 to

[t,x]=ode45(@cwfield,tt:0.01:tf,[0,0],opts);

(BTW, not going through 0 isn't a problem -- it just means the event won't ever be triggered, so ode45 will integrate until tf).

[Martin]

Yea I thought that it was because of using strings rather than handles...but when I use the @cwfield I get a different error! I have managed to get around this in the end using a bit of brute force, but I'll return to this when I have some time on Monday.

Thanks for all of the help, I am getting some really nice results from the code now.