如何实现两个耦合二阶微分方程组的六个初始条件？

Question

我对编码很陌生。关于我的问题的基本信息是: r1和r2是两个变量: u1 = dr1/dt，u2 = dr2/dt，和du1/dt = d^2r1/dt^2，du2/dt = d^2r2/dt^2。在Matlab代码中，r(1)表示r1，r(2) -> u1，r(3) -> r2，r(4) -> u2。rdot(2)是du1/dt的表达式，rdot(4)是du2/dt的表达式，理想情况下我只需要4个初始条件: r1( 0 )，u1(0)，r2(0)，u2(0)，它们是10d-6，0，5d-6，0。但在我的例子中，du1/dt依赖于du2/dt，反之亦然。见T1_1和T2_1的最后一项，而du1/dt和du2/dt的理想IC都是0。但是如何在我的代码中实现这一点呢？

我的密码在这里。

function rdot = f(t, r)

P_stat = 1.01325d5;
P_v = 2.3388d3;
mu = 1.002d-3;
sigma = 72.8d-3;
c_s = 1481d0;
poly_exp = 1.4d0;
rho = 998.2071d0;
f_s = 20d3;
P_s = 1.01325d5;
r1_eq = 10d-6;
r2_eq = 4d-6;
d = 1d-3;


rdot(1) = r(2);
P1_bw = ( (P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(1)) - (4.d0*mu*r(2)/r(1));
P1_ext = P_s*sin(2.d0*pi*f_s*(t + (r(1)/c_s)));
T2_1 = ((2.d0*r(3)*(r(4)^2.d0)) + ((r(3)^2.d0)*rdot(4)))/d;
T2_4 = (1.d0 - (r(2)/c_s))*r(1);
T2_5 = 1.5d0*(1.d0 - (r(2)/(3.d0*c_s)))*(r(2)^2.d0);
T2_6 = (1.d0 + (r(2)/c_s))*(P1_bw - P_stat + P_v - P1_ext)/rho;
T2_8 = ( (-3.d0*poly_exp*r(2)*(P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(2)/r(1)) - (4.d0*mu*(- ((r(2)^2.d0)/r(1)))) )/r(1);
T2_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(1)/c_s))))*(1.d0 + (r(2)/c_s) );
T2_7 = (r(1)/(rho*c_s))*(T2_8 - T2_9);
rdot(2) =  (T2_6 + T2_7 - T2_1 - T2_5)/(T2_4 + (4.d0*mu/(rho*c_s)));


rdot(3) = r(4);
P2_bw = ( (P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(3))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(3)) - (4.d0*mu*r(4)/r(3));
P2_ext = P_s*sin(2.d0*pi*f_s*(t + (r(3)/c_s)));
T1_1 = ((2.d0*r(1)*(r(2)^2.d0)) + ((r(1)^2.d0)*rdot(2)))/d;
T1_4 = (1.d0 - (r(4)/c_s))*r(3);
T1_5 = 1.5d0*(1.d0 - (r(4)/(3.d0*c_s)))*(r(4)^2.d0);
T1_6 = (1.d0 + (r(4)/c_s))*(P2_bw - P_stat + P_v - P2_ext)/rho;
T1_8 = ( (-3.d0*poly_exp*r(4)*(P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(3))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(4)/r(3)) - (4.d0*mu*(- ((r(4)^2.d0)/r(3)))) )/r(3);
T1_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(3)/c_s))))*(1.d0 + (r(4)/c_s) );
T1_7 = (r(3)/(rho*c_s))*(T1_8 - T1_9);
rdot(4) =  (T1_6 + T1_7 - T1_1 - T1_5)/(T1_4 + (4.d0*mu/(rho*c_s)));

rdot = rdot';

clc;
clear all;
close all;
time_range = [0 3000d-6];
initial_conditions = [10d-6 0.d0 5d-6 0.d0];
[t, r] = ode45('bubble', time_range, initial_conditions);
plot(t, r(:, 1), t, r(:, 3));

Answer 1

对于每一级的ODE，您都需要一个初始条件。这是由于您正在计算的函数的数量。在您的示例中，您得到了一个二阶ODE系统，该系统在求解后将提供以下四个功能：r(1), r(2), r(3), r(4)

为什么我们甚至需要初始条件？假设您有一个简单的导数：y'= y，我们知道函数y = exp(x) * C解决了这个问题，但是我们需要调整C才能得到“唯一的解决方案”。另一方面，给y'一个初始条件是没有意义的，因为一旦定义了y，它就被完全定义了。无论“外来”变量是以导数的形式出现还是以线性因子的形式出现，都不重要。它与集成电路的数量无关。

我希望我能澄清一点。在我看来，你的程序应该是那样的，但我还没有机会试一试。

Answer 2

从本质上说你的处境是

rdot(2) = a2 + b2*rdot(4)
rdot(4) = a4 + b4*rdot(2)

其中a2,b2,a4,b4包含表达式中的所有其他术语。

这是一个线性系统，你必须解决它才能得到正确的值来返回。您可以使用Matlab的线性求解器，或者在这个简单的二维情况下手工完成，

rdot(2) = a2 + b2*(a4 + b4*rdot(2))  ==>  rdot(2) = (a2 + b2*a4) / (1 - b2*b4)
rdot(4) = a4 + b4*(a2 + b2*rdot(4))  ==>  rdot(2) = (a4 + b4*a2) / (1 - b2*b4)

要应用这一点，您需要拆分

T1_1  as T1_1a + T1_1b*rdot(4), 

您可以根据给定的常量和状态变量计算T1_1a和T1_1b。最后你有了什么

 rdot(2)=(other + coeff*T1_1)/denom 

你得分头进去

 a2 =(other + coeff*T1_1a) / denom    and 
 b2 =         coeff*T1_1b  / denom 

然后对第二部分进行同样的处理，得到a4,b4，然后应用上面的解公式。

Answer 3

我已经解决了这个问题，通过引入两个带有初始条件的新变量，但是它们在代码运行时被更新。新代码包含两个变量: r2ddot和r1ddot；

function rdot = f(t, r)

P_stat = 1.01325d5;
P_v = 2.3388d3;
mu = 1.002d-3;
sigma = 72.8d-3;
c_s = 1481d0;
poly_exp = 1.4d0;
rho = 998.2071d0;
f_s = 20d3;
P_s = 1.01325d5;
r1_eq = 4d-6;
r2_eq = 5d-6;
d = 50*(r1_eq + r2_eq);

r2ddot = 0;
r1ddot = 0;

rdot(1) = r(2);
P2_bw = ( (P_stat - P_v + (2.d0*sigma/r2_eq))*((r2_eq/r(3))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(3)) - (4.d0*mu*r(4)/r(3));
P2_ext = P_s*sin(2.d0*pi*f_s*(t + (r(3)/c_s)));
T1_1 = ((2.d0*r(1)*(r(2)^2.d0)) + ((r(1)^2.d0)*r1ddot))/d;
T1_4 = (1.d0 - (r(4)/c_s))*r(3);
T1_5 = 1.5d0*(1.d0 - (r(4)/(3.d0*c_s)))*(r(4)^2.d0);
T1_6 = (1.d0 + (r(4)/c_s))*(P2_bw - P_stat + P_v - P2_ext)/rho;
T1_8 = ( (-3.d0*poly_exp*r(4)*(P_stat - P_v + (2.d0*sigma/r2_eq))*((r2_eq/r(3))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(4)/r(3)) - (4.d0*mu*(- ((r(4)^2.d0)/r(3)))) )/r(3);
T1_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(3)/c_s))))*(1.d0 + (r(4)/c_s) );
T1_7 = (r(3)/(rho*c_s))*(T1_8 - T1_9);
rdot(2) =  (T1_6 + T1_7 - T1_1 - T1_5)/(T1_4 + (4.d0*mu/(rho*c_s))) ;
r2ddot = rdot(2);

rdot(3) = r(4);
P1_bw = ( (P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(1)) - (4.d0*mu*r(2)/r(1));
P1_ext = P_s*sin(2.d0*pi*f_s*(t + (r(1)/c_s)));
T2_1 = ((2.d0*r(3)*(r(4)^2.d0)) + ((r(3)^2.d0)*r2ddot))/d;
T2_4 = (1.d0 - (r(2)/c_s))*r(1);
T2_5 = 1.5d0*(1.d0 - (r(2)/(3.d0*c_s)))*(r(2)^2.d0);
T2_6 = (1.d0 + (r(2)/c_s))*(P1_bw - P_stat + P_v - P1_ext)/rho;
T2_8 = ( (-3.d0*poly_exp*r(2)*(P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(2)/r(1)) - (4.d0*mu*(- ((r(2)^2.d0)/r(1)))) )/r(1);
T2_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(1)/c_s))))*(1.d0 + (r(2)/c_s) );
T2_7 = (r(1)/(rho*c_s))*(T2_8 - T2_9);
rdot(4) =  (T2_6 + T2_7 - T2_1 - T2_5)/(T2_4 + (4.d0*mu/(rho*c_s)));
r1ddot = rdot(4);
rdot = rdot';


clc;
clear all;
close all;
time_range = [0 1d-3];
initial_conditions = [4d-6 0.d0 5d-6 0.d0];
[t, r] = ode45('bubble', time_range, initial_conditions);
plot(t, r(:, 1), t, r(:, 3));

但我无法得到预期的结果。实际上，我正在尝试复制所附论文的结果，参见等式7。在这里输入链接描述。

​

​

第二组方程可通过互换性指标1和2得到。

重要注意事项:Mettin的论文中方程7的最后一项中有一个错误，可以通过检查各个项的尺寸来验证。正确的最后一项可以从另一篇论文中看出，https://journals.aps.org/pre/abstract/10.1103/PhysRevE.83.066313在下面的方程(1)中看到了最后一项。忽略方程中的其他附加项。重要的一点是方程在Rj中是二阶的，方程在Ri中有一个二阶项。这就是我想要的代码。

​

任何帮助都将不胜感激。