是這樣的，小弟最近在學(xué)數(shù)值方法，于是參考了以下鏈接的教程：

案例1：

是誰還在胡說有限差分法？_嗶哩嗶哩_bilibili

具體問題如下：

差分后形式如下：

python中，空間步長為0.1，時(shí)間步長為0.1時(shí)，結(jié)果為：

奇怪的來了，按照道理我縮小時(shí)間步長結(jié)果應(yīng)該更精確，但是我取空間步長為0.01，時(shí)間步長為0.01時(shí)，結(jié)果如下：

結(jié)果完全不對(duì)：

代碼如下：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D as ax
beta=0.01
deltax=0.01
deltat=0.01
lx=1
t=1
p=beta*deltat/(deltax*deltax)
N=lx/deltax
M=t/deltat
print(M,N,p)
M=int(M)
N=int(N)
ans=np.zeros((M,N))
def jisuan(x,y,z):
? ?u=(x+z-2*y)*p
? ?return u
for k in range(0,M):
? ?ans[k][0]=100
? ?ans[k][N-1]=100
print(ans)
for j in range(1,M):
? ?for i in range(1,N-1):
? ? ? ?ans[j][i]=ans[j-1][i]+jisuan(ans[j-1][i+1],ans[j-1][i],ans[j-1][i-1])

print(ans)
# plt.matshow(ans, cmap=plt.cm.Reds)#這里設(shè)置顏色為紅色，也可以設(shè)置其他顏色
# plt.title("matrix A")
# plt.show()
fig=plt.figure()
ax=ax(fig)
Y=np.arange(0,t,deltat)
X=np.arange(0,lx,deltax)
X,Y=np.meshgrid(X,Y)

plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# mpl.rcParams['legend.fontsize'] = 30 ?#線段標(biāo)簽 字體大小

ax.set_xlabel('空間長度',size=10,labelpad=10)
ax.set_ylabel('時(shí)間/s',size=10,labelpad=10)
ax.set_zlabel('溫度/C',size=10,labelpad=10)
ax.plot_surface(X, Y, ans, ?cmap=plt.get_cmap('rainbow'))
plt.show()
案例2：

來源于：微分方程數(shù)值求解——有限差分法 - 知乎 (zhihu.com)

具體問題如下：

邊界條件如下：

差分如下：

默認(rèn)初始場為0；

python中，x步長為0.01，y步長為0.01時(shí)：

x步長為0.025，y步長為0.025時(shí)：

很明顯最大值不同了！

代碼如下：

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D as ax
import matplotlib.animation as ma
pi=3.1415926
deltax=0.025
deltay=0.025
x=1
y=1
p=1/(deltax*deltax)
q=1/(deltay*deltay)
I=int(x/deltax)
J=int(y/deltay)
ans=np.zeros((J,I))
print(ans)
def f(x,y):
? ?u=(-2)*pi*pi*math.sin(pi*x)*math.sin(pi*y)
? ?return u
for j in range(1,J-1):

? ?for i in range(1,I-1):
? ? ? ?# print(i,j)
? ? ? ?ans[j][i]=(p*(ans[j][i+1]+ans[j][i-1])+q*(ans[j+1][i]+ans[j-1][i])-f(i*deltax,j*deltay))/(2*(p+q))
? ? ? ?if ans[j][i]<0:
? ? ? ? ? ?ans[j][i]=0

fig = plt.figure()
ax = ax(fig)
Y = np.arange(0, y, deltay)
X = np.arange(0, x, deltay)
X, Y = np.meshgrid(X, Y)
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# mpl.rcParams['legend.fontsize'] = 30 ?#線段標(biāo)簽 字體大小

ax.set_xlabel('x',size=10,labelpad=10)
ax.set_ylabel('y',size=10,labelpad=10)
ax.set_zlabel('',size=10,labelpad=10)
ax.plot_surface(X, Y, ans, cmap=plt.get_cmap('rainbow'))
plt.show()
# plt.pause(100)
# plt.close("all")
# print(ans)


# plt.matshow(ans, cmap=plt.cm.get_cmap("rainbow"))#這里設(shè)置顏色為紅色，也可以設(shè)置其他顏色
# plt.title("matrix A")
# plt.show()

希望大佬解惑！