Python火焰锋动力学和浅水表面波浪偏微分方程

🎯要点

🎯流图可视化正弦余弦矢量场 | 🎯解空间变化边界条件二维拉普拉斯方程 | 🎯解圆柱坐标系标量场 | 🎯解一维泊松方程 | 🎯解二维扩散方程 | 🎯解火焰锋的动力学偏微分方程 | 🎯解球对称几何偏微分方程 | 🎯解笛卡尔网格扩散方程 | 🎯解时间相关边界条件一维扩散方程 | 🎯解浅水表面波浪偏微分方程 | 🎯解空间耦合疾病传播偏微分方程模型 | 🎯化学自催化反应的理论模型 | 🎯解在空间扩散方程 | 🎯解量子力学波函数偏微分方程

📜微分方程 | 本文 - 用例

📜数学模型: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 Ax22u+Bxy2u+Cy22u+Dxu+Eyu+Fu=G
通过仅考虑前三个系数 A B 和 C,我们可以确定我们正在处理什么方程或我们正在解决什么问题。

如果 B 2 − 4 A C < 0 B^2-4 A C<0 B24AC<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 x22u+y22u=0
该方程的两个导数是空间 x 2 x^2 x2 y 2 y^2 y2的导数,没有时间导数。

如果 B 2 − 4 A C > 0 B^2-4 A C>0 B24AC>0,则我们有双曲偏微分方程,为波动方程:
∂ 2 u ∂ t 2 = ∂ 2 u ∂ y 2 \frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial y^2} t22u=y22u
如果 B 2 − 4 A C = 0 B^2-4 A C=0 B24AC=0,则我们有抛物线偏微分方程,为扩散方程:
∂ u ∂ t = ∂ 2 u ∂ y 2 \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial y^2} tu=y22u
尽管计算机非常擅长数学,但它们不懂微分方程。为了告诉计算机求解微分方程,我们需要对方程进行离散化。即:
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)+Δxxu+2(Δx)2x22u+
通过丢弃二阶项,我们最终将得到一个非常近似的公式:
∂ 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} xuΔ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) (xu)=h1[ui+1,jui,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) (xu)=h1[ui,jui1,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)=xxu+2(Δx)3x33u+

∂ 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} xuxu(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) (xu)i,j=2h1[ui+1,jui1,j]+O(h2)
一维扩散计算代码:

import numpy as np
from matplotlib import pyplot as plt

L = 100

x = np.linspace(0, 1, L)
u = np.zeros(L)


u[L//2] = 1  
u[0] = 0
u[L-1] = 0

for t in range(100):
    for i in range(1, L-1):
        u[i] += (u[i+1] - 2*u[i] + u[i-1])/4

fig, 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] = 1

fig, ax = plt.subplots(figsize=(20,10))
ax.axis('off')
plot = ax.contourf(u, cmap='jet')
def ans(f):
    global u, plot
    
    for 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])/4
    for c in plot.collections:
        c.remove()
    plot = ax.contourf(u, cmap='jet')
    return plot

anim = animation.FuncAnimation(fig, ans, frames=T)
plt.show()

👉参阅:计算思维 | 亚图跨际

### 使用Python求解二浅水方程 为了使用Python求解二浅水方程,可以借鉴有限差分方法来离散化微分方程,并利用`numpy`库来进行高效的矩阵运算。对于可视化方面,则可借助`matplotlib`及其扩展模块`mpl_toolkits.mplot3d`展示三形。 #### 定义基本变量与参数 首先定义网格大小、时间步长dt、空间步dxdy等基础配置项;接着初始化位高度\( \zeta \),平方向上的速u以及垂直方向上的速v这三个主要物理量的分布情况[^2]。 ```python import numpy as np from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt # 设置计算域尺寸及分辨率 nx, ny = 100, 50 # X轴Y轴格点数 Lx, Ly = 200., 100. # 计算区域范围(m) dx, dy = Lx/(nx-1), Ly/(ny-1) # 时间相关参数 nt = 100 # 总迭代次数 dt = .01 # 时间间隔(s) g = 9.8 # 重力加速度(m/s²) # 初始化场变量 ζ = np.zeros((ny,nx)) u = np.zeros((ny,nx)) v = np.zeros((ny,nx)) # 边界条件其他初始状态设定... ``` #### 更新规则设计 基于质量守恒定律动量守恒定律推导出来的更新公式用于逐步推进系统的演化过程: \[ \begin{aligned} &\frac{\partial h}{\partial t}+\nabla\cdot(h\vec v)=0 \\ &\frac{\partial (hu)}{\partial t}+f_{cor}\times hv-\nu\nabla^{2}(hu)+gh\nabla H=0\\ &\frac{\partial (hv)}{\partial t}-f_{cor}\times hu-\nu\nabla^{2}(hv)+gh\nabla H=0 \end{aligned} \] 其中\( f_{cor}=2Ωsinφ \)代表科里奥利效应的影响因子,而ν表示粘滞系数。上述表达式可以通过显式欧拉法或者其他更高级的时间积分算法转化为易于编程的形式[^4]。 具体到代码实现时,考虑到稳定性需求通常会选择较为保守的空间离散方案比如迎风格式或者中心差商形式。下面给出了种简化版的阶精度前向Euler方法作为例子说明如何编核心循环逻辑: ```python for n in range(nt): ζ_new = ... # 根据当前时刻的状态预测下刻的高度变化 u_new = ... # 同样地更新其他两个未知数 v_new = ... # 应用周期性边界条件或其他特定场景下的约束 # 将新值赋给旧值准备下轮迭代 ζ[:], u[:], v[:] = ζ_new.copy(), u_new.copy(), v_new.copy() ``` #### 可视化结果呈现 最后通过动画技术动态显示随时间演化的形态有助于直观理解整个物理过程。这步骤同样依赖于Matplotlib的强大绘功能完成: ```python fig = plt.figure(figsize=(8,6)) ax = fig.gca(projection='3d') def animate(i): ax.clear() surf = ax.plot_surface(X,Y,ζ,cmap="coolwarm",linewidth=0,antialiased=False) return [surf] ani = animation.FuncAnimation(fig,animate,interval=20,blit=True) plt.show() ``` 以上即为构建一个简单版本的二浅水方程数值模拟器所需的关键要素概述。当然实际应用中还需要考虑更多细节优化如提高时空分辨能力、引入复杂地形影响因素等等。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值