scipy一维热传导方程的求解

1.源码实现

from scipy.integrate import odeint
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
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 np.arange(0,M+1):
        U[0,k] = 0.0
        U[N,k] = 0.0

#初始条件
for i in np.arange(0,N):
        U[i,0]=4*i*h*(3-i*h)

#递推关系
for k in np.arange(0,M):
        for i in np.arange(1,N):
                U[i,k+1]=A*U[i+1,k]+(1-2*A)*U[i,k]+A*U[i-1,k]

def defun(i, k):
        return U[i,k]

x = np.arange(0,N)
t = np.arange(0,M)

x,t = np.meshgrid(x, t)
s = defun(x, t)

x = x * h
t = t * dt

#print(x)
#print(t)
#print(s)

# 绘图
fig = plt.figure()
ax = Axes3D(fig)
ax.invert_xaxis()
ax.plot_surface(x, t, s, cmap="viridis")

plt.savefig("test.png")

2.运行程序

$ python3 example.py

3.执行结果

test.png