Simultaneous differential equations, memory issues.
1 次查看(过去 30 天)
显示 更早的评论
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);
回答(1 个)
Walter Roberson
2011-5-12
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!