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

水字数用,哥以后转战其他网站了.这个评分系统太恶心了.太恶心了.
其实有个很有趣的现象.没考过计算机二级的人觉得二级Python弱爆了.一点用没有,备考计算机二级Python的人觉得难暴了.......
挺好的,说不出话来......
其实我个人认为非专业人士考一次计算机二级有助于学习的,还可以水创新学分.真是想不通哪个SX想出的创新学分.不是,这个人啊,他创新,创新这个是用指标逼出来的吗?

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