Having trouble solving two second order ODEs simultaneously

2023 2 14
1 个回答
回答已采纳
更新时间：2023 2 16
8 次查看（30 天）

0 个投票

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.

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

 采纳的回答

Torsten
Torsten 2023-2-14

0 个投票

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
Shishir Raut 2023-2-15
Sorry. I corrected it but the error is still there. What should I do?
Torsten
Torsten 2023-2-15
编辑：Torsten 2023-2-15
Ran in:
Ran in:
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
Shishir Raut 2023-2-15
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
Torsten 2023-2-15
编辑：Torsten 2023-2-15
Ran in:
Ran in:
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
Shishir Raut 2023-2-15
Thank you so much! You are a lifesaver!
Shishir Raut
Shishir Raut 2023-2-16
编辑：Shishir Raut 2023-2-16
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
Torsten 2023-2-16
Use
[t,y] = ode15s(M,tspan,y0);
instead of
[t,y] = ode45(M,tspan,y0);
But the result doesn't look promising.
Shishir Raut
Shishir Raut 2023-2-16
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
Torsten 2023-2-16
编辑：Torsten 2023-2-16
Ran in:
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

产品

版本

R2019b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by