Python微磁学磁倾斜和西塔规则算法

📜有限差分-用例

📜离散化偏微分方程求解器和模型定型 | 📜三维热传递偏微分方程解 | 📜特定资产期权价值偏微分方程计算 | 📜三维波偏微分方程空间导数计算 | 📜应力-速度公式一阶声波方程模拟二维地震波 | 📜微磁学计算磁化波动求解器、色散关系和能垒的弦法 | 📜磁倾斜导数数据平滑

📜指数衰减：🖊常微分方程数值求解器 | 🖊绘制衰减图 | 🖊绘制（正向欧拉、反向欧拉和克兰克-尼科尔森）西塔规则算法放大因子图 | 🖊泰勒级数展开符号计算三种算法误差 | 🖊模型误差、数据误差、离散化误差和舍入误差 | 🖊求解器泛化

📜Python热涨落流体力学求解算法和英伟达人工智能核评估模型

📜常微分方程用例：​Python机器人动力学和细胞酶常微分方程

✒️Python不同初始条件下热方程

有限差分法是获得偏微分和代数方程数值解的技术之一。在该方法中，解在有限网格点中以离散形式近似。

首先考虑一个偏微分方程：
$$u_t+a{u}_x=0$$
正向时间前向空间算法由下式给出：
$$\frac{V_m^{n+1}-V_m^n}{k}+a\frac{V_{m+1}^n-V_m^n}{h}=0$$
正向时间中心空间算法由下式给出：
$$\frac{V_m^{n-1}-V_m^n}{k}+a\cdot \frac{V_{m-}^{-}-V_{m-1}^n}{2h}-0$$
中心时间中心空间算法由下式给出
$$\frac{V_m^{n+1}-V_m^{n-1}}{2k}+a\cdot \frac{V_{m+1}^n-V_{m-1}^n}{2h}=0$$
让我们考虑另一个偏微分方程，
$$u_t=b{u}_{xx};b>0$$
正向时间中心空间算法由下式给出：
$$\frac{V_m^{n+1}-V_m^n}{k}=b\frac{V_{m+1}^n-2V_m^n+V_{m-1}^n}{h^2}=0$$
示例：数值求解
$$u_t=0.05u_{xx}$$

• $u$ 代表温度
• $x$ 表示 $0\le x\le L$​ 的位置
• $t$ 表示$t>0$的时间
• 边界条件为 $u(t,0)=0$ 和 $u(t,L)=0$（$t>0$）
• 初始条件为 $u(0,x)=\sin (\pi x)$ 对于 $0\le x\le L$
• $b$ 表示$b>0$ 的扩散系数

代码求解：

import numpy as np
import matplotlib.pyplot as plt


L = 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.C

for 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\le x\le L$

$$u(0,x)=\{\begin{array}{ll}2x& \text{ if }x<0.5\\ 2(1-x)& \text{ 否则 }\end{array}$$

代码数值解：

import numpy as np
import matplotlib.pyplot as plt


L = 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.C

for 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()

👉参阅一：计算思维

👉参阅二：亚图跨际