for loop for time dependent parameter values in ode solver only works for some values of t

Question

I'm using a simple if loop to change my parameter values within my ode script. Here is an example script I wrote that exhibits the same problem. So first the version which works:

function aah = al(t,x)

if (t>10000 && t<10300)
    ab = [0; 150];
else
    ab = [150; 0];
end

aah = [ab];

this can be run using

t = [0:1:10400];
x0 = [0,0];
[t,x] = ode23tb(@al, t,x0);

and visualised with

plot(t,x(:,1))
plot(t,x(:,2))

Ok that's the good version. Now if all you do is change t to

t = [0:1:12000];

the whole thing blows up. You might think it's just matlab averaging out the graph but it's not because if you look at

x(10300,2) 

the answer should be the same in both cases because the code hasn't changed. but this second version outputs 0, which is wrong!

What on earth is going on? Anyone got an idea?

Thank you so much for any help

Answer 1

Your function is constant (except 10000 < t < 10300), and therefore the internal solver starts to solve the system with very large time step, 10% of total time by default. (In the adaptive ODE solver, if the system is constant, higher order and lower order solution will give the same solution, and the (estimated) error will be zero. So the solver assumes that current time step is good enough.) You can see if you give tspan with just two element, start and end time.

t = [0 12000];

Usually the tspan does not affect to the internal time step of solver. The solvers solve the system with their internal time step, and then just interpolate at tspan given by the user. So if the internal time step unfortunately "leap over" the interval [10000, 10300], the solver won't know about the interval.

So you better set the maximum step size, relatively smaller than 300.

options = odeset('MaxStep', 10);
[t, x] = ode23tb(@al, t, x0, options);

If you don't want to solve with small step size whole time (and if you "know" when the function are not constant), you should solve separately.

t1 = [0 9990];
t2 = [9990 10310];
t3 = [10310 12000];

[T1, x1] = ode23tb(@al, t1, x0);
[T2, x2] = ode23tb(@al, t2, x1(end,:));
[T3, x3] = ode23tb(@al, t3, x2(end,:));

T = [T1; T2(2:end); T3(2:end)];
x = [x1; x2(2:end,:); x3(2:end,:)];