🎯要点
🎯二维热传导二阶偏微分方程 | 🎯调和函数和几何图曲率 | 🎯解潮汐波动方程 | 🎯解静止基态旋转球体流体运动函数 | 🎯水文空间插值 | 🎯流体流动模拟求解器 | 🎯随机算法解二维高压电容器电势 | 🎯解空心导电圆柱体交替电势 | 🎯稠密矩阵椭圆微分快速算法
📜拉普拉斯方程用例:Python火焰锋动力学和浅水表面波浪偏微分方程
📜泊松方程用例:Python低溫半导体电子束量子波算法计算
🍇Python椭圆微分-拉普拉斯方程
物理学中的许多问题与时间无关,但却具有丰富的物理意义:大质量物体产生的引力场、电荷分布的电势、拉伸膜的位移以及流体通过多孔介质的稳定流动……所有这些都可以用泊松方程建模:
∇
2
u
=
f
\nabla^2 u=f
∇2u=f
其中未知的
u
u
u 和已知的
f
f
f 是域
Ω
\Omega
Ω 中的空间函数。为了找到解,我们需要边界条件。边值问题包括在给定上述信息的情况下找到
u
u
u。在数字上,我们可以使用松弛方法来做到这一点,该方法从对
u
u
u 的初始猜测开始,然后迭代求解。
f
=
0
f=0
f=0(齐次情况)的特殊情况得出拉普拉斯方程:
∇
2
u
=
0
\nabla^2 u=0
∇2u=0
例如,稳定的二维热传导方程为:
∂
2
T
∂
x
2
+
∂
2
T
∂
y
2
=
0
\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=0
∂x2∂2T+∂y2∂2T=0
其中
T
T
T 是已达到稳定状态的温度。拉普拉斯方程对系统在所提供的边界条件下的平衡状态进行建模。研究拉普拉斯方程解的学科称为势理论,解本身通常就是势场。从现在开始,我们用
p
p
p 来表示我们的通用因变量,并再次写出拉普拉斯方程(二维):
∂
2
p
∂
x
2
+
∂
2
p
∂
y
2
=
0
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2}=0
∂x2∂2p+∂y2∂2p=0
与扩散方程一样,我们用中心差离散化二阶导数
p
i
+
1
,
j
−
2
p
i
,
j
+
p
i
−
1
,
j
Δ
x
2
+
p
i
,
j
+
1
−
2
p
i
,
j
+
p
i
,
j
−
1
Δ
y
2
=
0
\frac{p_{i+1, j}-2 p_{i, j}+p_{i-1, j}}{\Delta x^2}+\frac{p_{i, j+1}-2 p_{i, j}+p_{i, j-1}}{\Delta y^2}=0
Δx2pi+1,j−2pi,j+pi−1,j+Δy2pi,j+1−2pi,j+pi,j−1=0
当
Δ
x
=
Δ
y
\Delta x=\Delta y
Δx=Δy 时,我们最终得到以下等式:
p
i
+
1
,
j
+
p
i
−
1
,
j
+
p
i
,
j
+
1
+
p
i
,
j
−
1
−
4
p
i
,
j
=
0
p_{i+1, j}+p_{i-1, j}+p_{i, j+1}+p_{i, j-1}-4 p_{i, j}=0
pi+1,j+pi−1,j+pi,j+1+pi,j−1−4pi,j=0
这告诉我们,网格点
(
i
,
j
)
(i, j)
(i,j) 处的拉普拉斯微分算子可以使用该点处的
p
p
p 值(因子 -4 )和左右四个相邻点来离散计算,网格点
(
i
,
j
)
(i, j)
(i,j) 上方和下方。
假设我们想在一块计算机芯片上模拟稳态热传递,该芯片一侧绝缘(零诺伊曼边界层),两侧保持固定温度(狄利克雷条件),一侧接触具有正弦温度分布的组件。我们需要求解拉普拉斯方程,其边界条件如下:
p
=
0
at
x
=
0
∂
p
∂
x
=
0
at
x
=
L
p
=
0
at
y
=
0
p
=
sin
(
3
2
π
x
L
)
at
y
=
H
.
\begin{gathered} p=0 \text { at } x=0 \\ \frac{\partial p}{\partial x}=0 \text { at } x=L \\ p=0 \text { at } y=0 \\ p=\sin \left(\frac{\frac{3}{2} \pi x}{L}\right) \text { at } y=H . \end{gathered}
p=0 at x=0∂x∂p=0 at x=Lp=0 at y=0p=sin(L23πx) at y=H.
我们将
L
=
1
L=1
L=1 和
H
=
1
H=1
H=1 作为域在
x
x
x 和
y
y
y 方向上的大小。
from matplotlib import pyplot
import numpy
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
def plot_3D(x, y, p):
fig = pyplot.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
X,Y = numpy.meshgrid(x,y)
surf = ax.plot_surface(X,Y,p[:], rstride=1, cstride=1, cmap=cm.viridis,
linewidth=0, antialiased=False)
ax.set_xlim(0,1)
ax.set_ylim(0,1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
ax.view_init(30,45)
具有上述边界条件的拉普拉斯方程有一个解析解,由下式给出
p
(
x
,
y
)
=
sinh
(
3
2
π
y
L
)
sinh
(
3
2
π
H
L
)
sin
(
3
2
π
x
L
)
p(x, y)=\frac{\sinh \left(\frac{\frac{3}{2} \pi y}{L}\right)}{\sinh \left(\frac{\frac{3}{2} \pi H}{L}\right)} \sin \left(\frac{\frac{3}{2} \pi x}{L}\right)
p(x,y)=sinh(L23πH)sinh(L23πy)sin(L23πx)
def p_analytical(x, y):
X, Y = numpy.meshgrid(x,y)
p_an = numpy.sinh(1.5*numpy.pi*Y / x[-1]) /\
(numpy.sinh(1.5*numpy.pi*y[-1]/x[-1]))*numpy.sin(1.5*numpy.pi*X/x[-1])
return p_an
让我们尝试一下解析解,并用它来测试我们上面编写的plot_3D
函数。
nx = 41
ny = 41
x = numpy.linspace(0,1,nx)
y = numpy.linspace(0,1,ny)
p_an = p_analytical(x,y)
plot_3D(x,y,p_an)