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

55 篇文章 0 订阅
28 篇文章 0 订阅

🎯要点

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

📜微分方程 | 本文 - 用例

📜数学模型: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()

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值