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

5 views (last 30 days)
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.

Accepted Answer

Sam Chak
Sam Chak on 5 Feb 2024
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)

More Answers (2)

Walter Roberson
Walter Roberson on 4 Feb 2024
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.

Torsten
Torsten on 4 Feb 2024
Edited: Torsten on 4 Feb 2024
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.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!