计算科学 - ODE45：对结果的怀疑。正确与否？ - 吾爱随笔录

ODE45：对结果的怀疑。正确与否？

计算科学 matlab 颂

2021-12-08 22:20:03

我必须解决我不知道它是否有帮助，但请注意非常大，而时间跨度非常小（）。我得到的结果没有 MATLAB 返回任何错误或警告，但这个结果并不能说服我（我预计峰值在左右）。我已经尝试切换到 ode15s，但输出是一样的。我还能做些什么来查看输出是否发生变化，然后确保这个结果是正确的？

{\begin{matrix} x \ddot{x} + \frac{3}{2} {\dot{x}}^{2} = \frac{1}{ρ} (P_{v} - a_{1} t^{2} - a_{2} t - P_{0}) \\ x (0) = 5 \cdot 10^{- 6} \\ \dot{x} (0) = 0 \end{matrix}

$\left\{\begin{matrix} x\ddot{x}+\frac{3}{2}\dot{x}^2=\frac{1}{\rho}\left (P_v-a_1t^2-a_2t-P_0 \right )\\ x(0)=5\cdot10^{-6}\\ \dot{x}(0)=0 \end{matrix}\right.$

a_{1} = 2.488125 \cdot 10^{14}

$a_1=2.488125\cdot 10^{14}$

t_{f} = 70 \cdot 10^{- 6}

$t_f=70\cdot 10^{-6}$

t = 4 \cdot 10^{- 5}

$t=4\cdot 10^{-5}$

function vdot = no_thermal_effect_98(t,v)
rho = 959.78;
P0 = 101325;
Pvap = 94285.313;
a1 = (P0 - 1800) / (20e-6)^2;    %nondimensional
a2 = 2 * (1800 - P0) / (20e-6);  %nondimensional
%                                       
vdot = zeros(2,1);
vdot(1) = v(2);
vdot(2) = -1.5 * v(2) * v(2) / v(1) + 1 / (v(1) * rho) *...
    (Pvap - P0 - a1 * t^2 - a2 * t);

跑步

x0 = 5e-6; %meters
tf = 70e-6; %seconds
[t,v] = ode45(@no_thermal_effect_98,[0,tf],[x0,0]);
[t,v(:,1)];
plot(t,v(:,1))

编辑

我们首先考虑气泡外压力恒定的情况：。 $x\ddot{x}+\frac{3}{2}\dot{x}^2=\frac{1}{\rho}\left (P_v-P_0\right )$

如果，则气泡破裂（即半径单调减小）；相反，如果，则气泡会增长。 $P_0>P_v$ $P_0<P_v$

描述的气泡外的压力变化。 $P(t)=a_1t^2+a_2t+P_0$

在时，气泡外的压力大于气泡内的压力（在图中，用橙色线表示，等于）。这就是为什么我预计峰值，即泡沫开始崩溃的那一刻，将是。 $40 \mu s$ $Pv=94285.313 Pa$ $40 \mu s$

1个回答

我在 Python 中运行了相同的代码

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt


def func(v, t):
rho = 959.78
P0 = 101325
Pvap = 94285.313
a1 = (P0 - 1800) / (20e-6)**2
a2 = 2 * (1800 - P0) / (20e-6) 
v1 = v[1];
v2 = -1.5 * v[1] * v[1] / v[0] + 1 / (v[0] * rho) *\
(Pvap - P0 - a1 * t**2 - a2 * t)
return [v1, v2]


y0 = 5e-6
tf = 70e-6
t = np.linspace(0, tf, 1000)
y = odeint(func, [y0, 0], t)
plt.plot(1e5*t, 1e5*y[:, 0])

我得到了同样的结果

如果您担心常量的数值，也许您可以更改单位或无量纲化您的方程。例如，您可以更改时间 $\tau = t/t_0$ ，和 $t_0 = 20\times 10^{-6}$ . 并使用链式法则对左侧的导数进行更改，以获得类似

y y^{″} + \frac{3}{2} y^{' 2} = \frac{1}{ρ} [(P_{v} - P_{0}) - (P_{0} - P_{1}) τ^{2} - 2 (P_{1} - P_{0}) τ]

$y y''+\frac{3}{2}y'^2 = \frac{1}{\rho}\left [(P_v - P_0) - (P_0 - P_1)\tau^2 - 2(P_1 - P_0)\tau\right]$

其中、和。我尝试了这个版本并给了我相同的答案。 $y = x/t_0$ $y' = dy/d\tau$ $P_1 = 1800$

笔记

在这两种情况下，我都将轴缩放。 $10^5$
使用问题中的一些物理知识对方程进行无量纲化可能会更好。您可以查看参考 1 以获得一些想法。

参考

Kudryashov、Nikolay A. 和 Dmitry I. Sinelshchikov。“空泡和充气泡的瑞利方程的解析解。” 物理学杂志 A：数学与理论 47.40（2014）：405202。

其它你可能感兴趣的问题

上一篇两能级加法 Schwarz 中的全局问题是什么？下一篇如果步长太小，微分方程的数值求解会出错吗？