Is it Possible to Solve a Coupled System of ODEs with ODE45 in MatLAB?

[Niles]

Hi

I am trying to solve a simple set of coupled ODE's by ODE45. The coupled system is given by:

function xprime = eoms(t, x)

xprime = [
    1e9 + 5.0e4*x(3) - 50*x(1);
    4.0e1*x(1) - 3.3e3*x(2);
    2.0e3*x(2) - 5e4*x(3) + 3.5e7*heaviside(t-1)*x(4);
    1.0e3*x(2) - heaviside(t-1)*5.0e7*x(4)];

I solve it using the following command:

x0 = [0 0 0 0];
tspan = [0, 2];
[t, x] = ode45(@eoms, tspan, x0);

However when I compile MatLAB just keeps calculating, it doesn't give me a result. Maybe it is due to the very rapid rates in the equations. Do I have any options here, or am I not able to solve for the transient behavior?

Thanks in advance.


Niles.

 

[AlephZero]

Try mkaing the time interval much smaller, for example 0 to 1e-8, and see what happens. If that works, try increasing the interval (say by factors of 10) till it blows up.

Since you have a constant of 1e9 in your definition of xprime, it's likely the solution contains functions similar to exp(1e9 t) or exp(-1e9 t). A conditionally stable integration method like ODE45 may need need to take time steps of smaller than 1e-9 to avoid going unstable. It may not be able to find a suitable step length, and even if it does, more than 10⁹ steps may take a long time.

You mgiht have more luck changing ODE45 to one of the solver names ending in "s", which should work with a bigger step size without blowing up.

 

[Niles]

Try mkaing the time interval much smaller, for example 0 to 1e-8, and see what happens. If that works, try increasing the interval (say by factors of 10) till it blows up.

Since you have a constant of 1e9 in your definition of xprime, it's likely the solution contains functions similar to exp(1e9 t) or exp(-1e9 t). A conditionally stable integration method like ODE45 may need need to take time steps of smaller than 1e-9 to avoid going unstable. It may not be able to find a suitable step length, and even if it does, more than 10⁹ steps may take a long time.

You mgiht have more luck changing ODE45 to one of the solver names ending in "s", which should work with a bigger step size without blowing up.

Thanks for taking the time to reply. I have actually already tried your first suggestion, and it blows up as soon as the Heaviside step-function is different from 0. I also suspected that. I tried changing ODE45 to odes23, and now the solution pops up almost instantly! Wow, that is really good. Thanks!Niles.

 

[AlephZero]

I have actually already tried your first suggestion, and it blows up as soon as the Heaviside step-function is different from 0.

In that case, you can probably solve it with ODE45 by splitting it into two parts 0 to 1 and 1 to 2, so you force one time point to be exactly on the "edge" of the step function at t = 1.

You might need to make two versions of xprime with and without the step function included, so in effect you have heavisde(0) = 0 at the end of the first half of the solution, and heavisde(0) = 1 at the start of the second half.

That might give you a more accurate solution than odes23, which will "round off" the edge of the step a bit, in order to keep going.

 

[Niles]

Thanks for helping, that is kind of you. I tried extending the system of ODEs, but I get the message:

Warning: Failure at t=1.000000e+000. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-015) at time t.

Something tells me choosing a different solver won't help me here. And I can't even change the time-step. Do I have any alternatives left?


Niles.