Asked by Niles Martinsen
on 11 Sep 2012

Hi

I have a relatively simple set of coupled ODE's that I am trying to solve by ODE45. The following is the content of the .m-file, which shows the coupled system:

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. I think it has something to due with the fact that the Heaviside step-function, because if I take it out then it nicely solves the remaining part. I have also tried in a competing software, and the same thing happens (actually, there the solver just gives an error).

Maybe it is also due to the very rapid rates in the equations. Anyhow, I am really lost on what to do - unfortunately I need the transient behavior, so I can't just assume steady-state for some of the equations.

Do I have any options here?

Thanks in advance.

Best, Niles.

*No products are associated with this question.*

Answer by Star Strider
on 11 Sep 2012

Edited by Star Strider
on 11 Sep 2012

Accepted answer

You have a ‘stiff’ system, and `ode45` is not the best option for it, although it's an appropriate initial experiment. If `ode45` seems not to work (and the system otherwise seems to be solvable in some situations as yours was), switching to `ode15s` is appropriate, especially in systems with coefficients of widely-ranging magniudes such as yours are (from 4E+1 to 1E+9).

After experiencing the same problem you did with `ode45`, I ran your code with `ode15s` and it solved your `eoms` system in 1.2225 seconds.

So I suggest you simply change your ODE function call to:

[t, x] = ode15s(@eoms, tspan, x0);

and enjoy the results!

Also, I coded your function as an implicit function to run it, avoiding the need to save a separate function `m-file`. (This isn't always possible, but it works in many situations.)

The function:

eoms = @(t,x) [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)];

and the call to it changes to:

[t, x] = ode15s(eoms, tspan, x0);

Answer by Shashank
on 28 Jun 2013

Star Strider, could you please help me solve this set of four coupled ODEs?

tpcl = @(xi,y) [y(2); ((1 + (y(2)*y(2)))^(1.5))*(y(3) - (A/(y(1)^3)))/sigma; -3*mu_l*y(4)/(rho_l*h_lv*(y(1)^3)); (Twall - Tsat - (Tsat*y(3)/(rho_l*h_lv)))/((y(1)/k_l) + R_int)];

xi_step = 0:1e-8:xi_end;

[XI,Y] = ode15s(tpcl,xi_step,[delta_ad deltaprime_ad delta_p_ad Q_ad]);

Here, apart from the variables y(1), y(2), y(3) and y(4), everything is a constant. I tried ode45 and ode15 but it still did not work. The results look unrealistic to say the least.

Opportunities for recent engineering grads.

## 0 Comments