🎯要点
🎯流图可视化正弦余弦矢量场 | 🎯解空间变化边界条件二维拉普拉斯方程 | 🎯解圆柱坐标系标量场 | 🎯解一维泊松方程 | 🎯解二维扩散方程 | 🎯解火焰锋的动力学偏微分方程 | 🎯解球对称几何偏微分方程 | 🎯解笛卡尔网格扩散方程 | 🎯解时间相关边界条件一维扩散方程 | 🎯解浅水表面波浪偏微分方程 | 🎯解空间耦合疾病传播偏微分方程模型 | 🎯化学自催化反应的理论模型 | 🎯解在空间扩散方程 | 🎯解量子力学波函数偏微分方程
📜微分方程 | 本文 - 用例
📜数学模型:Python流感常微分方程房室数学模型
📜图计算和算法:Python流感传播感染康复图模型计算和算法
📜水文模型:Python流体数据统计模型和浅水渗流平流模型模拟
🍇Python三个维度扩散方程
偏微分方程是具有多个独立变量、依赖于这些变量的未知函数以及未知函数关于独立变量的偏导数的方程。
通常样式为:
A ∂ 2 u ∂ x 2 + B ∂ 2 u ∂ x ∂ y + C ∂ 2 u ∂ y 2 + D ∂ u ∂ x + E ∂ u ∂ y + F u = G A \frac{\partial^2 u}{\partial x^2}+ B \frac{\partial^2 u}{\partial x \partial y}+ C \frac{\partial^2 u}{\partial y^2}+ D \frac{\partial u}{\partial x}+ E \frac{\partial u}{\partial y}+ F u= G A∂x2∂2u+B∂x∂y∂2u+C∂y2∂2u+D∂x∂u+E∂y∂u+Fu=G
通过仅考虑前三个系数 A B 和 C,我们可以确定我们正在处理什么方程或我们正在解决什么问题。
如果 B 2 − 4 A C < 0 B^2-4 A C<0 B2−4AC<0,我们有一个椭圆偏微分方程,为拉普拉斯方程:
∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = 0 \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0 ∂x2∂2u+∂y2∂2u=0
该方程的两个导数是空间 x 2 x^2 x2和 y 2 y^2 y2的导数,没有时间导数。
如果 B 2 − 4 A C > 0 B^2-4 A C>0 B2−4AC>0,则我们有双曲偏微分方程,为波动方程:
∂ 2 u ∂ t 2 = ∂ 2 u ∂ y 2 \frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial y^2} ∂t2∂2u=∂y2∂2u
如果 B 2 − 4 A C = 0 B^2-4 A C=0 B2−4AC=0,则我们有抛物线偏微分方程,为扩散方程:
∂ u ∂ t = ∂ 2 u ∂ y 2 \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial y^2} ∂t∂u=∂y2∂2u
尽管计算机非常擅长数学,但它们不懂微分方程。为了告诉计算机求解微分方程,我们需要对方程进行离散化。即:
U = u ( x , y ) U = u ( x , y ) U=u(x,y)
以 x x x 表示的部分展开式为
u ( x + Δ x , y ) = u ( x , y ) + Δ x ∂ u ∂ x + ( Δ x ) 2 2 ∂ 2 u ∂ x 2 + … u(x+\Delta x, y)=u(x, y)+\Delta x \frac{\partial u}{\partial x}+\frac{(\Delta x)^2}{2} \frac{\partial^2 u}{\partial x^2}+\ldots u(x+Δx,y)=u(x,y)+Δx∂x∂u+2(Δx)2∂x2∂2u+…
通过丢弃二阶项,我们最终将得到一个非常近似的公式:
∂ u ∂ x ≈ u ( x + Δ x , y ) − u ( x , y ) Δ x \frac{\partial u}{\partial x} \approx \frac{u(x+\Delta x, y)-u(x, y)}{\Delta x} ∂x∂u≈Δxu(x+Δx,y)−u(x,y)
这是一个著名的欧拉方法,在常微分方程中可以看到。在某些有限差分中,它被称为前向差分。
📜流体力学电磁学运动学动力学化学和电路中欧拉法示例:Python重力弹弓流体晃动微分方程模型和交直流电阻电容电路
看一下下面的例子:
h = Δ x u i , j = u ( x i , y j ) u i + 1 , j = u ( x i + h , y j ) \begin{aligned} & h=\Delta x \\ & u_{i, j}=u\left(x_i, y_j\right) \\ & u_{i+1, j}=u\left(x_i+h, y_j\right) \end{aligned} h=Δxui,j=u(xi,yj)ui+1,j=u(xi+h,yj)
正向差分是:
( ∂ u ∂ x ) = 1 h [ u i + 1 , j − u i , j ] + O ( h ) \left(\frac{\partial u}{\partial x}\right)=\frac{1}{h}\left[u_{i+1, j}-u_{i, j}\right]+O(h) (∂x∂u)=h1[ui+1,j−ui,j]+O(h)
通过查看正向差分公式,我们注意到正向差分是通过当前步骤 u i , j u _{ i , j } ui,j 减去下一步 u i + 1 , j u _{ i +1, j} ui+1,j 得出的。直观上,我们可以想到反向差分。
( ∂ u ∂ x ) = 1 h [ u i , j − u i − 1 , j ] + O ( h ) \left(\frac{\partial u}{\partial x}\right)=\frac{1}{h}\left[u_{i, j}-u_{i-1, j}\right]+O(h) (∂x∂u)=h1[ui,j−ui−1,j]+O(h)
中间差分:
u ( x + Δ x , y ) − u ( x − Δ x , y ) = 2 Δ x ∂ u ∂ x + ( Δ x ) 3 2 ∂ 3 u ∂ x 3 + … u(x+\Delta x, y)-u(x-\Delta x, y)=2 \Delta x \frac{\partial u}{\partial x}+\frac{(\Delta x)^3}{2} \frac{\partial^3 u}{\partial x^3}+\ldots u(x+Δx,y)−u(x−Δx,y)=2Δx∂x∂u+2(Δx)3∂x3∂3u+…
∂ u ∂ x ≈ u ( x + Δ x , y ) − u ( x − Δ x , y ) 2 Δ x \frac{\partial u}{\partial x} \approx \frac{u(x+\Delta x, y)-u(x-\Delta x, y)}{2 \Delta x} ∂x∂u≈2Δxu(x+Δx,y)−u(x−Δx,y)
最终得到公式:
( ∂ u ∂ x ) i , j = 1 2 h [ u i + 1 , j − u i − 1 , j ] + O ( h 2 ) \left(\frac{\partial u}{\partial x}\right)_{i, j}=\frac{1}{2 h}\left[u_{i+1, j}-u_{i-1, j}\right]+O\left(h^2\right) (∂x∂u)i,j=2h1[ui+1,j−ui−1,j]+O(h2)
一维扩散计算代码:
import numpy as np
from matplotlib import pyplot as pltL = 100x = np.linspace(0, 1, L)
u = np.zeros(L)u[L//2] = 1
u[0] = 0
u[L-1] = 0for t in range(100):for i in range(1, L-1):u[i] += (u[i+1] - 2*u[i] + u[i-1])/4fig, ax = plt.subplots(figsize=(20,10))
ax.axis('off')
ax.scatter(x,u, linewidth=15, c=u, cmap='jet')
plt.show()
二维扩散计算代码:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
import seaborn as sns
sns.set()T = 400
u = np.zeros((100,100))
u[30,20] = 1fig, ax = plt.subplots(figsize=(20,10))
ax.axis('off')
plot = ax.contourf(u, cmap='jet')
def ans(f):global u, plotfor j in range(100):for i in range(1,99):u[i,j] += (u[i+1,j] + u[i,j+1] + u[i-1,j] + u[i,j-1] - 4*u[i,j])/4for c in plot.collections:c.remove()plot = ax.contourf(u, cmap='jet')return plotanim = animation.FuncAnimation(fig, ans, frames=T)
plt.show()