Simultaneous differential equations, memory issues.

Hi,

These are the equations I was trying to solve by Euler's method.

dM_q/dt=-gamma*M_q+N*B_q*(M_q+1)-c*k_q*M_q
dN/dt = P-A*N-N*sum(B_q*M_q)

Here q = may be from 100-500, and this is a block of the code I tried:

t1 = 0;                 %initial time
t2 = 1e-8;              %final time
n = 1e4;
dt = (t2-t1)/n;         %time step
%==============================================
N = 0;                              %laser inversion
for i = q1:mod_step:q2
  Mq((i-q1)/mod_step +1) = 1;     %initial photon number assignmnet in each mode
  Bq((i-q1)/mod_step+1) = B0/(1+((i-q0)/Q)^2);
  Kq((i-q1)/mod_step+1) = k0/((1+((i-q0)/w0)^2));
  mod((i-q1)/mod_step + 1,1) = i; %array of mode number
end
%==========================================
for t = 1:n
  TT = t1 + dt*(t-1);
  sumB = 0;

    for i = q1:mod_step:q2
        sumB = sumB + Bq((i-q1)/mod_step +1)*Mq((i-q1)/mod_step +1);
    end

    for i = q1:mod_step:q2
        Mq((i-q1)/mod_step+1) = Mq((i-q1)/mod_step+1) + dt*(-gamma*Mq((i-    q1)/mod_step+1) + Bq((i-q1)/mod_step+1)*N*(Mq((i-q1)/mod_step+1)+1) - Kq((i-q1)/mod_step+1)*c*Mq((i-q1)/mod_step+1)) ;

        if (rem(t,m1) == 0)   %extracting the output after m1 step
            M_q(t/m1,(i-q1)/mod_step+1) = Mq((i-q1)/mod_step+1);  
            T(t/m1,1) = TT;
        end

    end
    N = N + dt*(P - A*N - N*sumB);
end

[X1,Y1] = meshgrid(mod,T);
mesh(Y1,X1,M_q);
%============================================

The problem I am facing is: if I am decreasing 'dt' (to get a high accuracy), I got the error message, out if memory. I was told by one of my friend that he was able to solve this problem for 'dt = 10^-14' in Fortran. I want to run the program for a final time 't2 = a few millisec say, 10ms'. I could at the most go for 't2 = 1 microsec with dt = 10^-12', which eventually give me 'n=10^6', the iteration step of my for loop. But for 't2=10milisec', I have 'n=10^10 to 10^12', a huge number.

Does anybody know where am I doing the mistake, or how to speed up the problem ? I tried with ODE23 also, as suggested by one of the forum member, but ended up with the same problem.

Thanks for your time.

Edit: ODE23 solver:

Mq = ones(mode,1);
mod = [q1:mod_step:q2]';

N = 0;
X = [Mq;N];
%==========================================================================
TSPAN = [0 1.0E-7];
[T,X] = ode23(@ratequation,TSPAN,X);
[X1,Y1] = meshgrid(mod,T);
Mqq = X(1:end,1:end-1);
mesh(Y1,X1,Mqq);

[Andrew Newell]

There is no way that a do-it-yourself Euler method is going to outperform ode23. It would be better if you posted the code you used to call ode23.

[Graig]

The code to call ODE23 is posted in the question.

Here is my function file:

function laser = ratequation(T,X)

for i = q1:mod_step:q2

Bq((i-q1)/mod_step+1) = B0/(1+((i-q0)/Q)^2);

Kq((i-q1)/mod_step+1) = k0/((1+((i-q0)/w0)^2));

end

laser = [-gamma*X(1:end-1) + X(end)*Bq.*(X(1:end-1)+1) - c*Kq.*X(1:end-1); P - A*X(end) - X(end)*sum(Bq.*X(1:end-1))];

[Graig]

Any Idea?

Answer 1

[Walter Roberson]

It is not a good idea to overwrite the function "mod".

If you have any reasonably recent version of MATLAB, then you do not need to use meshgrid(): just pass the vectors in as the X and Y coordinates and the grid will be implied.

[Graig]

Right. Actually I can do that also. But I don't think it is the part making the program slow !!