Having trouble solving two second order ODEs simultaneously

Question

Shishir Raut 2023-2-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1912160-having-trouble-solving-two-second-order-odes-simultaneously

评论： Shishir Raut 2023-2-16

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-2-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1912160-having-trouble-solving-two-second-order-odes-simultaneously#answer_1171260

在 MATLAB Online 中打开

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 )

10 个评论
显示 8更早的评论隐藏 8更早的评论

Shishir Raut 2023-2-16

编辑：Shishir Raut 2023-2-16

在 MATLAB Online 中打开

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.

Shishir Raut 2023-2-16

在 MATLAB Online 中打开

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 2023-2-16

编辑：Torsten 2023-2-16

在 MATLAB Online 中打开

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 2023-2-16

Thank you, I greatly appreciate it!

请先登录，再进行评论。

Having trouble solving two second order ODEs simultaneously

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Having trouble solving two second order ODEs simultaneously

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论