Having trouble solving two second order ODEs simultaneously

I am trying to solve Keller Miksis equation for a double bubble system but I am having trouble breaking down the equation into two first order ODEs or solving these two equations simulatenously.

For ease, I am omitting the terms containing (dPsi/dt) and (dPsj/dt) from both equations.

The code I used is:

syms Ri(t) Rj(t) x(t) DRi(t) DRj(t) Dx(t)
c=350
Psi = 1000
Psj = 1000
Dij = 1e-6
rho = 997

ode1= (1-(1/c)*diff(Ri,t))*Ri*(diff(Ri,t,2))+ 1.5*(1-(1/(3*c))*diff(Ri,t))*((diff(Ri,t))^2) == (1/rho)*(1+(1/c)*diff(Ri,t))*Psi - (1/Dij)*(2*diff(Rj,t)*diff(Rj,t,2)^2+ diff(Rj,t,2)*diff(Rj,t)^2 )
[ode1, var1] = reduceDifferentialOrder(ode1,Ri);

ode2 = (1-(1/c)*diff(Rj,t))*Rj*(diff(Rj,t,2))+ 1.5*(1-(1/(3*c))*diff(Rj,t))*((diff(Rj,t))^2) == (1/rho)*(1+(1/c)*diff(Rj,t))*Psj - (1/Dij)*(2*diff(Ri,t)*diff(Ri,t,2)^2+ diff(Ri,t,2)*diff(Ri,t)^2 )
[ode2, var2] = reduceDifferentialOrder(ode2,Rj);

ode = [ode1; ode2];
var = [var1; var2];

[odes, vars] = odeToVectorField(ode);

ode_fun = matlabFunction(odes, 'Vars',{'t','Y'});

y0 = [0 0];
tspan = [0 5];
[t,y] = ode45(ode_fun,tspan,y0);

plot(t,R)

The error I am getting is:

I would really appreciate if you could help me figure out to solve these simultaenous second order ODEs.

Answer 1

[Torsten]

According to your equations in the graphics, these terms in the MATLAB code are wrong:

- (1/Dij)*(2*diff(Rj,t)*diff(Rj,t,2)^2+ diff(Rj,t,2)*diff(Rj,t)^2 )
- (1/Dij)*(2*diff(Ri,t)*diff(Ri,t,2)^2+ diff(Ri,t,2)*diff(Ri,t)^2 )

[Shishir Raut]

Sorry. I corrected it but the error is still there. What should I do?

[Torsten]

syms Ri(t) Rj(t) 
c=350;
Psi = 1000;
Psj = 1000;
Dij = 1e-6;
rho = 997;

ode1= (1-(1/c)*diff(Ri,t))*Ri*(diff(Ri,t,2))+ 1.5*(1-(1/(3*c))*diff(Ri,t))*((diff(Ri,t))^2) == (1/rho)*(1+(1/c)*diff(Ri,t))*Psi - (1/Dij)*(2*Rj*diff(Rj,t)^2 + Rj^2*diff(Rj,t,2) );
ode2 = (1-(1/c)*diff(Rj,t))*Rj*(diff(Rj,t,2))+ 1.5*(1-(1/(3*c))*diff(Rj,t))*((diff(Rj,t))^2) == (1/rho)*(1+(1/c)*diff(Rj,t))*Psj - (1/Dij)*(2*Ri*diff(Ri,t)^2 + Ri^2*diff(Ri,t,2) );
ode = [ode1; ode2];
[V,S] = odeToVectorField(ode)

V = 

S = 

[Shishir Raut]

Thank you so much, I was able to get the result when I ran it in my professor's laptop, I have on older version of MATLAB or maybe there is some issue I am unaware. Unrelated to my original query, I wanted to ask if it was possible to run ode45 function in the newly obtained system of first order ODEs and how I could obtain it? I would really appreciate your help.

[Torsten]

Yes, supply initial values for Rj, dRj/dt, Ri and dRi/dt in the order fixed in S and use ode45 or another numerical integrator to solve:

syms Ri(t) Rj(t) 
c=350;
Psi = 1000;
Psj = 1000;
Dij = 1e-6;
rho = 997;
ode1= (1-(1/c)*diff(Ri,t))*Ri*(diff(Ri,t,2))+ 1.5*(1-(1/(3*c))*diff(Ri,t))*((diff(Ri,t))^2) == (1/rho)*(1+(1/c)*diff(Ri,t))*Psi - (1/Dij)*(2*Rj*diff(Rj,t)^2 + Rj^2*diff(Rj,t,2) );
ode2 = (1-(1/c)*diff(Rj,t))*Rj*(diff(Rj,t,2))+ 1.5*(1-(1/(3*c))*diff(Rj,t))*((diff(Rj,t))^2) == (1/rho)*(1+(1/c)*diff(Rj,t))*Psj - (1/Dij)*(2*Ri*diff(Ri,t)^2 + Ri^2*diff(Ri,t,2) );
ode = [ode1; ode2];
[V,S] = odeToVectorField(ode);
M = matlabFunction(V,'vars',{'t','Y'});
y0 = [1 1 0.5 1];
tspan = [0 5];
[t,y] = ode45(M,tspan,y0);
plot(t,[y(:,1),y(:,3)])
grid on

[Shishir Raut]

Thank you so much! You are a lifesaver!

[Shishir Raut]

Hello, sorry to bother you again. For the equations in the question, I tried to insert the definitions of Psi and Psj into the equations.

This time I am not omitting the dPsi/dt and dPsj/dt terms and thus wrote the code as follows:

syms Ri(t) Rj(t) 
c = 1480;                                              
pv = 2338;                                                         
p2 = 101325;                                                                                                                   
rho = 1000 ;                                                         
sigmaL = 0.072 ;                                                       
muL = 1.005e-3 ;                                                                                                               
gamma = 1.4 ;    
Dij = 1e-6;
P0 = 101325
Roi = 100*1e-6
Roj = 50*1e-6
tau = 5

ode1= (1-(1/c)*diff(Ri,t))*Ri*(diff(Ri,t,2))+ 1.5*(1-(1/(3*c))*diff(Ri,t))*((diff(Ri,t))^2) == (1/rho)*(1+(1/c)*diff(Ri,t))*((P0+2*sigmaL/Roi)*((Roi/Ri)^(3*gamma))-2*sigmaL/Ri - 4*muL*diff(Ri,t)/Ri - (pv + ((p2-pv)*t)/tau)) +(Ri/(rho*c))*diff(((P0+2*sigmaL/Roi)*((Roi/Ri)^(3*gamma))-2*sigmaL/Ri - 4*muL*diff(Ri,t)/Ri - (pv + ((p2-pv)*t)/tau)),t) - (1/Dij)*(2*Rj*diff(Rj,t)^2 + Rj^2*diff(Rj,t,2) );
ode2 = (1-(1/c)*diff(Rj,t))*Rj*(diff(Rj,t,2))+ 1.5*(1-(1/(3*c))*diff(Rj,t))*((diff(Rj,t))^2) == (1/rho)*(1+(1/c)*diff(Rj,t))*((P0+2*sigmaL/Roj)*((Roj/Rj)^(3*gamma))-2*sigmaL/Rj - 4*muL*diff(Rj,t)/Rj - (pv + ((p2-pv)*t)/tau)) +(Rj/(rho*c))*diff(((P0+2*sigmaL/Roj)*((Roj/Rj)^(3*gamma))-2*sigmaL/Rj - 4*muL*diff(Rj,t)/Rj - (pv + ((p2-pv)*t)/tau)),t) - (1/Dij)*(2*Ri*diff(Ri,t)^2 + Ri^2*diff(Ri,t,2) );
ode = [ode1; ode2];
[V,S] = odeToVectorField(ode)

M = matlabFunction(V,'vars',{'t','Y'})
y0 = [Roj 0 Roi 0];
tspan = [0 5];
[t,y] = ode45(M,tspan,y0);
plot(t,[y(:,1),y(:,3)])
grid on

% psi = (P0+2*sigmaL/Roi)*((Roi/Ri)^(3*gamma))-2*sigmaL/Ri - 4*muL*diff(Ri,t)/Ri - (pv + ((p2-pv)*t)/tau)
% psj = (P0+2*sigmaL/Roj)*((Roj/Rj)^(3*gamma))-2*sigmaL/Rj - 4*muL*diff(Rj,t)/Rj - (pv + ((p2-pv)*t)/tau)

I am supposed to set dRi/dt and dRj/dt as 0 as the initial conditions. When I do that, the code runs infinitely and I can't seem to resolve that.I would really appreciate your help.

PS: This is an attempt to replicate the work done by a specific paper and I have used their equations as such. I have been able to obtain similar nature of graph for single bubble case.

This is the paper for your reference. I would really appreciate if I could get some help.

[Torsten]

Use

[t,y] = ode15s(M,tspan,y0);

instead of

[t,y] = ode45(M,tspan,y0);

But the result doesn't look promising.

[Shishir Raut]

Yeah, the graphs get flatlined completely.

I used the following code for single bubble case and the graph is very similar to that of the paper.

clear,clc
T0      = 298 ;                                                         % Ambient temperature, K
freq    = 26500;                                                        % Driving frequency, Hz
r0      = 5 ;                                                           % Initial bubble radius, microns
p0      = 1 ;                                                           % Ambient pressure, atm
pA      = 1.2 ;                                                         % Supplied oscillating pressure, atm
rho0    = 997 ;                                                         % liquid density, kg m-3
sigmaL  = 0.072 ;                                                       % surface tension, N m-1
muL     = 1e-3 ;                                                        % liquid viscosity, kg m-1 s-1
c       = 1485 ;                                                        % speed of sound, m s-1
gamma   = 1.8 ;                                                         % air specific heats ratio, -
tspan   = [0 1/freq] ;                                               
IC      = [r0*1E-6,0,p0*101325+2*sigmaL/(r0*1E-6),T0] ;                 
odefun  = @(t,w)keller_miksis(t,w,r0*1E-6,T0,rho0,freq,p0*101325,pA*101325,sigmaL,muL,c,gamma) ;
[time,SOL]  = ode45(odefun,tspan,IC) ;
t = time ;                                                             
R = SOL(:,1)*1E+6 ;                                                         
S = SOL(:,2) ;                                                          %Local velocity                                                    
p = SOL(:,3)/101325 ;                                                      
T = SOL(:,4) ;    


plot(t,R,'LineWidth',1.5)
title('Bubble radius')
xlabel('t','FontSize', 12, 'FontWeight', 'bold')
ylabel('R, \mu m','FontSize', 12, 'FontWeight', 'bold')

function f = keller_miksis(t,w,r0,T0,rho0,freq,p0,pA,sigmaL,muL,c,gamma)
R   = w(1) ; 
dR  = w(2) ;
p   = w(3) ;
T   = w(4) ;
dpdt = -3*gamma*p*dR/R ;
dTdt = 3*(1-gamma)*T*dR/R ;
pinf = (p0-pA)*t ;
pL = (p0+2*sigmaL/r0)*(r0/R)^(3*gamma)-2*sigmaL/R-4*muL*dR/R ;
ddR = ((1+dR/c)*(pL-pinf)/(rho0)-1.5*(1-dR/(3*c))*(dR^2)+(R/(rho0*c))*(dpdt+(2*sigmaL*dR+4*muL*dR^2)/(R^2)))/((1-dR/c)*R+4*muL/(rho0*c)) ;
f = [dR ; ddR ; dpdt ; dTdt];                                                      
end

I wonder how I can get the double bubble case to work. Thanks for the help nonetheless. I am able to obtain a set of first order ODEs at the very least, I will try to see how I can proceed from there.

[Torsten]

Maybe better to read.

syms Ri(t) Rj(t) psi(t) psj(t)
c = 1480;                                              
pv = 2338;                                                         
p2 = 101325;                                                                                                                   
rho = 1000 ;                                                         
sigmaL = 0.072 ;                                                       
muL = 1.005e-3 ;                                                                                                               
gamma = 1.4 ;    
Dij = 1e-6;
P0 = 101325;
Roi = 100*1e-6;
Roj = 50*1e-6;
tau = 5;
psi = (P0+2*sigmaL/Roi)*(Roi/Ri)^(3*gamma) -2*sigmaL/Ri - 4*muL*diff(Ri,t)/Ri - (pv + ((p2-pv)*t)/tau);
psj = (P0+2*sigmaL/Roj)*(Roj/Rj)^(3*gamma) -2*sigmaL/Rj - 4*muL*diff(Rj,t)/Rj - (pv + ((p2-pv)*t)/tau);
ode1= (1-(1/c)*diff(Ri,t))*Ri*(diff(Ri,t,2))+ 1.5*(1-(1/(3*c))*diff(Ri,t))*((diff(Ri,t))^2) == 1/rho * (psi + 1/c*diff(Ri*psi,t)) - (1/Dij)*(2*Rj*diff(Rj,t)^2 + Rj^2*diff(Rj,t,2) );
ode2 = (1-(1/c)*diff(Rj,t))*Rj*(diff(Rj,t,2))+ 1.5*(1-(1/(3*c))*diff(Rj,t))*((diff(Rj,t))^2) == 1/rho *(psj + 1/c*diff(Rj*psj,t)) - (1/Dij)*(2*Ri*diff(Ri,t)^2 + Ri^2*diff(Ri,t,2) );
ode = [ode1; ode2];
[V,S] = odeToVectorField(ode);
M = matlabFunction(V,'vars',{'t','Y'});
y0 = [Roj 0 Roi 0];
tspan = [0 5];
[t,y] = ode15s(M,tspan,y0);
plot(t,[y(:,1),y(:,3)])
grid on

[Shishir Raut]

Thank you, I greatly appreciate it!