📜有限差分-用例
📜离散化偏微分方程求解器和模型定型 | 📜三维热传递偏微分方程解 | 📜特定资产期权价值偏微分方程计算 | 📜三维波偏微分方程空间导数计算 | 📜应力-速度公式一阶声波方程模拟二维地震波 | 📜微磁学计算磁化波动求解器、色散关系和能垒的弦法 | 📜磁倾斜导数数据平滑
📜指数衰减:🖊常微分方程数值求解器 | 🖊绘制衰减图 | 🖊绘制(正向欧拉、反向欧拉和克兰克-尼科尔森)西塔规则算法放大因子图 | 🖊泰勒级数展开符号计算三种算法误差 | 🖊模型误差、数据误差、离散化误差和舍入误差 | 🖊求解器泛化
📜Python热涨落流体力学求解算法和英伟达人工智能核评估模型
📜常微分方程用例:Python机器人动力学和细胞酶常微分方程
✒️Python不同初始条件下热方程
有限差分法是获得偏微分和代数方程数值解的技术之一。在该方法中,解在有限网格点中以离散形式近似。
首先考虑一个偏微分方程:
u t + a u x = 0 u_t+a u_x=0 ut+aux=0
正向时间前向空间算法由下式给出:
V m n + 1 − V m n k + a V m + 1 n − V m n h = 0 \frac{V_m^{n+1}-V_m^n}{k}+a \frac{V_{m+1}^n-V_m^n}{h}=0 kVmn+1−Vmn+ahVm+1n−Vmn=0
正向时间中心空间算法由下式给出:
V m n − 1 − V m n k + a ⋅ V m − − − V m − 1 n 2 h − 0 \frac{V_m^{n-1}-V_m^n}{k}+a \cdot \frac{V_{m-}^{-}-V_{m-1}^n}{2 h}-0 kVmn−1−Vmn+a⋅2hVm−−−Vm−1n−0
中心时间中心空间算法由下式给出
V m n + 1 − V m n − 1 2 k + a ⋅ V m + 1 n − V m − 1 n 2 h = 0 \frac{V_m^{n+1}-V_m^{n-1}}{2 k}+a \cdot \frac{V_{m+1}^n-V_{m-1}^n}{2 h}=0 2kVmn+1−Vmn−1+a⋅2hVm+1n−Vm−1n=0
让我们考虑另一个偏微分方程,
u t = b u x x ; b > 0 u_t=b u_{x x} ; \quad b>0 ut=buxx;b>0
正向时间中心空间算法由下式给出:
V m n + 1 − V m n k = b V m + 1 n − 2 V m n + V m − 1 n h 2 = 0 \frac{V_m^{n+1}-V_m^n}{k}=b \frac{V_{m+1}^n-2 V_m^n+V_{m-1}^n}{h^2}=0 kVmn+1−Vmn=bh2Vm+1n−2Vmn+Vm−1n=0
示例:数值求解
u t = 0.05 u x x u_t=0.05 u_{x x} ut=0.05uxx
- u u u 代表温度
- x x x 表示 0 ≤ x ≤ L 0 \leq x \leq L 0≤x≤L 的位置
- t t t 表示 t > 0 t>0 t>0的时间
- 边界条件为 u ( t , 0 ) = 0 u(t, 0)=0 u(t,0)=0 和 u ( t , L ) = 0 u(t, L)=0 u(t,L)=0( t > 0 t>0 t>0)
- 初始条件为 u ( 0 , x ) = sin ( π x ) u(0, x)=\sin (\pi x) u(0,x)=sin(πx) 对于 0 ≤ x ≤ L 0 \leq x \leq L 0≤x≤L
- b b b 表示 b > 0 b>0 b>0 的扩散系数
代码求解:
import numpy as np
import matplotlib.pyplot as pltL = 1
T = 1
m = 5
n = 5
h = L / m
k = T / n
b = 0.05
mu = k / h**2 c = b * mu
if c <= 0 or c >= 0.5:print('Scheme is unstable')v = np.zeros((m + 1, n + 1))
ic1 = lambda x: np.sin(np.pi * x)for j in range(1, m + 2):v[0, j - 1] = ic1((j - 1) * h)b1 = lambda t: 0 # L.B.C
b2 = lambda t: 0 # R.B.Cfor i in range(1, n + 2):v[i - 1, 0] = b1((i - 1) * k)v[i - 1, n] = b2((i - 1) * k)for i in range(n):for j in range(1, m):v[i + 1, j] = (1 - 2 * b * mu) * v[i, j] + b * mu * v[i, j + 1] + b * mu * v[i, j - 1]x = np.linspace(0, L, m + 1)
t = np.linspace(0, T, n + 1)
X, T = np.meshgrid(x, t)fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, T, v, cmap='viridis')
ax.set_xlabel('Space X')
ax.set_ylabel('Time T')
ax.set_zlabel('V')
plt.title('Python for Heat')
plt.show()
接上例,初始条件改为:
对于 0 ≤ x ≤ L 0 \leq x \leq L 0≤x≤L
u ( 0 , x ) = { 2 x if x < 0.5 2 ( 1 − x ) 否则 u(0, x)= \begin{cases}2 x & \text { if } x<0.5 \\ 2(1-x) & \text { 否则 }\end{cases} u(0,x)={2x2(1−x) if x<0.5 否则
代码数值解:
import numpy as np
import matplotlib.pyplot as pltL = 1
T = 1
m = 5
n = 5
h = L / m
k = T / n
b = 0.05
mu = k / h ** 2 c = b * mu
if c <= 0 or c >= 0.5:print('Scheme is unstable')v = np.zeros((m + 1, n + 1))
ic1 = lambda x: 2 * x
ic2 = lambda x: 2 * (1 - x)x = np.linspace(0, L, m + 1)
x = np.linspace(0, L, m + 1)
for j in range(1, m + 2):if x[j - 1] < 0.5:v[0, j - 1] = ic1(x[j - 1]) else:v[0, j - 1] = ic2(x[j - 1]) b1 = lambda t: 0 # L.B.C
b2 = lambda t: 0 # R.B.Cfor i in range(1, n + 2):v[i - 1, 0] = b1((i - 1) * k)v[i - 1, n] = b2((i - 1) * k)for i in range(n):for j in range(1, m):v[i + 1, j] = (1 - 2 * b * mu) * v[i, j] + b * mu * v[i, j + 1] + b * mu * v[i, j - 1]# Visualization
x = np.linspace(0, L, m + 1)
t = np.linspace(0, T, n + 1)
X, T = np.meshgrid(x, t)fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, T, v, cmap='viridis')
ax.set_xlabel('Space X')
ax.set_ylabel('Time T')
ax.set_zlabel('V')
plt.title('Python for Heat ')
plt.show()