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

椭圆型方程

§1
差分逼近的基本概念
考虑二阶微分方程边值问题
d 2u Lu 2 qu f , a x b, dx u (a) , u (b) , (1.1) (1.2)
其中 q，f 为 [ a , b ] 上的连续函数， q 0, , 为给定常数. 将其分成等分，分点为
称
uh 收敛到边值问题的解 u .
对于差分方程
Lhvi fi , i 1, 2,3,L , N 1,
定义1.3
v0 vN 0 , 如果存在与网格 I h 及右端 fh 无关的常数
数 M 和 h0 , 使 || vh || M || f h ||R ,
0 h h0
称差分方程关于右端稳定.
第二章
椭圆形方程的有限差分法
有限差分法和有限元方法是解偏微分方程的两种主要数值
方法.
有限差分法：从定解问题的微分或积分形式出发，用数值微商或数值积分导出相应的线性代数方程组. 有限元方法：从定解问题的変分形式出发，用RitzGalerkin 方法导出相应的线性代数方程组，但基函数要按
特定方式选取.
取 x(1) x0 a, x(2) x1 , 得
2
(2.9) (2.10)
W (a) W ( x1 ) 2 qudx
d2 du hi 1 hi dx 2 ( p dx ) 12 i
d 3u 2 p O ( h ) dx 3 i
于是得逼近方程 (2.1)~(2.2) 的差分方程：
ui 1 ui ui ui 1 2 p 1 Lhui pi 1 i h h h h i i 1 i 1 i 2 2 i i 1, 2,, N 1 ui 1 ui qiui fi , hi hi 1 u0 , uN

五点差分格式求解椭圆型偏微分方程(解线性方程组方法)

五点差分格式求解椭圆型偏微分方程（解线性方程组方法）五点差分格式是一种常用的数值方法，用于求解椭圆型偏微分方程。

该方法将偏微分方程中的二阶导数项用差分近似替代，并将偏微分方程转化为一个线性方程组。

本文将介绍五点差分格式的推导过程，并使用该方法求解一个简单的椭圆型偏微分方程。

假设我们要求解的偏微分方程为：∂²u/∂x²+∂²u/∂y²=f(x,y)其中，u是未知函数，f(x,y)是已知函数。

我们将该方程离散化，将坐标(x,y)分别用h表示，将u(x,y)用U(i,j)表示，其中i和j分别表示x和y的离散位置。

我们可以使用中心差分近似来计算二阶导数，得到：∂²u/∂x²≈(U(i+1,j)-2U(i,j)+U(i-1,j))/h²∂²u/∂y²≈(U(i,j+1)-2U(i,j)+U(i,j-1))/h²将上述近似代入原方程，得到：(U(i+1,j)-2U(i,j)+U(i-1,j))/h²+(U(i, j+1) - 2U(i, j) + U(i, j-1)) / h² = f(ih, jh)整理上述方程，得到：U(i+1, j) + U(i-1, j) + U(i, j+1) + U(i, j-1) - 4U(i, j) = h² * f(ih, jh)该方程表示了U(i,j)与其相邻四个点的关系。

我们可以将整个区域离散化为一个网格，每个网格点都满足类似的方程。

离散化后的方程可以写成一个线性方程组的形式。

例如，在一个矩形区域内，我们将x轴和y轴的区间划分为n个小区间，即x轴上的取值为0, h, 2h, ..., nh；y轴上的取值为0, h,2h, ..., nh。

在该区域内，一共有(n-1)²个内部网格点。

我们可以将这些网格点按照其中一种顺序依次编号，从而将线性方程组表示为一个矩阵方程。

有限差分法求解偏微分方程MATLAB

南京理工大学课程考核论文课程名称：高等数值分析论文题目：有限差分法求解偏微分方程*名：**学号： 1成绩：有限差分法求解偏微分方程一、主要内容1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：22(,)()u uf x t t xαα∂∂-=∂∂其中为常数具体求解的偏微分方程如下：22001(,0)sin()(0,)(1,)00u u x t x u x x u t u t t π⎧∂∂-=≤≤⎪∂∂⎪⎪⎪=⎨⎪⎪==≥⎪⎪⎩2.推导五种差分格式、截断误差并分析其稳定性；3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析；4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。

有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：()2100000000()()()()()()()......()(())1!2!!n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+- (2-1)求解区域的网格划分步长参数如下：11k k k kt t x x h τ++-=⎧⎨-=⎩ (2-2) 2.1 古典显格式2.1.1 古典显格式的推导由泰勒展开公式将(,)u x t 对时间展开得 2,(,)(,)()()(())i i k i k k k uu x t u x t t t o t t t∂=+-+-∂ (2-3) 当1k t t +=时有21,112,(,)(,)()()(())(,)()()i k i k i k k k k k i k i k uu x t u x t t t o t t tuu x t o tττ+++∂=+-+-∂∂=+⋅+∂ (2-4)得到对时间的一阶偏导数1,(,)(,)()=()i k i k i k u x t u x t uo t ττ+-∂+∂ (2-5) 由泰勒展开公式将(,)u x t 对位置展开得223,,21(,)(,)()()()()(())2!k i k i k i i k i i u uu x t u x t x x x x o x x x x∂∂=+-+-+-∂∂ (2-6)当11i i x x x x +-==和时，代入式(2-6)得2231,1,1122231,1,1121(,)(,)()()()()(())2!1(,)(,)()()()()(())2!i k i k i k i i i k i i i i i k i k i k i i i k i i i iu uu x t u x t x x x x o x x x x u u u x t u x t x x x x o x x x x ++++----⎧∂∂=+-+-+-⎪⎪∂∂⎨∂∂⎪=+-+-+-⎪∂∂⎩(2-7) 因为1k k x x h +-=，代入上式得2231,,22231,,21(,)(,)()()()2!1(,)(,)()()()2!i k i k i k i k i k i k i k i ku uu x t u x t h h o h x xu u u x t u x t h h o h x x +-⎧∂∂=+⋅+⋅+⎪⎪∂∂⎨∂∂⎪=-⋅+⋅+⎪∂∂⎩ (2-8) 得到对位置的二阶偏导数2211,22(,)2(,)(,)()()i k i k i k i k u x t u x t u x t u o h x h+--+∂=+∂ (2-9) 将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得21112(,)(,)(,)2(,)(,)(,)()i k i k i k i k i k i k u x t u x t u x t u x t u x t f x t o h h αττ++---+⎡⎤-=++⎢⎥⎣⎦(2-10)为了方便我们可以将式(2-10)写成11122k kk k k k i i i i i i u u u u u f h ατ++-⎡⎤--+-=⎢⎥⎣⎦(2-11) ()11122k k k k k k i i i i i i u u uu u f h τατ++----+= (2-12)最后得到古典显格式的差分格式为()111(12)k k k k k i i i i i u ra u r u u f ατ++-=-+++ (2-13)2r h τ=其中，古典显格式的差分格式的截断误差是2()o h τ+。

有限差分法解偏微分方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

五点差分格式求解poisson方程

五点差分格式求解poisson方程五点差分格式是一种常用的数值求解偏微分方程的方法，特别适用于求解Poisson方程。

本文将介绍五点差分格式的原理和求解Poisson方程的具体步骤。

1. 五点差分格式原理：五点差分格式是一种离散化方法，将求解区域离散成网格，通过近似代替微分方程中的导数项，从而转化为代数方程组。

该方法利用了函数在离散的点上的函数值和导数值之间的关系，通过求解差分方程组来获得解的近似值。

2. Poisson方程的离散化：Poisson方程是一个二阶偏微分方程，可表示为：∇^2 u = f其中，u是未知函数，f是已知函数。

我们要在给定的边界条件下求解u。

为了使用五点差分格式，我们需要对方程进行离散化。

将求解区域离散为网格，假设网格的步长为h，则u的离散点可表示为(u_{i, j})，其中i和j分别表示网格点的行坐标和列坐标，通过对函数u在(u_{i, j})点上的二阶导数进行近似，Poisson方程可以转化为：(u_{i+1, j} + u_{i-1, j} + u_{i, j+1} + u_{i, j-1} - 4u_{i, j}) / h^2 = f_{i, j}这就是离散化后的差分方程。

3. 求解步骤：(1) 确定求解区域和边界条件：首先需要确定方程的求解区域和边界条件。

根据实际问题设定，将求解区域划分为适当大小的网格，并给出边界条件。

(2) 建立差分方程：将Poisson方程离散化为差分方程，根据离散点上的函数值和导数值之间的关系，建立差分方程组。

(3) 求解差分方程组：利用代数方法，将差分方程组转化为线性方程组，并求解该方程组获得离散点上的近似解。

(4) 进行后处理：根据求解的离散点上的近似解，可以进行相应的后处理操作，如可视化等，以获得更直观的结果。

4. 参考内容：(1) "Numerical Methods for Engineers and Scientists" by Amos Gilat and Vish Subramaniam：该书介绍了五点差分格式等常用的数值方法，并给出了详细的计算步骤和示例。

五点差分格式

《微分方程数值解》大作业（一）——椭圆型方程编程计算：采用五点差分格式求如下椭圆型方程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）需要对数值解与精确解作比较，以及不同步长选取下的误差比较。

椭圆型方程的有限差分法

第4章椭圆型方程的有限差分法§ 2 一维差分格式1、用积分插值法导出逼近微分方程的差分格式。

, d / du 、 du 上 , Lu=- (p )+r +qu=f,a<x<b, dx dx dxu(a)= a ,u(b)=如果系数p,q,r 以及右端f 光滑，则可用中矩形公式计算得解：考虑在［a,b ］内任一小区间［x ⑴，x (2)］，将上式在此区间上积分得 X ⑵d dU曲 d ；(p(x)d ；)dxW(x ⑴)W(x ⑵)x(2)dU(1) r ——dxxdx x ⑵dU 」 xr dxx (2) x (1) qudx X (2)X (1) qudxX (2)x (1)fdXX ⑵X (1) fdX其中，W(x) p吧特别地，取［X ⑴，X ⑵］为对偶单元［x i 1/2 , x i 1/2 ］，则W (x i 1/2)W(x i 1/2)x1/2dUr ——dxxi 1/2 qudxxi 1/ 2xi 1/ 2 fdx 。

Xi 1/ 2-J..将(1.2)改写成—dxW(x)再沿［x i 1/2, x i 1/2］积分，得U i U iX iW(x) dx , Xi 1p(x)矩形公式，得W 1/2a iU i U ih i a iNdx x1 p(x)］xi 1/ 2又qUdxXi 1/ 2h ih i2：d i U i , d ixi 1/ 2 x q(x)dxX1/2x1/2dU , r ——dx.U i 1 U i 1bi "T-2 h i h i 1x 1/2 xr(x)dxxi1/22ih i h i 1x1/ 2 f (x)dxxi 1/2(1.1) (1.2)利用中(1.3)(1.4)(1.5)(1.6)将(1.3) ~ (1.5)代入1.1),即得微分方程的差分格式a iU i U i 1h*h i h i 1)d i U b^Ui21-(h i h i 1) i则逼近阶为O(h 2)。

椭圆型微分方程

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

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

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

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

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

概括整个实验过程。

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

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

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

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

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

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

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

二阶椭圆偏微分方程实例求解(附matlab代码)

−2
�� 2 �� 2 �� 2 �� 2
�� +1/2,�� + ��−1/2,�� + �� ,�� −1/2 + �� ,�� +1/2 + ��
而 f 是（n − 1）维的列向量，具体如下： ��11 ��12 ⋮ �� = ��1,��−1 ��21 ⋮ ⋮ �� −1,��−1
2
三、求解过程
3.1 对系数矩阵的分析对上述模型的求解就是对线性方程组的求解。通过观察，我发现 P 是一个对角占优的矩阵，这不仅确定了解的唯一性，还保证了迭代法的收敛性。此外，还可以确定进行 LU 分解，若使用高斯消去法还可以省去选主元的工作。
二、问题分析与模型建立
2.1 微分方程上的符号说明 �� , �� = �� , �� = �� , �� = �� + �� , �� = �� − �� , �� = 1�� , �� = y 2 e y x 2 e x e xy y 2 x 2 1 e xy 2.2 课本上差分方程的缺陷课本上的差分方程为： �� − ��−1,�� −1,�� + �� ,�� −1 �� ,�� −1 + �� +1,�� +1,�� + �� ,�� +1 �� ,�� +1 = ��

五点差分格式课件

03
预测和决策
通过五点差分法对时间序列数据进行处理和分析，可以用于预测未来数
据的变化趋势，为决策提供依据和支持。
金融领域中的应用
股票价格分析
五点差分法在金融领域中可用于股票价格分析，通过对股票价格数据进行差分处理，分析股票价格的波动性和趋势性，从而预测股票价格的未来走势。
风险管理五点差分法在风险管理领域中可用于分析金融市场的波动性和相关性，从而评估投资组合的风险水平，为投资者提供决策依据。
与有限差分法（FDM）
五点差分法适用于求解二阶偏微分方程，而FDM适用于一阶偏微分方程。
02
五点差分法的原理
差分法的数学原理
差分法是一种数值分析方法，通过将连续函数离散化，将微分问题转化为差分问题，从而能够用离散的数值方法进行计算。
在差分法中，将连续的时间或空间离散成若干个小区间，并取其离散的代表点，分别用数值函数表示，然后构造离散的数值关系式近似代替微分方程或差分方程，从而可以求解微分方程或差分方程的数值解。
2
有限差分法
五点差分法是有限差分法的一种，通过在时间和空间上离散化偏微分方程，将连续的物理量离散化为有限个离散点上的数值，从而将微分问题转化为代数问题。
3
精度和稳定性
五点差分法具有较高的精度和稳定性，能够提供较为精确的数值解，并且对于不同的初边值条件和边界条件具有较强的适应性。
数据分析中的应用
与计算机科学技术的结合
五点差分法可以与计算机科学技术结合，如并行计算、云计算等，提高计算效率和精度。
在新领域的应用
金融领域
五点差分法可以应用于金融领域中的期权定价、风险管理等问题。
环境科学领域
五点差分法可以应用于环境科学领域中的流体动力学、气候模拟等问题。

五点差分格式

矩阵方程AU=K，K由边界条件所确定，解得
U = [U1 U2 U3 U4 U5 U6 U7 U8 U9]’=A-1K
= [55.7143 43.2143 27.1429 79.6429 70.0000 45.3571 112.8571 111.7857 84.2857]T
椭圆型方程的五点差分格式
椭圆型方程的五点差分格式
U2 4U2 U2
U2
U3 4U3
U3
U4
4U4 U4 U4
U5
U5 4U5 U5
U5
U6 U6 4U6
U6
U7
4U7 U7
U8
U8 4U8 U8
U9
U9 4U9
100 20 20 80 0 0 260 180 180
2u 2u x2 y2 0
在区域 (x, y) | 0 x 4,0 y 4
内的近似解，边界值为：
uy (x,0) 0,u(x,4) 180,0 x 4 u(0, y) 80,u(4, y) 0,0 y 4
取 h x y 1。

B

I

2 4

1
2I B
g 2g1 g2 2g3 g4 0 g6 2g7 g8 2g9 T
U U1 U2 U3 U4 U5 U6 U7 U8 U9 T
则矩阵方程为 AU 2hg
椭圆型方程的五点差分格式
加密网格，取 h = 0.5
椭圆型方程的五点差分格式
椭圆型方程的五点差分格式
定义向量
U [U1,1,U 2,1,,U M 1,1;U1,2 ,U 2,2 ,,U M 1,2; U1,M 1,U 2,M 1,,U M 1,M 1]T

椭圆型方程的差分解法

椭圆型方程的差分解法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方程五点格式，方程以较快速度迭代收缩。

泊松方程-椭圆型方程的五点格式

称为半整数点，则由节点 a x0 x 1 x 3 x
2 2 i 1 2
x
N
1 2
又构成[a, b]的一个网格剖分，称为对偶剖分.
其次用差商代替微商将方程(2.1)在节点xi离散化，为此，对充分光滑的解u，由Taylor展式可得 u ( xi 1 ) u ( xi 1 ) hi hi 1 hi 1 hi d 2u du [ ]i [ 2 ]i o(h 2 ), (2.3) dx 2 dx 其中[ ]i 表示括号内函数xi点取值。
hi 1 hi 2 du du d 3u ([ p ] 1 [ p ] 1 ) [ p 3 ]i o(h 2 ), hi hi 1 dx i 2 dx i 2 12 dx hi 1 hi d 2 hi 1 hi d du du d 3u [ ( p )]i [ 2 ( p )]i [ p 3 )]i o(h 2 ), (2.6) dx dx 4 dx dx 12 dx
两族直线的交点(ih1 , jh2 )称为网点或节点，记为( xi , y j )或(i, j ). 说两个节点( xi , y j )和( xi , y j )是相邻的, 如果 xi xi yi yi 1或 i i j j 1. h1 h2
xi
利用中矩形公式，得 W ai [ 又 1 hi
i
1 2
ai
ui ui 1 , hi
(2.16) (2.17 ) (2.18)

xi
xi 1
dx 1 ] . p( x)

x
i
1 2
x
qudx
i
1 2
hi hi 1 d i ui , 2

椭圆方程差分格式

由Guass公式有
D
D
u u udxdy D n ds ( l l l ) n ds D l1 2 3 4
从上面的公式把边界积分转化为在四条边上的积分，于是积分分成段进行，也即转化为定 4 积分。
对于在边l1上，因为在矩形边1上的外法向就是的负 l y 方向，弧长的微分 dx，于是对此定积分用矩 ds 形公式近似计算，并且用差商代替微商，得到
u u x (b( x, y) x )dxdy [b( xi , y j 12 ) x ( xi , y j 12 ) Dij u b( xi , y 1 ) ( xi , y 1 )]h j j 2 x 2
c( x, y)udxdy c u hk, f ( x, y)dxdy f
１：直接差分方法
1 1 (aij xuij ) 2 y (aij y uij ) cijuij f ij 2 x h h
２：有限体积法（积分差分方法）
u x (a( x, y) x )dxdy Dij
y
u u [a( x 1 , y ) ( x 1 , y ) a( x 1 , y ) ( x 1 , y )]dy i i x i 2 x i 2 2 2 y 1
j 2
j
1 2
对上面定积分利用梯形公式有
u u x (a( x, y) x )dxdy [a( xi 12 , y j ) x ( xi 12 , y j ) Dij u a( x 1 , y j ) ( x 1 , y j )]k i x i 2 2
1.1:五点差分格式
对于poisson方程，考虑在内部节点 xi，y j）取值，（于是有 2u 2u [ 2 ]ij [ 2 ]ij [ f ( x, y )]ij x y

椭圆型偏微分方程的求解及其应用[文献综述]

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

早在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 是关于这些自变量的未知数。

偏微(13)椭圆型方程的差分方法

2
S1 e u0 e u0 e u0 e u0
1 4 2 4u0 h u0 h 2 D 4 u0 h6 12
2

1 4 2 S2 4u0 2h u0 h 4 D 4 u0 h6 6
其中D x , y | 0 x , y 1 .
1 取特殊的网格 h k 。此时网格点分布见图5.3， 3
1.1 五点差分格式
在内点P , P2 , P3和P4 1 上用差分格式（）， 1.6
在其余点，即边界点取边界条件
1 x 2 y 2 a x, y log
1 u 1 u r u r r r 2 2 f r , r r
2
(1.8)
r
y x y , tan . 域0 r ,0 2 . x
2 2
方程( 8) 1. 的系数当r 0时具有奇异性，因此为了选出我们感兴趣的解，需补充u在r 0处有界的条件，可设u满足
1 u xi h, y j 2u xi , y j u xi h, y j （1.2） h2 Nhomakorabea

2u 2u u 2 2 f ( x , y（） ) 1.1 x y

1 u x i , y j k 2u x i , y j u x i , y j k 2 k （1.3） 2 2 4 4 u k 2 4 u xi ,1 4 u xi ,2 y y ij 24 y
h uij ui 1, j 2uij ui 1, j h2 ui , j 1 2uij ui , j 1 k2 f ij ,

五点差分格式求解椭圆型偏微分方程(解线性方程组方法)

五点差分格式求解椭圆型偏微分方程华南理工大学数创班张华富【1】简述采用五点差分格式求解椭圆型偏微分方程，通过求解线性方程组AU=F 得到数值解；对数值解与精确解进行误差分析，分析收敛性及收敛阶。

难点：求出块三对角系数矩阵A ，以及右端项F ；【2】例题编程计算：采用五点差分格式求如下椭圆型方程2222uu x y f (x,y),(x,y);∂∂∂∂--=∈Ω其中f (x,y)、Ω及边值条件为：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 xy =+【3】代码(1)q1.m%求解例题clearclc M=[5 10 20 30 40 50 60 70 80 90]; %网格剖分步数N=M;f_jingque=@(x,y)log(x^2+y^2); %方程的精确解f=@(~,~)0; %方程的右端项u_0=@(y)log(1+y^2); %边界条件u_M=@(y)log(4+y^2);v_0=@(x)2*log(x);v_N=@(x)log(x^2+1);for i=10:-1:1[h(i),n(i)]=Question(1,2,0,1,f_jingque,f,u_0,u_M,v_0,v_N,M(i),N(i));%对每一个剖分，返回其误差向量的无穷范数n与步长hendn=log(n);h=log(h);p=polyfit(h,n,1) %拟合曲线表达式figureplot(h,n,'*',h,polyval(p,h),'b-') %作图，画出拟合曲线legend('原始数据点','拟合曲线');（2）调用函数 Question.mfunction [h,n]=Question(a,b,c,d,f_jingque,f,u_0,u_M,v_0,v_N,M,N) %a,b,c,d为边界点，f_jingque为偏微分方程的精确解，f为方程右端项，u及v为边界条件，M、N为网格剖分步数h1=(b-a)/M;h2=(d-c)/N; %x,y方向的步长h=sqrt(h1^2+h2^2);for i=1:M-1x(i)=a+i*h1;V_0(i)=v_0(x(i));V_N(i)=v_N(x(i));endfor j=1:N-1y(j)=c+j*h2;U_0(j)=u_0(y(j));U_M(j)=u_M(y(j));endalpha_1=1/(h1^2); %五点差分格式的系数alpha_2=1/(h2^2);alpha_3=1/(h1^2);alpha_4=1/(h2^2);alpha_0=alpha_1+alpha_2+alpha_3+alpha_4;B=zeros(M-1,M-1); %求块三对角矩阵A中的主对角线元素B，B为三对角矩阵B(1,1)=alpha_0;B(1,2)=-alpha_1;B(end,end)=alpha_0;B(end,end-1)=-alpha_3;for i=2:M-2B(i,i)=alpha_0;B(i,i-1)=-alpha_3;B(i,i+1)=-alpha_1;endD=-alpha_4*eye(M-1);C=-alpha_2*eye(M-1);A=release(B,D,C,M,N); %生成块三对角矩阵A%线性方程组AU=F的右端项F的计算F=zeros((M-1)*(N-1),1);%四个顶角的FF(1)=f(x(1),y(1))+alpha_3*U_0(1)+alpha_4*V_0(1);F(M-1)=f(x(M-1),y(1))+alpha_1*U_M(1)+alpha_4*V_0(M-1);F((M-1)*(N-2)+1)=f(x(1),y(N-1))+alpha_2*V_N(1)+alpha_3*U_0(N-1);F((M-1)*(N-1))=f(x(M-1),y(N-1))+alpha_1*U_M(N-1)+alpha_2*V_N(M-1); for k=2:M-2 %第一层F(k)=f(x(k),y(1))+alpha_4*V_0(k);endfor k=(M-1)*(N-2)+2:(M-1)*(N-1) -1 %第N-1层F(k)=f(x(k-(M-1)*(N-2)),y(N-1))+alpha_2*V_N(k-(M-1)*(N-2));endl=1;k=1+(M-1)*l; %第一列while k<(M-1)*(N-2)+1F(k)=f(x(1),y(l+1))+alpha_3*U_0(l+1);l=l+1;k=1+(M-1)*l;endl=1;k=M-1+(M-1)*l; %第M-1列while k<(M-1)*(N-1)F(k)=f(x(M-1),y(l+1))+alpha_1*U_M(l+1);l=l+1;k=M-1+(M-1)*l;endU_shuzhi2=A\F; %左除求数值解向量U_shuzhi1=reshape(U_shuzhi2,[M-1,N-1]); %将向量转换为矩阵Accurate1=zeros(N-1,M-1);for j=1:N-1for i=1:M-1Accurate1(i,j)=f_jingque(x(i),y(j)); %精确解矩阵endendAccurate2=reshape(Accurate1,[(M-1)*(N-1),1]); %精确解矩阵化为按自然顺序排序的向量error=Accurate2-U_shuzhi2; %精确解与数值解的误差向量n=norm(error,inf); %误差向量的无穷范数plot_123(x,y,U_shuzhi1,Accurate1); %作图，画出数值解、精确解与误差矩阵的图像（3）调用函数 release.m%生成块三对角矩阵Afunction A=release(B,D,C,M,N)A=zeros((M-1)*(N-1));A(1:M-1,1:M-1)=B;A(1:M-1,M:2*(M-1))=C;A((M-1)*(N-2)+1:(M-1)*(N-1),(M-1)*(N-3)+1:(M-1)*(N-2))=D;A((M-1)*(N-2)+1:(M-1)*(N-1),(M-1)*(N-2)+1:(M-1)*(N-1))=B;for k=2:N-2A((M-1)*(k-1)+1:(M-1)*k,(M-1)*(k-1)+1:(M-1)*k)=B;A((M-1)*(k-1)+1:(M-1)*k,(M-1)*(k-2)+1:(M-1)*(k-1))=D;A((M-1)*(k-1)+1:(M-1)*k,(M-1)*k+1:(M-1)*(k+1))=C;end（4）调用函数 plot_123.mfunction plot_123(x,y,U_shuzhi,Accurate)figuresubplot(1,3,1)surf(x,y,U_shuzhi'); %按自然顺序排序输出，需将矩阵转置xlabel('x')ylabel('y')title('数值解图像')subplot(1,3,2)surf(x,y,Accurate');xlabel('x')ylabel('y')title('精确解图像')subplot(1,3,3)surf(x,y,(U_shuzhi-Accurate)');xlabel('x')ylabel('y')title('误差图像')【4】运行m文件q1.m，得到取不同步长时的数值解图像与精确解图像，以及误差图像；通过对步长与误差的拟合，可以得到收敛性以及收敛阶；比如，当M=N=30时，对应的图像如下图所示：拟合h与n，得到图像如下：得到的拟合函数为一次函数，且其系数为：可见，该方程的五点差分格式收敛，且收敛阶约为二阶，这与理论收敛阶数相一致。

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

function z=f(x,y)
z=(pi^2-1)*exp(x)*sin(pi*y);
function z=g1x(x)
z=0;
function z=g2x(x)
z=0;
function z=g1y(y)
z=sin(pi*y);
function z=g2y(y)
z=exp(2)*sin(pi*y);
m=length(x);
n=length(y);
c1=zeros((n-2)*(m-2),1);
c1(1)=-1/h1^2*g1y(y(2))-1/h2^2*g1x(x(2));
for r=2:n-2
c1(r)=-1/h1^2*g1y(y(r+1));
end
for r=2:m-3
c1((r-1)*(n-2)+1)=-1/h2^2*g1x(x(r+1));
c1(r*(n-2))=-1/h2^2*g2x(x(r+1));
end
c1((m-3)*(n-2)+1)=-1/h2^2*g1x(x(m-1))-1/h1^2*g2y(y(2));
for r=2:n-2
c1((m-3)*(n-2)+r)=-1/h1^2*g2y(y(r+1));
end
c1((n-2)*(m-2))=-1/h1^2*g2y(y(n-1))-1ቤተ መጻሕፍቲ ባይዱh2^2*g2x(x(m-1));
x=a:h1:b;
y=c:h2:d;
m=length(x);
n=length(y);
C=tri1(-1/h1^2,2*(1/h1^2+1/h2^2),-1/h1^2,n-2);
D=-1/h1^2*eye(n-2);