用五点有限差分格式求解椭圆型方程(偏微分方程) 程序2

有限差分法求解偏微分方程摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。

关键词：计算力学，偏微分方程，有限差分法Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method.Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Method1 引言机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。

通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。

求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。

另一方面，计算机技术的发展使得计算更精确、更迅速。

五点差分格式求解poisson方程在数学和科学领域中，Poisson方程是一种常见的偏微分方程，通常用于描述在某个区域内的物理现象。

求解Poisson方程在科学计算和数值模拟中具有重要的应用。

本文将介绍一种常用的数值求解方法——五点差分格式，用于求解Poisson方程。

1. 引言Poisson方程是二阶偏微分方程，形式如下：∇²u = f(x, y)其中，u是未知函数，f(x, y)是已知函数，∇²是拉普拉斯算子。

求解Poisson方程的目标是找到满足方程的函数u。

2. 五点差分格式五点差分格式是一种常用的离散化方法，用于将连续的Poisson方程转化为离散的数值方程。

在二维情况下，我们可以将区域划分为若干个格点，然后利用差分近似来求解。

假设我们在二维区域(0, 1) x (0, 1)内选择了一组均匀的格点(xi, yj)，其中i和j分别表示x和y方向上的索引。

则方程中的拉普拉斯算子可以被近似为：(∂²u/∂x²) ≈ (u(i+1,j) - 2u(i,j) + u(i-1,j))/Δx²(∂²u/∂y²) ≈ (u(i,j+1) - 2u(i,j) + u(i,j-1))/Δy²其中Δx和Δy是格点的空间间隔。

将上述近似代入Poisson方程，我们可以得到离散的数值方程：(u(i+1,j) - 2u(i,j) + u(i-1,j))/Δx² + (u(i,j+1) - 2u(i,j) + u(i,j-1))/Δy² = f(i,j)3. 离散方程求解根据五点差分格式，我们可以得到离散的数值方程，进而求解Poisson方程。

为了方便计算，我们可以定义一个N x N的矩阵A，其中N表示在x和y方向上的格点数目。

矩阵A中的每个元素对应于方程中的一个未知数。

我们可以将方程表示为矩阵形式：Au = b其中，u是未知函数的向量，b是已知函数f的向量。

有限差分法解偏微分方程综述绪论有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。

在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。

有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。

该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。

有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。

该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。

对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。

从差分的空间形式来考虑，可分为中心格式和逆风格式。

考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。

目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。

差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。

构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。

其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。

通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。

有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。

《微分方程数值解》大作业（一）——椭圆型方程编程计算：采用五点差分格式求如下椭圆型方程2222uu x y f (x,y),(x,y);∂∂∂∂--=∈Ω其中f (x,y)、Ω及边条件为：1． f (x,y)0,= (1,2)(0,1)Ω=⨯， 且边条件如下：222u(x,0)2ln x,u(x,1)ln(x 1)1x 2;u(1,y)ln(1y ),u(2,y)ln(4y ),0y 1.⎧==+<<⎪⎨=+=+<<⎪⎩， 问题存在精确解为： 22(,)ln()u x y x y =+2．f (x,y)4,=- (0,1)(0,2)Ω=⨯，且边条件如下：2222u(x,0)x ,u(x,2)(x 2)0x 1;u(0,y)y ,u(1,y)(y 1),0y 2.⎧==-<<⎪⎨==-<<⎪⎩， 问题存在精确解为： 2(,)()u x y x y =-3．f (x,y)cos(x y)cos(x y),=++- (0,)(0,)2πΩ=π⨯，且边条件如下：u(x,0)cos x,u(x,)00x ;2u(0,y)cos y,u(,y)cos y,0y .2π⎧==<<π⎪⎪⎨π⎪=π=-<<⎪⎩， 问题存在精确解为： (,)cos cos u x y x y =.代码：主函数1，差分解function g=fivepoints(x1,x2,y1,y2,M,N)%变步长法h=(x2-x1)/M; %横轴步长k=(y2-y21/N; %纵轴步长m=M-1;n=N-1;h1=h^2;r=h1/k^2; %五点中的上下两个点的系数t=2+2*r; %五点中的中心点的系数x=x1+(x2-x1)*(0:M)/M; %x，y向量表示横纵坐标y=y1+(y2-y1)*(0:N)/N;a=zeros(m*n,m*n);b=zeros(m*n,1);%初始化a，b矩阵，a为系数矩阵%内部的（m-2）*（n-2）个点for i=2:m-1for j=2:n-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m) -r zeros(1,m-2) -1 t -1 zeros(1,m-2) -r zeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1));endend%下边缘j=1;for i=2:m-1a(i+(j-1)*m,:)=[zeros(1,i-2) -1 t -1 zeros(1,m-2) -r zeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*bottom(x(i+1));end;%右边缘i=m;for j=2:n-1a(i+(j-1)*m,:)=[zeros(1,(j-1)*m-1) -r zeros(1,m-2) -1 t zeros(1,m-1) -r zeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+right(y(j+1));end%上边缘j=n;for i=2:m-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m) -r zeros(1,m-2) -1 t -1 zeros(1,m-i-1)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1));end%左边缘i=1;for j=2:n-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m) -r zeros(1,m-1) t -1 zeros(1,m-2) -rzeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+left(y(j+1));end;%左下角的那个点i=1;j=1;a(1,:)=[t -1 zeros(1,m-2) -r zeros(1,(n-1)*m-1)];b(1)=h1*f(x(2),y(2))+r*bottom(x(2))+left(y(2));%右下角的那个点i=m;j=1;a(i+(j-1)*m,:)=[zeros(1,m-2) -1 t zeros(1,m-1) -r zeros(1,(n-2)*m)]; b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*bottom(x(i+1))+right(y(j+1)); %左上角的那个点i=1;j=n;a(i+(j-1)*m,:)=[zeros(1,(n-2)*m) -r zeros(1,m-1) t -1 zeros(1,m-2)]; b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1))+left(y(j+1));%右上角的那个点i=m;j=n;a(i+(j-1)*m,:)=[zeros(1,(n-1)*m-1) -r zeros(1,m-2) -1 t];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1))+right(y(j+1));u=a\bab2，精确解：function g=ni(x1,x2,y1,y2,M,N)m=M-1;n=N-1;x=x1+(x2-x1)*(0:M)/M;y=y1+(y2-y1)*(0:N)/N;for i=1:mfor j=1:nu1(i+(j-1)*m)=f1(x(i+1),y(j+1))endend（1）辅助函数function g=f(x,y)g=0;function g=bottom(x)g=2*log(x);function g=right(y)g=log(4+y^2);function g=top(x)g=log(x^2+1);function g=left(y)g=log(1+y^2);function g=f1(x,y)g=log(x^2+y^2);运行fivepoints(1,2,0,1,4,4)u =数值解0.4847467147016780.8376456266975491.1390195099193150.5944295076643080.9158860659528741.1974022894530100.7539416986884711.0340668399966291.287784599003526a =4 -1 0 -1 0 0 0 0 0 -1 4 -1 0 -1 0 0 0 0 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 0 0 0 -1 0 -1 4 -1 0 -1 0 0 0 -1 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 0 0 0 0 -1 0 -1 4 -1 0 0 0 0 0 -1 0 -1 4b =0.5069117244448540.8109302162163292.5210301235267010.2231435513142101.4469189829363251.3872704470929461.1786549963416462.919669266564466运行ni(1,2,0,1,4,4)u1 =精确解Columns 1 through 30.485507815781701 0.838329190404443 1.139434283188365 Columns 4 through 60.594707107746693 0.916290731874155 1.197703191312341 Columns 7 through 90.753771802376380 1.034073767530539 1.287854288306638 误差很小（2）辅助函数function g=f(x,y)g=-4;function g=bottom(x)g=x^2;function g=right(y)g=(y-1)^2;function g=top(x)g=(x-2)^2;function g=left(y)g=y^2;function g=f1(x,y)g=(x-y)^2;fivepoints(1,2,0,1,4,4)fivepoints(0,1,0,2,4,4)u =0.062500000000000-0.0000000000000000.0625000000000000.5625000000000000.2500000000000000.0625000000000001.5625000000000001.0000000000000000.562500000000000a =Columns 1 through 32.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000-1.0000000000000000 -1.000000000000000 2.500000000000000 -0.250000000000000 0 00 -0.250000000000000 00 0 -0.2500000000000000 0 00 0 00 0 0Columns 4 through 6-0.250000000000000 0 00 -0.250000000000000 00 0 -0.2500000000000002.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000-1.0000000000000000 -1.000000000000000 2.500000000000000 -0.250000000000000 0 00 -0.250000000000000 00 0 -0.250000000000000Columns 7 through 90 0 00 0 00 0 0-0.250000000000000 0 00 -0.250000000000000 00 0 -0.2500000000000002.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000-1.0000000000000000 -1.000000000000000 2.500000000000000b =0.015625000000000-0.1875000000000000.1406250000000000.750000000000000-0.250000000000000-0.2500000000000002.7656250000000000.3125000000000000.390625000000000精确解ni(0,1,0,2,4,4)u1 =u1 =Columns 1 through 30.062500000000000 0 0.062500000000000 Columns 4 through 60.562500000000000 0.2500000000000000.062500000000000Columns 7 through 91.562500000000000 1.0000000000000000.562500000000000误差很小（3）辅助函数function g=f(x,y)g=cosd(x+y)+cosd(x-y);function g=bottom(x)g=cosd(x);function g=right(y)g=-cosd(y);function g=top(x)g=0;function g=left(y)g=cosd(y);function g=f1(x,y)g=cosd(x)*cosd(y);数值解Pi=3.1415926fivepoints(0,pi,0,pi/2,4,4)u =0.6578183624886530.000000024999241-0.6578183271343870.5049807980892560.000000019229497-0.5049807708946410.2736443626241530.000000010432161-0.273644347870850a =10 -1 0 -4 0 0 0 0 0 -1 10 -1 0 -4 0 0 0 0 0 -1 10 0 0 -4 0 0 0 -4 0 0 10 -1 0 -4 0 0 0 -4 0 -1 10 -1 0 -4 0 0 0 -4 0 -1 10 0 0 -4 0 0 0 -4 0 0 10 -1 0 0 0 0 0 -4 0 -1 10 -1 0 0 0 0 0 -4 0 -1 10b =4.5582604075302670.000000137720159-4.5582602127645491.323957*********0.000000023374742-1.3239570281549570.7165204234523470.000000012650320-0.716520405562093精确解ni(0,pi,0,pi/2,4,4)u1 =Columns 1 through 30.653281493003155 0.000000024755257-0.653281457993935Columns 4 through 60.500000013397448 0.000000018946853-0.499999986602551Columns 7 through 90.270598066826879 0.000000010253963-0.270598052325585误差很小注：（1）需要对数值解与精确解作比较，以及不同步长选取下的误差比较。

数学与计算科学学院实验报告
实验项目名称椭圆型方程数值解
所属课程名称微分方程数值解法
实验类型验证
实验日期
班级信计0902
学号
姓名
成绩
附录1：源程序
附录2：实验报告填写说明
1．实验项目名称：要求与实验教学大纲一致。

2．实验目的：目的要明确，要抓住重点，符合实验教学大纲要求。

3．实验原理：简要说明本实验项目所涉及的理论知识。

4．实验环境：实验用的软、硬件环境。

5．实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容。

概括整个实验过程。

对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。

对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设
计思路和设计方法，再配以相应的文字说明。

对于创新性实验，还应注明其创新点、特色。

6．实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。

7．实验结论（结果）：根据实验过程中得到的结果，做出结论。

8．实验小结：本次实验心得体会、思考和建议。

9．指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。

用差分法解椭圆型偏微分方程-（Uxx+Uyy)=(pi*pi-1)e^xsi n(pi*y) 0<x<2; 0<y<1U(0,y)=sin(pi*y),U(2,y)=e^ 2sin(pi*y); 0=<y<=1 U(x,0)=0, U(x,1)=0; 0=<x<=2先自己去看一下关于五点差分法的理论书籍Matlab程序：unction[p e u x y k]=wudianchafenfa(h,m, n,kmax,ep)% g-s迭代法解五点差分法问题%kmax为最大迭代次数%m,n为x,y方向的网格数，例如（2-0）/0.01=200;%e为误差，p为精确解syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h;y=0+(0:n)*h;for(i=1:n+1)u(i,1)=sin(pi*y(i));u(i,m+1)=exp(1)*exp(1) *sin(pi*y(i));endfor(i=1:n)for(j=1:m)f(i,j)=(pi*pi-1)*exp(x (j))*sin(pi*y(i));endendt=zeros(n-1,m-1);for(k=1:kmax)for(i=2:n)for(j=2:m)temp=h*h*f(i,j)/4+(u(i ,j+1)+u(i,j-1)+u(i+1,j )+u(i-1,j))/4;t(i,j)=(temp-u(i,j))*( temp-u(i,j));u(i,j)=temp;endendt(i,j)=sqrt(t(i,j));if(k>kmax)break;endif(max(max(t))<ep)break;endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j))*sin(p i*y(i));e(i,j)=abs(u(i,j)-exp( x(j))*sin(pi*y(i)));endEnd在命令窗口中输入：[p e u x y k]=wudianchafenfa(0.1,20, 10,10000,1e-6) k=147 surf(x,y,u) ;xlabel(‘x’);ylabel(‘y’); zlabel(‘u’);Title(‘五点差分法解椭圆型偏微分方程例1’)就可以得到下图ｓｕｒｆ（ｘ，ｙ，ｐ）ｓｕｒｆ（ｘ，ｙ，ｅ）[p e u x y k]=wudianchafenfa(0.05,4 0,20,10000,1e-6)[p e u x y k]=wudianchafenfa(0.025,80,40,10000,1e-6)为什么分得越小，误差会变大呢？我们试试运行：[p e u x y k]=wudianchafenfa(0.025, 80，40,10000,1e-8)K=2164surf(x,y,e)误差变小了吧还可以试试[p e u x y k]=wudianchafenfa(0.025, 80，40,10000,1e-10)K=3355误差又大了一点再试试[p e u x y k]=wudianchafenfa(0.025, 80,40,10000,1e-11)k=3952误差趋于稳定总结：最终的误差曲面与网格数有关，也与设定的迭代前后两次差值（ep,看程序）有关；固定网格数，随着设定的迭代前后两次差值变小，误差由大比变小，中间有一个最小值，随着又增大一点，最后趋于稳定。

有限差分法解偏微分方程综述绪论有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。

在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。

有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。

该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。

有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。

该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。

对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。

从差分的空间形式来考虑，可分为中心格式和逆风格式。

考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。

目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。

差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。

构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。

其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。

通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。

有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。

matlab有限差分法求解椭圆型偏微分方程
有限差分法是一种求解偏微分方程的经典数值方法，它将连续的
偏微分方程转化为离散的代数方程，从而能够使用计算机进行计算。

在 MATLAB 中，我们可以使用有限差分法来求解椭圆型偏微分方程。

椭圆型偏微分方程通常用来描述有稳态的空间分布的物理现象，
如稳态的温度分布。

其通用的数学形式为：
∇·(a(x,y)∇u(x,y)) + f(x,y) = 0
其中，u(x,y) 是要求解的函数，a(x,y) 是定义在区域Ω上的
函数，它代表了该区域内各点的材料特性，f(x,y) 是特定的源项函数。

有限差分法将区域Ω划分为离散的点集，然后通过对这些点之
间的差分运算进行逐点计算，得到离散式。

例如，可以使用中心差分
法对 u(x,y) 在某个点(x0,y0) 的二阶偏导数进行离散化，得到：(u(x0+Δx,y0) - 2u(x0,y0) + u(x0-Δx,y0)) / Δx^2
同样，对于 a(x,y)在点(x0,y0)的取值，我们也可以使用中心差
分法进行离散化：
(a(x0+Δx,y0) + a(x0,y0)) / 2
经过离散化后，我们可以将偏微分方程变为一个线性代数方程组，使用 MATLAB 的矩阵运算功能进行求解。

需要注意的是，在实际计算中，由于矩阵求逆时存在数值不稳定的问题，因此需要对矩阵进行一
定的处理，如使用迭代法或预处理技术等。

总之，有限差分法是一种常用的求解偏微分方程的数值方法，在MATLAB 中也有相应的实现。

通过离散化连续的偏微分方程，我们能够
在计算机上高效地求解椭圆型偏微分方程，提高计算效率，解决实际
问题。

matlab编程解Laplace偏微分方程(五点差分法)应用五点差分法解Laplace方程,帮你一劳永逸的解决问题!functionu=Laplace(D1,b某0,b某f,by0,byf,n某,ny,M,delta) min某=D1(1);ma某某=D1(2);miny=D1(3);ma某y=D1(4);h某=(ma某某-min某)/(n某-1);hy=(ma某y-miny)/(ny-1);某=min某:h某:ma某某;y=miny:hy:ma某y;forj=1:n某u0(1,j)=by0(某(j));u0(ny,j)=byf(某(j));endfori=2:(ny-1)u0(i,1)=b某0(y(i));u0(i,n某)=b某f(y(i));endfori=1:n某um=u0(1,n某)+u0(ny,n某);endforj=2:ny-1um=u0(j,1)+u0(j,n某);endfori=2:n某-1forj=2:ny-1u0(j,i)=um/(2某(n某+ny-2));endendu1=u0;k=0;form=1:Mforj=2:(n某-1)fori=2:(ny-1)u0(i,j)=1.7某(u0(i-1,j)+u0(i,j-1)+u1(i+1,j)+u1(i,j+1)-4某u1(i,j))/4+u1(i,j);endendifma某(ma某(ab(u0-u1)))<deltabreak;eleu1=u0;endk=k+1;endu0;printf('迭代次数：k=%d',k) ubplot(2,2,1)meh(某,y,u0)ubplot(2,2,2)urf(某,y,u0)hadinginterp某label('某')ylabel('y')zlabel('u0')ubplot(2,2,3)[c,h]=contour3(某,y,u0); clabel(c,h)ubplot(2,2,4)[m,n]=contour(某,y,u0); clabel(m,n)while(1)应用五点差分法解Laplace方程,帮你一劳永逸的解决问题!dip('是否求某点函数值是请输入1否请输入0')a=input('a='); ifa==1dip('输入坐标:')某=input('某=');y=input('y=');if某<min某|某>ma某某|y<miny|y>ma某ydip('超出坐标范围') elep=round(1+某/h某);q=round(1+y/hy);printf('u(某,y)=%f',u0(q,p))endeledip('bye-bye')break;endend程序中的边界条件详见程序代码e某1401Laplace方程边界条件。

椭圆型方程的差分解法1.引言考虑问题①二维Poisson 方程2222(,)u u f x y x y ⎛⎫∂∂-+= ⎪∂∂⎝⎭， (,)x y ∈Ω 其中Ω为2R 中的一个有界区域，其边界Γ为分段光滑曲线。

在Γ上u 满足下列边界条件之一：⑴(,)u x y αΓ=（第一边值条件）， ⑵(,)ux y n βΓ∂=∂（第二边值条件）， ⑶(,)uku x y n γΓ∂+=∂（第三边值条件）， (,),(,),(,),(,),(,)f x y x y x y x y k x y αβγ都是连续函数，0k ≥.2.差分格式将区间[,]a b 作m 等分，记为11()/,,0i h b a m x a ih i m =-=+≤≤；将区间[,]c d 作n 等分，记为22()/,,0i h d c n y c jh j n =-=+≤≤.称1h 为x 方向的步长，2h 为y 方向的步长。

2.1 Poisson 方程五点差分格式参考单如图所示：以(,)i j x y 为中心沿y 方向Taylor 展开：41)(),j u y o h +①41)(),j u y o h +②41(),u h21(),o h ③22(),o h ④(,),i j ij f x y R -=+(,),i j f x y -=○6 j+1考虑到边值条件(,)(,)u x y x y αΓ=，构成差分格式：11112212(,)2(,)(,)(,)2(,)(,)(,),(,)(,),i j i j i j i j i j i j i j u x y u x y u x y u x y u x y u x y f x y h h u x y x y α+-+-Γ⎧-+-+⎛⎫-+=⎪ ⎪⎨⎝⎭⎪=⎩○72.2 Poisson 方程九点差分格式由上式 ③ + ④ 得：11112212442221244222222122222(,)2(,)(,)(,)2(,)(,)(,)1(,)()12(,)(,)1(,)12i j i j i j i j i j i j h i j i j iji j i j i j u x y u x y u x y u x y u x y u x y u x y h h u u u x y h h o h x y u x y u x y u x y h h x y x y +-+--+-+=+⎡⎤∂∂=∆+++⎢⎥∂∂⎣⎦⎛⎫∂∂⎛⎫∂∂=∆+++- ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭422212222242222212122222(,)()12(,)(,)(,)1(,)()1212i j i j i j i j i j u x y h h o h x y f x y f x y u x y h h f x y h h o h x y x y ∂++∂∂⎛⎫∂∂∂+=--+-+ ⎪ ⎪∂∂∂∂⎝⎭○8 又()41122222211111112212311111(,)(,)2(,)(,)()1[(,)2(,)(,)2(,)2(,)(,)(,)2(,)(,)]()i j xx i j xx i j xx i j i j i j i j i j i j i j i j i j i j u x y u x y u x y u x y o h x y h u x y u x y u x y u x y u x y u x y h h u x y u x y u x y o h +-+++-++-+----∂-+=+∂∂=-+--++-++ 则得到：222222121121112112222221211212122222221112111211()(,)(210)(,)()(,)(210)(,)20()(,)(210)(,)(210)(,)()(,)()(,)i j i j i j i j i j i j i j i j i j h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y ---+--++-+++-++--++-+++-+--+-+2212222241222,12(,)(,)1(,)()12i j i j i j h hf x y f x y f x y h h o h x y ⎛⎫∂∂=--++ ⎪ ⎪∂∂⎝⎭○9 舍去截断误差得到逼近Poisson 方程的九点差分方程○10：()()2212,11,,11,1,11,11,11,122122212(,)[42]121(,)(,),12i j i j i j i j i j i j i j i j i j i j ij xx i j yy i j h h u x y u u u u u u u u u h h f h f x y h f x y -++-+---++-++-∆--+++++++''''=++考虑到边值条件(,)(,)u x y x y αΓ=，构成差分格式○11：()()2212,11,,11,1,11,11,11,122122212(,)[42]121(,)(,),12(,)(,),i j i j i j i j i j i j i j i j i j i j ijxx i j yy i j h h u x y u u u u u u u u u h h f h f x y h f x y u x y x y α-++-+---++-+Γ⎧+-∆--+++++++⎪⎪⎪''''=++⎨⎪⎪=⎪⎩3.格式求解3.1 Poisson 方程五点差分格式记122,1,j j j m j m j u u u u u --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，0.j n ≤≤ 矩阵格式改写为：11,11j j j j Du Cu Du f j m -+++=≤≤-，其中2221212222112122221121222112(1)111211112111121112m h h h h h h h C h h h h h h h -⎡⎤⎛⎫+-⎢⎥ ⎪⎝⎭⎢⎥⎢⎥⎛⎫⎢⎥-+- ⎪⎢⎥⎝⎭⎢⎥=⎢⎥⎢⎥⎛⎫⎢⎥-+- ⎪⎢⎥⎝⎭⎢⎥⎛⎫⎢⎥-+ ⎪⎢⎥⎝⎭⎣⎦，22222222(1)1111m h h D h h -⎡⎤-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦，10212212111(,)(,)(,)(,)1(,)(,)j j j j m j m j m j m f x y x y h f x y f f x y f x y x y h ---⎡⎤+Φ⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥+Φ⎢⎥⎣⎦， 可进一步写为：110222211(1)*(1).n n n n n n m u f Du C D u f D C D u f DC D u f Du D C -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦3.2 Poisson 方程九点差分格式记122,1,j j j m j m j u u u u u --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，0.j n ≤≤ 矩阵格式改写为：11,11j j j j Du Cu Du f j m -+++=≤≤-，其中2222121222222212121222222212121222221212(1)20()(210)(210)20()(210)(210)20()(210)(210)20()m h h h h h h h h h h C h h h h h h h h h h -⎡⎤+-⎢⎥-+-⎢⎥⎢⎥=⎢⎥-+-⎢⎥⎢⎥-+⎣⎦， 2222211222222212211222222212211222221221(1)(210)()()(210)()()(210)()()(210)m h h h h h h h h h h D h h h h h h h h h h -⎡⎤--+⎢⎥-+--+⎢⎥⎢⎥=⎢⎥-+--+⎢⎥⎢⎥-+-⎣⎦，22121022221211(,)(210)(,)(,)(,)(,)(210)(,)j j j j m j m j m j m f x y h h x y f x y f f x y f x y h h x y ---⎡⎤--Φ⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥+-Φ⎣⎦， 可进一步写为：110222211(1)*(1).n n n n n n m u f Du C Du f D C D u f DC D u f Du D C -------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦4.数值例子4.1 Poisson 方程五点差分格式计算如下问题：22220,01,01,(0,)sin cos ,(2,)(sin cos ),01,(,0),(,1)(sin1cos1),0 1.x x u u x y x y u y y y u y e y y y u x e u x e x ⎛⎫∂∂-+=<<<< ⎪∂∂⎝⎭=+=+≤≤==+<<其精确解为：(,)(sin cos ).x u x y e y y =+,11,1,,1,222222122112112()(,),i j i j i j i j i j i j u u u u u f x y h h h h h h -+-++=++++ 考虑到本例中h1=h2，则有2,11,1,,1,(,),4i j i j i j i j i j i j u u u u h f x y u -+-+++++=利用Gauss-Seidel 迭代方法对k=0,1,2,……,计算112,11,1,,11(,),41,2,....,1;1,2,...., 1.k k k k i j i j i j i j i j k ij u u u u h f x y u i m j n ++--+++++++==-=-表1 部分结点处的精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差绝对值图1 取h=1/4时所得的数值解曲线图2 取h=1/4时所得的误差曲线图3 取h=1/16时所得的数值解曲线图4 取h=1/16时所得的误差曲线图5 取h=1/64时所得的数值解曲线图6 精确解曲线图7 取h=1/64时所得的误差曲线4.2 Poisson 方程九点差分格式计算如下问题：22220,01,01,(0,)sin cos ,(2,)(sin cos ),01,(,0),(,1)(sin1cos1),0 1.x x u u x y x y u y y y u y e y y y u x e u x e x ⎛⎫∂∂-+=<<<< ⎪∂∂⎝⎭=+=+≤≤==+<<其精确解为(,)(sin cos ).x u x y e y y =+222222221212121112122222222121112111211211222211120()(,)12(,)()(,)(102)(,)()(,)()(,)()(,)(102)(,)(102)(,)(10i j i j i j i j i j i j i j i j i j h h u x y h h f x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h h u x y h ----++++--++=+++-+++++++-+-+2212)(,)i j h u x y +-考虑到本例中h1=h2，则有,11,1,,11,11,11,11,1,4(),20i j i j i j i j i j i j i j i j i j u u u u u u u u u -+-+--++-++-+++++++=利用Gauss-Seidel 迭代方法对k=0,1,2,……,计算1111,11,1,,11,11,11,11,11,4(),201,2,....,1;1,2,...., 1.k k k k k k k k i j i j i j i j i j i j i j i j k i j u u u u u u u u u i m j n ++++-+-+--++-++-++++++++==-=-表1 部分结点处的精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差绝对值表3 取不同步长时部分结点处数值解的最大误差图1 取h=1/4时所得的数值解曲线图2 取h=1/16时所得的数值解曲线图3 取h=1/64时所得的数值解曲线图4 取h=1/4时所得的误差曲线图5 取h=1/16时所得的误差曲线图6 取h=1/64时所得的误差曲线5.结论观察Poisson方程五点格式，方程以较快速度迭代收缩。

计算流体力学大作业——有限差分法解Poisson 方程五点格式解区域内Poisson 方程摘要：本文结合计算流体力学课上所学知识，采用数值解法中的有限差分法求解Poisson 方程（偏微分方程中椭圆型方程的一种），并用其五点格式采用高斯—塞德尔（Gauss-Seidel ）迭代求解。

并比较了数值近似解与真实解，以及不同步长情况下误差的大小，得到了一定的结论。

关键词：Poisson 方程 有限差分法 五点格式一、计算流体流体力学的特点计算流体力学中许多问题求解最终都会变成偏微分方程的求解，而在数学上，除了几种极少数情况外，要求出它们精确解是很难的。

计算机技术的发展使得这一难题的一很好地解决。

二、偏微分方程的种类2.1、 椭圆型偏微分方程椭圆型偏微分方程的一般形式为()(,)div c u au f x t -∇+= 其中：若12(,,,,)(,)n u u x x x t u x t ==，u ∇为u 的梯度，则其定义为 12,,,n u u x x x ⎡⎤∂∂∂∇=⎢⎥∂∂∂⎣⎦ 散度()div v 的定义为12()n div v v x x x ⎛⎫∂∂∂=+++ ⎪∂∂∂⎝⎭这样，()div c u ∇可以更明确地表示为1122()n n u u u div c u c c c x x x x x x ⎡⎤⎛⎫⎛⎫⎛⎫∂∂∂∂∂∂∇=+++⎢⎥ ⎪ ⎪ ⎪∂∂∂∂∂∂⎝⎭⎝⎭⎝⎭⎣⎦若c 为常数，则进一步化简为 22222212()n div c u c u c u x x x ⎛⎫∂∂∂∇=+++=∆ ⎪∂∂∂⎝⎭其中，∆又称为Laplace 算子。

这样椭圆型偏微分方程可以简单地写为22222212(,)n c u au f x t x x x ⎛⎫∂∂∂-++++= ⎪∂∂∂⎝⎭2.2、抛物型偏微分方程抛物型偏微分方程的一般形式为 ()(,)u d div c u au f x t t∂-∇+=∂ 根据上面叙述，若c 为常数，则该方程可以更简单地写为22222212(,)n u d c u au f x t t x x x ⎛⎫∂∂∂∂-++++= ⎪∂∂∂∂⎝⎭ 2.3、双曲型偏微分方程双曲型偏微分方程的一般形式为22()(,)u d div c u au f x t t∂-∇+=∂ 若c 为常数，则可以将该方程简化为2222222212(,)n u d c u au f x t t x x x ⎛⎫∂∂∂∂-++++= ⎪∂∂∂∂⎝⎭三类方程的直接的区别在于u 对t 的导数的阶次。

毕业论文文献综述信息与计算科学椭圆型偏微分方程的求解及其应用一、 前言部分微积分产生以后，人们就开始把力学中的一些问题，归结为偏微分方程进行研究。

早在18世纪初，人们已经将弦线振动的问题归结为弦振动方程，并开始探讨了它的解法。

随后，人们又陆续了解了流体的运动、弹性体的平衡和振动、热传导、电磁相互作用、原子核和电子的相互作用、化学反应过程等等自然现象的基本规律，把它们写成偏微分方程的形式，并且求出了典型问题的解答，从而能通过实践，验证这些基本规律的正确性，显示了数学物理方程对于认识自然界基本规律的重要性。

有了基本规律，人们还要利用这些基本规律来研究复杂的自然现象和解决复杂的工程技术问题，这就需要求出数学物理方程中的许多特定问题的解答。

随着电子计算机的出现及计算技术的发展，即使是相当复杂的问题，也有可能计算出解得足够精确的数值来，这对于预测自然现象的变化（如天气预报）和进行各种工程设计（如机械强度的计算）都有着很重要的作用[1]。

许多复杂的自然现象，其运动规律、过程和状态都是通过微分方程这种数学形式来描述的。

当我们研究只有一个自变量的运动过程时出现的微分方程称为常微分方程。

当一个微分方程除了含有几个自变量和未知数外，还含有未知数的偏导数时，称为偏微分方程[2]-[6]。

在偏微分方程中，偏导数自然是不可缺少的。

例如： ()(),,u ua x y f x y x y∂∂+=∂∂ （1.1.1） 拉普拉斯方程22232220u u uu x y z∂∂∂∆=++=∂∂∂（1.1.2） 热传导方程()222,,u u a f x t u t x ∂∂=+∂∂（1.1.3） 波动方程()2222,,u a u f t x y t∂=∆+∂（1.1.4）等都是偏微分方程。

其中，u 为未知数，a 为常数，(),a x y 、f 为已知函数。

偏微分方程的一般形式为()112,,,,,,,,0n n x x F x x x u u u ⋅⋅⋅⋅⋅⋅⋅⋅⋅= （1.1.5） 其中：F 为已知函数；12,,,n x x x ⋅⋅⋅为自变量；u 是关于这些自变量的未知数。