二维Poisson方程五点差分格式与Python实现_用有限差分方法(五点差分格式)求解正方形域上的poisson方程边值问题

• 最近没怎么写新文章，主要在学抽象代数
• 下学期还有凸分析
• 好累的一学期
• 哦对，我不是数学系的，我是物理系的。而且博主需要澄清一下，博主没有对象，至少现在还没有。

---

---

• 好，兄弟们，好习惯，先上代码后说话！

Python 实现

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.sparse as ss
from scipy.sparse.linalg import spsolve
 
class PDE2DModel:
    #均应该传入一个一维二元数组，表示起止值
    def __init__(self,x,y):
        assert len(x)==len(y)==2,"ERROR:UNEXPECTED SV and EV!!"
        self.x = x
        self.y = y
 
    #hx表示X上的步长
    #hy表示Y上的步长
    def space_grid(self,hx,hy):
        M = int(round((self.x[1]-self.x[0])/hx,0))
        N = int(round((self.y[1]-self.y[0])/hy,0))
        assert M==N>=3,"至少网格数是合理的"
        X = np.linspace(self.x[0],self.x[1],M+1)
        Y = np.linspace(self.y[0],self.y[1],N+1)
        return M,N,X,Y
 
    def f(self,X,Y):
        return 6*X*Y**3+6*X**3*Y+np.e**X*np.sin(Y)-np.e**X*np.sin(Y)
 
    def solution(self,X,Y):
        return np.e**X*np.sin(Y)-X**3*Y**3
 
    #左边界
    def left(self,Y):
        return np.sin(Y)
    #右边界
    def right(self,Y):
        return np.e**3*np.sin(Y)-27*Y**3
    #上边界
    def up(self,X):
        return np.sin(1)*np.e**X-X**3
    #下边界
    def down(self,X):
        return 0*X
 
#解算核心
def NDM5_2D(PDE2DModel,hx,hy):
    M,N,X0,Y0 = PDE2DModel.space_grid(hx,hy)
    Y,X = np.meshgrid(Y0,X0)
##    print("X0",X0)
##    print("Y0",Y0)
##    #数值结果保存在U中 从0到N共N+1个
##    print("M",M)
##    print("N",N)
    U = np.zeros((M+1,N+1))
    U[0,:]  = PDE2DModel.left(Y0)
    U[-1,:] = PDE2DModel.right(Y0)
    U[:,0]  = PDE2DModel.down(X0)
    U[:,-1] = PDE2DModel.up(X0)
 
    D = np.diag([-1/(hy**2) for i in range(M-1)])
    C = np.zeros((N-1,N-1),dtype="float64")
    for i in range(N-1):
        C[i][i] = 2*(1/hx**2+1/hy**2)
        if i<N-2:
            C[i][i+1] = -1/hx**2
            C[i+1][i] = -1/hx**2
 
 
    u0 = np.array([[PDE2DModel.down(X0[i])] for i in range(1,M)])
    un = np.array([[PDE2DModel.up(X0[i])] for i in range(1,M)])
    
    F = np.zeros((M-1)*(N-1)) 
    for j in range(1,N):
        #for i in range(1,M):
        F[(N-1)*(j-1):(N-1)*(j)] = PDE2DModel.f(X0[1:M],np.array([Y0[j] for i in range(N-1)]))
 
        F[(N-1)*(j-1)] += PDE2DModel.left(Y0[j])/hx**2
        F[(N-1)*(j)-1] += PDE2DModel.right(Y0[j])/hx**2
 
    F[:N-1] -= np.dot(D,u0).T[0]
    F[(M-1)*(N-2):] -= np.dot(D,un).T[0]
    F = np.mat(F).T
##    print(F)
 
    Dnew = np.zeros(((M-1)*(N-1),(N-1)*(M-1)))
    for i in range((N-1)*(N-1)):
        Dnew[i][i] = 2*(1/hx**2+1/hy**2)
 
        if i<(N-2)*(N-1):
            Dnew[i][i+N-1] = -1/hy**2
            Dnew[i+N-1][i] = -1/hy**2
 
    for i in range(N-1):
        for j in range(N-2):
            Dnew[(N-1)*i+j][(N-1)*i+j+1] = -1/hx**2
            Dnew[(N-1)*i+j+1][(N-1)*i+j] = -1/hx**2
 
    print("差分方程构造完成！解算开始！")  
    Unew = np.linalg.solve(Dnew,F)
    #print(Unew)
    U[1:-1,1:-1] = Unew[:,0].reshape((N-1,N-1)).T
    return X,Y,U
 
#数据可视化
def Visualized():
    x = np.array([0,3])
    y = np.array([0,1])
 
    pde = PDE2DModel(x,y)
    X,Y,U = NDM5_2D(pde,0.03,0.01)
    u = pde.solution(X,Y)
 
    print("解算完成！绘图已开始！")
    plt.figure(figsize=(15,5))
    ax1 = plt.subplot(131,projection="3d")
    ax2 = plt.subplot(132,projection="3d")
    ax3 = plt.subplot(133,projection="3d")
 
    ax1.set_title("Numeric Solution")
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.plot_surface(X,Y,U,cmap="gist_ncar")
 
    ax2.set_title("Exact Solution")
    ax2.set_xlabel("x")
    ax2.set_ylabel("y")
    ax2.plot_surface(X,Y,u,cmap="gist_ncar")
 
    e = np.abs(U-u)
    ax3.set_title("Error")
    ax3.set_xlabel("x")
    ax3.set_ylabel("y")
    ax3.plot_surface(X,Y,e,cmap="gist_ncar")
 
    plt.show()
    return U,u,X,Y
 
U,u,X,Y = Visualized()

---

• 好习惯2,效果图

很好,接下去接着讲

五点差分格式

•  这里我要用截图,截的我的小论文.没办法,CSDN写作太麻烦了,没有Latex好用.我记得这个好像可以调整,但我没学会

 

 很好,结束了.

注意事项

• 你在写代码的时候,容易犯几个错误.
• 在纸上图画出来了,然后呢,纸上是右手系,写进数组了,就忘了是啥系了!
• 比如:

    Y,X = np.meshgrid(Y0,X0)

• 我们需要重点看一下np.meshgrid问题.这才是大家想要的数据结构.....

import numpy as np
x = np.array([1,2,3,4,5,6,7,8,9,10])
y = np.array([11,12,13,14,15,16,17,18,19,20])
 
Y,X = np.meshgrid(y,x)

•  其实这个还涉及到稀疏矩阵的问题,稀疏矩阵解决不了的话,你的运算量会变得超大!一定要用系数矩阵优化才能解决问题!!

所以,优化之后!

用稀疏矩阵优化

import numpy as np
import scipy as sp
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix,csc_matrix
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.sparse as ss
from scipy.sparse.linalg import spsolve
import  time
 
class PDE2DModel:
    def __init__(self):
        self.x = np.array([0,3])
        self.y = np.array([0,1])
 
    def space_grid(self,hx,hy):
        M = int(round((self.x[1]-self.x[0])/hx,0))
        N = int(round((self.y[1]-self.y[0])/hy,0))
        assert M==N>=3,"ERROR:UNECPECTED GRIDS M:"+str(M)+" N:"+str(N)
        X = np.linspace(self.x[0],self.x[1],M+1)
        Y = np.linspace(self.y[0],self.y[1],N+1)
        return M,N,X,Y
 
    def f(self,X,Y):
        return 6*X*Y**3+6*X**3*Y+np.e**X*np.sin(Y)-np.e**X*np.sin(Y)
 
    #左边界
    def left(self,Y):
        return np.sin(Y)
    #右边界
    def right(self,Y):
        return np.e**3*np.sin(Y)-27*Y**3
    #上边界
    def up(self,X):
        return np.sin(1)*np.e**X-X**3
    #下边界
    def down(self,X):
        return 0*X
 
 
def NDM5_2D(PDE2DModel,hx,hy):
 
    M,N,X0,Y0 = PDE2DModel.space_grid(hx,hy)
    Y,X = np.meshgrid(Y0,X0)
 
    U = np.zeros((M+1,N+1))
    U[0,:]  = PDE2DModel.left(Y0)
    U[-1,:] = PDE2DModel.right(Y0)
    U[:,0]  = PDE2DModel.down(X0)
    U[:,-1] = PDE2DModel.up(X0)
 
    D = np.diag([-1/(hy**2) for i in range(M-1)])
    C = np.zeros((N-1,N-1),dtype="float64")
    for i in range(N-1):
        C[i][i] = 2*(1/hx**2+1/hy**2)
        if i<N-2:
            C[i][i+1] = -1/hx**2
            C[i+1][i] = -1/hx**2
 
    u0 = np.array([[PDE2DModel.down(X0[i])] for i in range(1,M)])
    un = np.array([[PDE2DModel.up(X0[i])] for i in range(1,M)])
    
    F = np.zeros((M-1)*(N-1)) 
    for j in range(1,N):
        F[(N-1)*(j-1):(N-1)*(j)] = PDE2DModel.f(X0[1:M],np.array([Y0[j] for i in range(N-1)]))
 
        F[(N-1)*(j-1)] += PDE2DModel.left(Y0[j])/hx**2
        F[(N-1)*(j)-1] += PDE2DModel.right(Y0[j])/hx**2
 
    F[:N-1] -= np.dot(D,u0).T[0]
    F[(M-1)*(N-2):] -= np.dot(D,un).T[0]
    F = np.mat(F).T
 
    Dnew = np.zeros(((M-1)*(N-1),(N-1)*(M-1)))
    for i in range((N-1)*(N-1)):
        Dnew[i][i] = 2*(1/hx**2+1/hy**2)
 
        if i<(N-2)*(N-1):
            Dnew[i][i+N-1] = -1/hy**2
            Dnew[i+N-1][i] = -1/hy**2
 
    for i in range(N-1):
        for j in range(N-2):
            Dnew[(N-1)*i+j][(N-1)*i+j+1] = -1/hx**2
            Dnew[(N-1)*i+j+1][(N-1)*i+j] = -1/hx**2
 
    print("差分方程构造完成！解算开始！")
##    start = time.time()
##    Unew = np.linalg.solve(Dnew,F)
##    U[1:-1,1:-1] = Unew[:,0].reshape((N-1,N-1)).T
##    end = time.time()
##    t = end-start
 
    start = time.time()
    SDnew = csr_matrix(Dnew)
    SF = csr_matrix(F)
    SUnew = spsolve(SDnew,SF)
    U[1:-1,1:-1] = SUnew.reshape((N-1,N-1)).T
    end = time.time()
    t = end-start
    return X,Y,U,t
 
def Visualized(X,Y,U):
 
    print("解算完成！绘图已开始！")
    plt.figure(figsize=(15,5))
    ax1 = plt.subplot(111,projection="3d")
    ax1.set_title("Numeric Solution")
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.plot_surface(X,Y,U,cmap="gist_ncar")
    plt.show()
 
pde = PDE2DModel()
X,Y,U,t = NDM5_2D(pde,0.03,0.01)
print(t)
Visualized(X,Y,U)

Python技巧

map

f = lambda x:x+5
X = [1,2,3,4,5]
map(f,X)
<map object at 0x00000246B5CE6110>
list(map(f,X))
[6, 7, 8, 9, 10]

Python自制警告

• Python 自制警告是很有必要的.
• 自制的警告与自制的错误不一样,不会打断程序运行(显然废话)

import warnings
a = 0
if a==0:
    warnings.warn("IGNORE!")
print(a+3)

Warning (from warnings module):
  File "C:\Users\LX\Desktop\新建 文本文档.py", line 34
    warnings.warn("IGNORE!")
UserWarning: IGNORE THIS WARNING!
3

Python matplotlib技巧

##import matplotlib
##import matplotlib.pyplot as plt
##from mpl_toolkits.mplot3d import Axes3D
##matplotlib.rcParams["font.sans-serif"] = ["SimHei"]
##matplotlib.rcParams["axes.unicode_minus"] = False
##
##import numpy as np
##
##x = np.linspace(1,2,6)
##y = np.linspace(3,9,13)
##Y,X = np.meshgrid(y,x)
##X2,Y2 = np.meshgrid(x,y)
##
##plt.figure(figsize=(15,5))
##ax1 = plt.subplot(121,projection="3d")
##ax2 = plt.subplot(122,projection="3d")
##
##ax1.set_title("Inverse order")
##ax2.set_title("Ordinary order")
##ax1.set_xlabel("x")
##ax2.set_xlabel("x")
##
##f = lambda x,y:x*2+y*3
##ax1.plot_surface(X,Y,f(X,Y))
##
##ax2.plot_surface(X2,Y2,f(X2,Y2))
##
##plt.show()