利用python求解微分方程

案例詳情：轉載自 (用Python數值求解偏微分方程 - yubr的文章 - 知乎)
[https://zhuanlan.zhihu.com/p/81488678]

import numpy as np
import matplotlib.pyplot as plt

h = 0.1#空間步長
N =30#空間步數
dt = 0.0001#時間步長
M = 10000#時間的步數
A = dt/(h**2) #lambda*tau/h^2
U = np.zeros([N+1,M+1])#建立二維空數組
Space = np.arange(0,(N+1)*h,h)#建立空間等差數列，從0到3，公差是h

#邊界條件
for k in range(0,M+1):
    U[0,k] = 0.0
    U[N,k] = 0.0

#初始條件
for i in range(0,N):
    U[i,0]=4*i*h*(3-i*h)

#遞推關系
for k in range(0,M):
    for i in range(1,N):
        U[i,k+1]=A*U[i+1,k]+(1-2*A)*U[i,k]+A*U[i-1,k]

#不同時刻的溫度隨空間坐標的變化
plt.plot(Space,U[:,0], 'g-', label='t=0',linewidth=1.0)
plt.plot(Space,U[:,3000], 'b-', label='t=3/10',linewidth=1.0)
plt.plot(Space,U[:,6000], 'k-', label='t=6/10',linewidth=1.0)
plt.plot(Space,U[:,9000], 'r-', label='t=9/10',linewidth=1.0)
plt.plot(Space,U[:,10000], 'y-', label='t=1',linewidth=1.0)
plt.ylabel('u(x,t)', fontsize=20)
plt.xlabel('x', fontsize=20)
plt.xlim(0,3)
plt.ylim(-2,10)
plt.legend(loc='upper right')
plt.show()

#溫度等高線隨時空坐標的變化，溫度越高，顏色越偏紅
extent = [0,1,0,3]#時間和空間的取值范圍
levels = np.arange(0,10,0.1)#溫度等高線的變化范圍0-10，變化間隔為0.1
plt.contourf(U,levels,origin='lower',extent=extent,cmap='jet')
plt.ylabel('x', fontsize=20)
plt.xlabel('t', fontsize=20)
plt.show()