MATLAB: Solving a system of ODEs

MATLABnanodeode45system of ode

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

Best Answer

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.
Related Question