Solving a system of ODEs

Hello

I have a problem when trying solve this system of ODEs. I have tried both with the ode45 solver and by constructing a RK4 for loop as well. Both codes results in NaN answers.. I have checked the equations etc. a number of times but can't seem to find the issue. I would if you could help to see if I have missed something? Below is seen the system of two ODEs that I'm trying to solve along with the input and ICs. The system is to be solved simultaneously for R(t) and Pg(t).

Code 1 with ode solver(function file and script file):

%Function file
function diffeqs=ode_sys(t,var)
R=var(1); %dependent variable 1
Pg=var(2); %dependent variable 2


Pg0 = 5.6*10^6; %Saturation pressure(Pa)
Pa = 1.013*10^5; %Ambient pressure(Pa) 
Sigma = 0.045; %Surface tension(kg/s^2)
Rho = 80; %Density (kg/m^3)
Eta = 1.5*10^2; %Viscosity(Pa*s)
MCO2 = 0.044; %Molarmass (kg/mol)
D = 6.317e-10; %Diffusion coefficient(m^2/s)
KH = 1.839e-5; %Henry's konstant
T = 373.15; %Temperature(K)
Rmin = 2*Sigma/(Pg0-Pa); %Minimum radius
R0 = 1.8e-8; %R0 
a = 7.99e-3; %kg^2/s*m^4 


%dR/dt
diffeqs(1,1) = R*((Pg-Pa-2*Sigma/R)/4*Eta);
%dPg/dt
diffeqs(2,1) = a*((Pg0-Pg)^2/(Pg*R^3-Pg0*R0^3))-3*Pg*(diffeqs(1,1)/R); 
end

%Script file
dt=0.01;
trange=0:dt:300;
ICs=[1.8e-8,5.6*10^6]; %I.C

[tsol,varsol]=ode45(@ode_sys,trange,ICs);


figure
plot(tsol,varsol(:,1)*1e9); %plot af R
title('Radius vs t')
xlabel('t(s)')
xlim([0 300])
ylim([0 10000])
ylabel('radius(nm)')

figure
plot(tsol,varsol(:,2)); %plot af Pg
title('Pg vs t')
xlabel('t(s)')
xlim([0 300])
ylabel('Pg(Pa)')

Alternative method with RK4 approach:

Pg0 = 5.6*10^6; %Pressure saturation(Pa)
Pa = 1.013*10^5; %TAmbient pressure(Pa)
Sigma = 0.045; %Surface tension(kg/s^2)
Rho = 1120; %Density(kg/m^3)
Eta = 1.5*10^2; %Viskosity(Pa*s)
MCO2 = 0.044; %M of CO2 (kg/mol)
D = 6.317e-10; %Diffusion constant(m^2/s)
KH = 1.839e-5; %Henry's konstant
T = 373.15; %Temperature(K)

Rmin = 2*Sigma/(Pg0-Pa); %Minimum radius (m)
R0 = Rmin*1.1; %R0 (m)

a = 7.99e-3; %kg^2/s*m^4 (se mathcad)


%Function handle for R
fR=@(t,R,Pg) R*((Pg-Pa-2*Sigma/R)/4*Eta);
%Function handle for Pg
fPg=@(t,R,Pg) a*((Pg0-Pg)^2/(Pg*R^3-Pg0*R0^3))-3*Pg*(R*((Pg-Pa-2*Sigma/R)/(4*Eta))/R); 

%I.C
t(1) = 0;
R(1) = R0;
Pg(1) = Pg0;

dt=0.1;
tfinal = 300;
N=ceil(tfinal/dt);

%Update loop
for i=1:N
    %Update time
    t(i+1)=t(i)+dt;
    %Update R og Pg
    k1R =   fR(t(i),         R(i),           Pg(i)); 
    k1Pg =  fPg(t(i),        R(i),           Pg(i)); 
    k2R =   fR(t(i)+dt/2,    R(i)+dt/2*k1R,  Pg(i)+dt/2*k1Pg); 
    k2Pg =  fPg(t(i)+dt/2,   R(i)+dt/2*k1R,  Pg(i)+dt/2*k1Pg); 
    k3R =   fR(t(i)+dt/2,    R(i)+dt/2*k2R,  Pg(i)+dt/2*k2Pg); 
    k3Pg =  fPg(t(i)+dt/2,   R(i)+dt/2*k2R,  Pg(i)+dt/2*k2Pg); 
    k4R =   fR(t(i)+dt,      R(i)+dt*k3R,    Pg(i)+dt*k3Pg); 
    k4Pg =  fPg(t(i)+dt,     R(i)+dt*k3R,    Pg(i)+dt*k3Pg); 
    
    R(i+1) = R(i) + dt/6*(k1R + 2*k2R + 2*k3R + k4R);
    Pg(i+1) = Pg(i) + dt/6*(k1Pg + 2*k2Pg + 2*k3Pg + k4Pg);
end

Regards

Kaare

Answer 1

[Ameer Hamza]

There seems to be an issue with the way you implemented the ODE function. The initial conditions are

ICs=[1.8e-8,5.6*10^6]; %I.C: R = 1.8e-8, Pg = 5.6*10^6

and when you call ode45 with these initial conditions, then in first time-step the value in ode45 functions are defined like this

R=1.8e-8; %dependent variable 1
Pg=5.6*10^6; %dependent variable 2
..
..
..

Putting these values and other constants defined in your function into this term

((Pg0-Pg)^2/(Pg*R^3-Pg0*R0^3)) % 0/0 because Pg0 == Pg and R == R0

will result in a 0/0 condition, which creates NaN. You need to check if the ODE equations are correctly implemented in MATLAB.

[Kaare Hounisen]

Hello Ameer

Thanks for your answer. I also noticed that and changed the Pg0 in order to insure that it would not be 0/0. However, the result seems very wrong as the Pg(row 2) goes straight to the value of Pg0 and the R(row 1) elavates too quickly in its value?

Regards

Kaare

[Ameer Hamza]

There might also be the issue of numerical stability. The value of the states is quite extreme. The value of r is factional, whereas Pg is of the order of 10^5. Somehow, scaling the state values to bring them at a comparable range might improve the numerical stability of the ode solver.

[Kaare Hounisen]

I've looked it over again and there seems to be a problem with the units that I can't fix no matter what expression I use for Kh. Both expressions for dPg/dt should be Pa/s, however, one side seems to end up being Pa/m*s if I use kg/kg*Pa for Kh. I will have to revise the article again to see if any mistakes were made in the derivation. Thanks for you help.

[Ameer Hamza]

Yes, revisiting the derivation seems to be a good strategy. Also, make sure that you have an idea about the theoretical solution so you can verify whether the numerical solution is correct.