Coupled Differential Equation

Dear Friend,

I am solving a coupled differential equation in matlab to simulate the laser rate equation. The number of differential equation depends on the number of modes I will put through input, usually it is a very high number say 500-600.

I have two for loops inside another for loop.

I use the usual trick: dx/dt = Ax (say this is the differential equation)

so x2 = x1 + dt(Ax) (I solve it this way giving an initial condition on x1)

The problem is, this equation will be valid as long as abs(dx/x)<<1. And we need 'dt' for this purpose very small, which will eventually increases the iteration of my for loop. Now when I do that, I got an out of memory error.

Is there a way to get rid of it. I was thinking of extracting few outputs from the iteration (not all the output), but also it didn't work.

Thanks for your time.

-Graig

[Jarrod Rivituso]

Forgive me if I misunderstand, but are you not using one of the ODE solvers?

[Graig]

No, I was not using ODE solution. I am not sure, whether I can use ODE for such a high number of differential equation. Will glad to know about it.

Answer 1

[Andrew Newell]

Here is an example where someone is using ode45 to solve an even bigger problem. To avoid out-of-memory errors, make A into a sparse matrix if you can.

EDIT: If I understand your notation, M is actually a vector with components M_q, and ditto for K and B. So you should create a new vector

X = [M; N];

and then define your function as follows:

f = @(x) [-gamma*x(1:end-1) + x(end)*B.*(x(1:end-1)+1) - c*K.*x(1:end-1); P - A*x(end) - x(end)*sum(B.*x(1:end-1))];

(or you could create an easier-to-read function in a separate file).

Answer 2

[Jarrod Rivituso]

I would recommend giving these a try:

http://www.mathworks.com/help/techdoc/ref/ode23.html

I've never used them with the number of states you have. You may end up finding that your system is "stiff", depending on the system dynamics. In this case, you may prefer ode15s or other stiff solvers to the commonly used ode45.

Answer 3

[Graig]

Thanks guys. But I don't know how to put this solver in a loop. Need to think about it.

[Jarrod Rivituso]

You shouldn't need to put this in a loop just to solve a system of equations. The idea is that you tell it

- the initial conditions

- a function that returns the state derivatives

- a time span

and then it will go solve the equations for you.

Is there something else that is forcing you to use a loop?

[Graig]

The equations are like these:

dM_q/dt = -gamma*M_q + B_q*N*(M_q+1) - K_q*c*M_q

dN/dt = P - A*N - N*sum(B_q*M_q)

where B_q = B0/1+((q-q0)/Q)^2

and K_q = K0/1+((q-q0)/w)^2

Few details:

q: no of modes, so M_q will be the photon numbers in mode q. In this case q is more than 500.

Term giving me trouble in ODE solver is sum(B_q*M_q). I think this should be inside a loop, as we are adding more and more modes, this term is keeping considering these modes.

[Andrew Newell]

I have edited my answer to take this information into account.

[Graig]

Rather than Euler's method, I can at least go to a higher order in time using ODE23. But above that I still have the 'out of memory' error. I tried a few other solver also without success.

[Andrew Newell]

What do you mean by "above that"?

[Andrew Newell]

This question is starting to sprawl and it is not clear what is going on. I recommend you accept an answer for this question and submit a new one in which you show us the code you used and describe the problem you are having.

[Graig]

Above that I meant 't_final'. I could solve the equation when I am limited to TSPAN = [0 1.0E-6] when the other parameters e.g., gamma and so on are set to some value. If I increase the final time, meaning 1.0E-6 to say 1.0E-4 it went to the error message. I am bound to keep the input parameters (gamma, P etc) fix. Physically this should work for all time interval.

Answer 4

[Graig]

Here is a discussion I found recently on this topic.