🍍Python物理学有限差分微分求解器和动画波形传播

关键词

Python | 物理 | 有限差分 | 常微分方程 | 偏微分方程 | 求解器 | 三维 | 二维 | 动画 | 可视化 | 数值 | 符号 | 泰勒序列 | 显性/隐性方法 | 欧拉法 | 离散 | 傅里叶 | 龙格-库塔 | 收敛 | 交错网格 | 线性 | 二次函数 | 阻尼 | 自由体图 | 模拟 | 弹性 | 波形 | 反射 | 网格 | 频率 | 数值色散关系 | 矢量化 | 高斯函数 | 扩散 | 平流 | 非线性 | 数学

Python物理学有限差分微分求解器和动画波形传播 | Notion亞圖跨際 on Notion

🏈page指点迷津 | Brief

🎯要点

1. Python数值和符号计算：

   1. 振动常微分方程：🎯中心差分求解器，绘制移动窗口研究长时间序列。🎯符号计算离散方程量化误差。🎯Python数值对比正向欧拉方法，反向欧拉方法，克兰克-尼科尔森方法和龙格-库塔法解。🎯数值计算振动能量数值收敛。🎯欧拉-克罗默方法求解器。🎯交错网格法求解器。🎯使用线性/二次函数和符号计算验证常微分方程。🎯线性/二次阻尼求解器和验证函数，绘制阻尼正弦函数图。🎯绘制钟摆动态自由体图模拟，创建求解器计算。🎯模拟弹性摆。

   2. 波形运动偏微分方程：🎯弹簧波形模拟：方程算式转换，制定递归算法，代码实现数值求解器，结果动画演示，二次解和敛率的验证函数。示例：创建吉他弦一维波形方程数值求解器，动画模拟弦振动波形。🎯反射边界的诺依曼条件，代码实现一维波形方程，验证函数检测结果。🎯可变波速，代码实现可变系数。🎯代码实现给定条件下，通用一维波形方程求解器。🎯差分方程正弦和余弦波分析，代码实现网格，离散傅里叶变换，计算对应频率。🎯对离散化参数进行泰勒级数展开，代码计算数值色散关系，代码计算二维数值色散关系并绘制二维数值色散误差结果图。🎯代码解二维线性波形方程，矢量化计算函数代码，高斯函数绘制波形在正方形中传播。

   3. 扩散偏微方程：🎯代码实现具有边界条件的一维扩散求解器，代码验证求解器，动画可视化正向欧拉结果。🎯实现前向欧拉算法和稳定条件下一维求解器，实现求解器矢量化函数 🎯一维扩散方程隐性方法：代码实现反向欧拉法求解器，及其矢量化函数。🎯二维和三维扩散方程求解器。🎯异质介质中的扩散离散化后：代码实现求解器，稀疏矩阵处理。🎯分段恒定介质中扩散的代码实现，绘图可视化恒定扩散系数结果。🎯实现稠密系数矩阵求解器，数值验证。🎯代码实现稀疏稀疏矩阵计算。🎯雅可比方法代码实现，解线性系统。

   4. 平流方程，非线性问题

🍇Python热方程有限差分计算示例

热方程基本上是一个偏微分方程，即：

$\frac{\partial u}{\partial t}-\alpha \nabla u=0$

如果我们想用二维（笛卡尔）求解，我们可以这样写上面的热方程

$\frac{\partial u}{\partial t}-\alpha\left(\frac{\partial^2 u}{\partial x}+\frac{\partial^2 u}{\partial y}\right)=0$

其中u是我们想要知道的数量，t是时间变量，x和y是空间变量，$\alpha$是扩散常数。所以基本上我们希望在 x 和 y 中的任何地方以及随着时间的 t 找到解决方案 u。现在让我们简单地了解一下有限差分法。有限差分法是一种通过有限差分逼近导数来求解微分方程的数值方法。请记住导数的定义是：

$f^{\prime}(a)=\lim _{h \rightarrow 0} \frac{f(a+h)-f(a)}{h}$

在有限差分法中，我们对其进行近似并去除限制。 因此，我们不使用微分和极限符号，而是使用 delta 符号，即有限差分。 请注意，这过于简单化了，因为我们必须使用泰勒级数展开，并通过假设某些项足够小来推导出它，但我们得到了该方法背后的粗略想法。

$f^{\prime}(a) \approx \frac{f(a+h)-f(a)}{h}$

在有限差分方法中，我们将“离散化”空间域和时间间隔 x、y 和 t。我们可以这样写

$\begin{aligned} &x_i=i \Delta x\\ &y_j=j \Delta y\\ &t_k=k \Delta t \end{aligned}$

正如我们所看到的，i、j 和 k 分别是 x、y 和 t 的每个差异的步骤。我们想要的是解u，即

$u(x, y, t)=u_{i, j}^k$

请注意，k 是上标，表示 u 的时间步长。我们可以使用有限差分法写出上面的热方程，如下所示

$\frac{u_{i, j}^{k+1}-u_{i, j}^k}{\Delta t}-\alpha\left(\frac{u_{i+1, j}^k-2 u_{i, j}^k+u_{i-1, j}^k}{\Delta x^2}+\frac{u_{i, j+1}^k-2 u_{i, j}^k+u_{i, j-1}^k}{\Delta y^2}\right)=0$

如果我们通过 $\Delta x=\Delta y$ 排列上面的方程，我们得到最终的方程

$u_{i, j}^{k+1}=\gamma\left(u_{i+1, j}^k+u_{i-1, j}^k+u_{i, j+1}^k+u_{i, j-1}^k-4 u_{i, j}^k\right)+u_{i, j}^k$

其中，$\gamma=\alpha \frac{\Delta t}{\Delta x^2}$

对于我们的模型，我们取 $\Delta x=1$ 和$\alpha=2.0$。现在我们可以使用 Python 代码以数值方式解决这个问题，以查看各处的温度（用 i 和 j 表示）和随时间变化的温度（用 k 表示）。我们首先导入所有必要的库，然后设置边界和初始条件。

 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.animation as animation
 from matplotlib.animation import FuncAnimation
 ​
 plate_length = 50
 max_iter_time = 1000
 ​
 alpha = 2.0
 delta_x = 1
 ​
 delta_t = (delta_x ** 2)/(4 * alpha)
 gamma = (alpha * delta_t) / (delta_x ** 2)
 ​
 u = np.empty((max_iter_time, plate_length, plate_length))
 ​
 u_initial = 0.0
 ​
 u_top = 100.0
 u_left = 0.0
 u_bottom = 0.0
 u_right = 0.0
 ​
 u.fill(u_initial)
 ​
 u[:, (plate_length-1):, :] = u_top
 u[:, :, :1] = u_left
 u[:, :1, 1:] = u_bottom
 u[:, :, (plate_length-1):] = u_right
 ​

我们已经设置了初始条件和边界条件，让我们根据上面推导的有限差分方法编写计算函数。

 def calculate(u):
   for k in range(0, max_iter_time-1, 1):
     for i in range(1, plate_length-1, delta_x):
       for j in range(1, plate_length-1, delta_x):
         u[k + 1, i, j] = gamma * (u[k][i+1][j] + u[k][i-1][j] + u[k][i][j+1] + u[k][i][j-1] - 4*u[k][i][j]) + u[k][i][j]
   
   return u

让我们准备绘图函数，以便我们可以将解（对于每个 k）可视化为热图。我们使用Matplotlib库，它很容易使用。