How to solve 'out of memory' problem while solving coupled differential equation ?

Hello Everyone

I am trying to solve a set of coupled differential equations . It has four variables that are coupled together. The code is correct to my knowledge. The code is as follows--

clc
clear 
close all

T_ar=300;                                         
k_B=1.380649e-23;
e=1.60217663e-19;                                 
m_e=9.109e-31;                                     
epsilon=8.854e-12;                                 
A_sp=.005024;                                      
A_sg=.007854;                                      
l_B=0.075;                                         
n=1.96e16;                                        

upsilon=0.08;                                      
T_e=3;                                             

m_i=40;                                            
amu=1.660539e-27;                               
v_B=sqrt((e*T_e)/(m_i*amu));                      
v_e = sqrt(8*k_B*T_ar/(pi*m_e));                  

p=1;                                               
n_g=(p/(k_B*T_ar));                                                                               
Km=10^-13;
v_m=Km*n_g;                                       
v_eff=v_m+(v_e/l_B);                               


J_e=e*n*v_e;                                       
J_i=e*n*v_B;                                      

Jip=J_i*A_sp;                                     
Jig=J_i*A_sg;                                      
Jep=J_e*A_sp;
Jeg=J_e*A_sg;

                                    
L_pl=(l_B*m_e)/(e^2*n*A_sp);                       
R_pl=v_eff*L_pl ;                                 
Lpl_inv=1/L_pl;

A=(1/(2*e*epsilon*n*A_sp^2));
xx=(e*A/(T_e));
B=(1/(2*e*epsilon*n*A_sg^2));

V0=300;                                           
C_B=3.34e-6;                                      


f = @(t,x) [-x(4)-Jip+Jep*exp(-(e*A/(T_e))*(x(1)^2));
            x(4)-Jig+Jeg*exp(-(e*B/(T_e)).*(x(2)^2));
            -(1/C_B).*x(4);
            (A.*(x(1)^2)-B.*(x(2)^2)+(V0*cos(t))+x(3))-((v_eff/Lpl_inv).*x(4))];
[t,xa] = ode45(f,[0, 0.1],[0,0,-62,0]);


figure(1)
plot(t,xa(:,1),'Color','blue'),grid on;
title('xa1')
xlabel('t'), ylabel('Qsp (Powered Electrode Sheath Charge)')

figure(2)
hold off
plot(t,xa(:,2),'Color','blue'),grid on;
title('xa2')
xlabel('t'), ylabel('Qsp (Grounded Electrode Sheath Charge)')

figure(3)
plot(t,xa(:,3),'Color','blue')
grid on;
title('xa3')
xlabel('t'), ylabel('Self Bias Voltage')

figure(4)
plot(t,xa(:,4),'Color','blue'),grid on;
title('xa4')
xlabel('t'), ylabel('Plasma Current')

Now while the timespan is set to 0 to 0.1 the code runs. But out of 4 variables ony two of them reaches steady state. If I set the timespan longer like 0 to 1 or 0 to 5 , after a while of calculation matlab shows out of memory error. one such error notification is--

out of memory

Error in ode45 (line 475) yout = [yout, zeros(neq, chunk, 'like', yout)]; %#ok<AGROW>

I have a 16 GB ram pc. I am more or less certain adding 16 GB more ram wont remedy the problem. I am confident that after 5 seconds all of the variabls would reach equilibrium. Is there any way that the calculation data upto 5 seconds wont be saved on my ram and after 5 seconds it would be saved in my ram. Or maybe any other effective method which can solve my problem.

Thank You.

Answer 1

[Sam Chak]

Hi @Golam

I have a suspicion that the main issue lies in the state equation of , possibly due to the excessively large values of A and B. Unfortunately, my simulation settings on the MATLAB forum platform are limited to a maximum of 0.06 seconds. However, if you utilize the Approximated ODEs, which exclude the dynamics of and , you should be able to run the simulation smoothly. The reason for this omission is that depends solely on , and converges to zero. Additionally, the time responses of and appear to be increasing linearly.

%% Parameters
T_ar    = 300;                                         
k_B     = 1.380649e-23;
e       = 1.60217663e-19;                                 
m_e     = 9.109e-31;                                     
epsilon = 8.854e-12;                                 
A_sp    = .005024;                                      
A_sg    = .007854;                                      
l_B     = 0.075;                                         
n       = 1.96e16;                                        
upsilon = 0.08;                                      
T_e     = 3;                                             
m_i     = 40;                                            
amu     = 1.660539e-27;                               
v_B     = sqrt((e*T_e)/(m_i*amu));                
v_e     = sqrt(8*k_B*T_ar/(pi*m_e));              
p       = 1;                                               
n_g     = (p/(k_B*T_ar));                                                                               
Km      = 10^-13;
v_m     = Km*n_g;                                       
v_eff   = v_m + (v_e/l_B);                               
J_e     = e*n*v_e;                                       
J_i     = e*n*v_B;                                      
Jip     = J_i*A_sp;
Jig     = J_i*A_sg;
Jep     = J_e*A_sp;
Jeg     = J_e*A_sg;                 
L_pl    = (l_B*m_e)/(e^2*n*A_sp);                       
R_pl    = v_eff*L_pl;
Lpl_inv = 1/L_pl;
A       = (1/(2*e*epsilon*n*A_sp^2))

A = 7.1247e+17

xx      = (e*A/(T_e));
B       = (1/(2*e*epsilon*n*A_sg^2))

B = 2.9153e+17

V0      = 300;
C_B     = 3.34e-6;

%% Original ODEs (x3 dynamics depends on x4 only, and x4 converges to 0)
f       = @(t, x) [- x(4) - Jip + Jep*exp(- (e*A/T_e)*x(1)^2);
                     x(4) - Jig + Jeg*exp(- (e*B/T_e)*x(2)^2);
                   - (1/C_B)*x(4);
                   A*x(1)^2 - B*x(2)^2 + V0*cos(t) + x(3) - v_eff/Lpl_inv*x(4)];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, xa] = ode45(f, [0:0.0001:0.06], [0, 0, -62, 0], options);

figure(1)
subplot(2,2,1), plot(t, xa(:,1)), grid on, xlabel('t'), title('x_{1}')
subplot(2,2,2), plot(t, xa(:,2)), grid on, xlabel('t'), title('x_{2}')
subplot(2,2,3), plot(t, xa(:,3)), grid on, xlabel('t'), title('x_{3}')
subplot(2,2,4), plot(t, xa(:,4)), grid on, xlabel('t'), title('x_{4}')

%% Approximated ODEs (the dynamics of x3 and x4 are omitted)
f       = @(t, x) [- Jip + Jep*exp(- (e*A/T_e)*x(1)^2);
                   - Jig + Jeg*exp(- (e*B/T_e)*x(2)^2)];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, xa] = ode45(f, [0:0.0001:0.1], [0, 0], options);

figure(2)
subplot(2,1,1), plot(t, xa(:,1)), grid on, xlabel('t'), title('x_{1}')
subplot(2,1,2), plot(t, xa(:,2)), grid on, xlabel('t'), title('x_{2}')

%% Check initial value of the state equations
% t       = 0;
% x       = [0, 0, -62, 0];
% f(t, x)

Answer 2

[Walter Roberson]

You have what is probably a stiff system, but you are using ode45() for it instead of using a stiff solver such as ode15s()

You are using a tspan that is two elements long, so you are implicitly asking for it to output whenever it feels like it.

In practice "whenever it feels like it" will be (default settings) four points every time that a step is "accepted". Because it is a stiff system, lots and lots and lots of points will be accepted...

You should probably switch to ode15s() instead of ode45()

But failing that, you should switch to a tspan vector that is more than two elements; the output will only be at those particular locations.

[Golam]

@Walter Roberson Thanks for your valuable suggestion. I will let you know if this works.

Answer 3

[Torsten]

ode45 collects too many outputs if it saves the solution for every successful step (as it is done if you set tspan = [0 0.1], e.g.)

Use something like

tspan = linspace(0,1,500);
[t,xa] = ode45(f,tspan,[0,0,-62,0]);

instead.

And since your system deals with reaction rates, I guess your equations are stiff. So use ode15s instead of ode45.

[Golam]

@Torsten Thanks for your valuable suggestion. I will let you know if this works.