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 spsolveclass PDE2DModel:#均应该传入一个一维二元数组，表示起止值def __init__(self,x,y):assert len(x)==len(y)==2,"ERROR:UNEXPECTED SV and EV!!"self.x = xself.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,Ydef 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**2C[i+1][i] = -1/hx**2u0 = 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**2F[(N-1)*(j)-1] += PDE2DModel.right(Y0[j])/hx**2F[: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**2Dnew[i+N-1][i] = -1/hy**2for i in range(N-1):for j in range(N-2):Dnew[(N-1)*i+j][(N-1)*i+j+1] = -1/hx**2Dnew[(N-1)*i+j+1][(N-1)*i+j] = -1/hx**2print("差分方程构造完成！解算开始！")  Unew = np.linalg.solve(Dnew,F)#print(Unew)U[1:-1,1:-1] = Unew[:,0].reshape((N-1,N-1)).Treturn 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,YU,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  timeclass 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,Ydef 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*Xdef 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**2C[i+1][i] = -1/hx**2u0 = 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**2F[(N-1)*(j)-1] += PDE2DModel.right(Y0[j])/hx**2F[:N-1] -= np.dot(D,u0).T[0]F[(M-1)*(N-2):] -= np.dot(D,un).T[0]F = np.mat(F).TDnew = 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**2Dnew[i+N-1][i] = -1/hy**2for i in range(N-1):for j in range(N-2):Dnew[(N-1)*i+j][(N-1)*i+j+1] = -1/hx**2Dnew[(N-1)*i+j+1][(N-1)*i+j] = -1/hx**2print("差分方程构造完成！解算开始！")
##    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-startstart = 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)).Tend = time.time()t = end-startreturn X,Y,U,tdef 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 34warnings.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()

二维Poisson方程五点差分格式与Python实现相关推荐

二维Poisson方程五点差分格式及简单求解方法Python实现
二维Poisson方程简介给出二维 PoissonPoissonPoisson 方程 DirichletDirichletDirichlet 边值问题: {−Δu=f(x,y)(x,y)∈Ωu=φ ...
二维有限元方程matlab,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维 Poisson 方程的 MATLAB 实现陈莲a ,郭元辉b ,邹叶童a ( 西华师范大学 a．数学与信息学院; b．教育信息技术中心,四川南充 6437009) 摘要: ...
galerkin有限元法matlab实现,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维Poisson方程的MATLAB实现陈莲a,郭元辉b,邹叶童a [摘要]文章讨论了圆形区域上的三角形单元剖分.有限元空间,通过变分形式离散得到有限元方程. 用MATLAB编程求得数值 ...
二维泊松方程数值解-五点差分法-共轭梯度法-python实现
松方程有很多现成的工具可以用,这里主要是为了加深对算法的理解.题目如下题目的要点在于找到泊松方程的系数矩阵.在五点法里面,系数矩阵一共五条对角线,一条主对角线,四条副对角线.碰到边界的时候有的对角线 ...
Poisson方程五点差分格式例题及解答
写在前面例题解题过程总结参考写在前面本文针对偏微分方程数值解中出现的一道例题进行分析,详细介绍了五点差分格式的公式推导及应用. 例题在单位正方形Ω‾:0⩽x⩽1,0⩽y⩽1\overli ...
二维声波方程的有限差分法数值模拟
二维声波方程的有限差分法数值模拟文章目录二维声波方程的有限差分法数值模拟一.实现效果二.matlab代码分享三.python代码分享一.实现效果二.matlab代码分享 close al ...
二维burgers方程_二维Burgers方程的RKDG有限元解法
二维 Burgers 方程的 RKDG 有限元解法 ∗ 马艳春 1, 张寅虎 2, 冯新龙 1 [摘要] 摘要 : 本文应用 RKDG 有限元方法求解具有周期边界条件的二维非粘性 Burgers ...
二维burgers方程_用格子Boltzmann方法研究二维Burgers方程
用格子 Boltzmann 方法研究二维 Burgers 方程张伟 ; 李文杰 [期刊名称] <天津城市建设学院学报> [年 ( 卷 ), 期] 2012(018)001 [ 摘要 ] ...
有限元方法基础-以二维拉普拉斯方程为例(附程序)
文章目录前言问题描述问题区域偏微分方程变分问题(弱形式) 有限元离散二维线性有限元 : 参考基函数 2D linear finite element : affine mapping 有限 ...

二维Poisson方程五点差分格式与Python实现