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

Python | 数学 | 物理 | 离散化 | 偏微分方程 | 常微分方程 | 三维 | 热方程 | 资产 | 空间导数 | 地震波 | 微磁学 | 磁倾斜导数 | 西塔规则 | 算法

Python微磁学磁倾斜和西塔规则算法 | Notion亚图跨际 on Notion

🏈指点迷津 | Brief

📜有限差分-用例

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

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

📜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}{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$$

让我们考虑另一个偏微分方程，

$$u_t=b u_{x x} ; \quad b>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$$

示例：数值求解

$$u_t=0.05 u_{x x}$$

• $u$ 代表温度

• $x$表示 $0 \leq x \leq L$​ 的位置

• $t$表示$t>0$的时间

• 边界条件为 $u(t, 0)=0$ 和 $u(t, L)=0（t>0）$

• 初始条件为 $u(0, x)=\sin (\pi x)$对于 $0 \leq x \leq 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 \leq x \leq L$

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

代码数值解：

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