五点差分法解椭圆型偏微分方程
偏微分方程(椭圆型)数值解2.4-5
ai−1, j
=
⎛ ⎜⎝
Ai
−
1 2
,
j
−
h1 2
Cij
⎞ ⎟⎠
h1−2
,
ai, j+1
=
⎛ ⎜⎝
Ai
∂n q1q2
∂n q2q3
∂n q3m4
• = ∫∫ f ( x, y)dxdy m1
G0
p1
（4.6）
对于后四项仿照公式（4.3）的方法离散化，例如
∫ ∫ ( ) m1q1
∂u ds ≈ ∂n
m1q1 ⋅ p0 p1
u1 − u0
,
( ) q2q1
∂u ds ≈ ∂n
q1q2 ⋅ p0 p2
u2 − u0
考虑矩形网格：h1和h2分别为x和y方向的步长，Gh为网格 节点集合，Gh为网格界点集合，Gh = Gh ∪ Γh 。 总假设Gh是联通
的，即 ∀P, P ∈ Gh 必有一串 Pi ∈Gh (i = 1, 2, L, m −1) 可将 P, P
排列成顺序序列：P, P1, P2 , L, Pm , P 使前后两点互为相邻节点。
G0的面积:
m
q7
(G
=
0)
q1 ，i
⎛ = 12× ⎜⎜⎝
= 1,
1× 2
2,
3 3
L, 6。
h
×
1 2
h
⎞ ⎟⎟⎠
五点差分法解椭圆型偏微分方程
用差分法解椭圆型偏微分方程-(Uxx+Uyy)=(pi*pi-1)e^xsin(pi*y) 0<x<2; 0<y<1U(0,y)=sin(pi*y),U(2,y)=e^2sin(pi*y); 0=<y<=1U(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(pi*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,40,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,瞧程序)有关;固定网格数,随着设定的迭代前后两次差值变小,误差由大比变小,中间有一个最小值,随着又增大一点,最后趋于稳定。
椭圆型方程的差分方法
椭圆型方程的差分方法差分方法是一种数值计算方法，使用近似的差商来表示微分方程。
椭圆型方程是一类常见的偏微分方程，具有重要的数学和物理应用。
在本文中，我们将介绍椭圆型方程的差分方法，并讨论其优点和缺点。
一、椭圆型方程的差分近似L[u]=-∂(p∂u/∂x)/∂x-∂(q∂u/∂y)/∂y+r(x,y)u=f(x,y)其中，L[u]是一个偏微分算子，u(x,y)是未知函数，p(x,y)，q(x,y)，r(x,y)，f(x,y)是已知函数。
椭圆型方程的解通常在一个区域Ω上求解。
差分方法的主要思想是用网格来离散化区域Ω，将连续的偏微分方程转化为离散的代数方程。
对于椭圆型方程，我们可以选择矩形网格，其中Ω可以被划分为N*M个小矩形，并且网格的步长为Δx和Δy。
假设我们要在网格点(xi, yj)处求解未知函数的值uij，其中i和j分别表示网格的行索引和列索引。
我们可以使用中心差分法来近似x和y方向的偏导数，从而得到离散形式的椭圆型方程:L[u] ≈ -(p(xi+1/2, yj)(ui+1,j - ui-1,j)/Δx^2 + p(xi,yj+1/2)(ui,j+1 - ui,j-1)/Δy^2) + q(xi,yj)uij = f(xi,yj)其中，p(xi+1/2, yj)和p(xi, yj+1/2)分别表示在(xi+1/2, yj)和(xi, yj+1/2)处的系数。
可以通过有限差分方式计算出这些系数。
将上述公式在每个网格点(xi, yj)处形成一个方程，从而得到一个线性方程组。
通过求解这个线性方程组，我们可以得到网格点上的未知函数值。
二、椭圆型方程差分方法的优点和缺点差分方法是一种简单有效的数值计算方法，具有以下优点：1.可以处理任意形状的区域Ω：差分方法可以适应不规则网格和复杂区域，因此适用于各种几何形状的椭圆型方程求解。
2.数值稳定性：差分方法可以确保数值解的稳定性，避免数值上的不稳定问题。
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
第九章偏微分方程差分方法汇总
第9章 偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。
由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。
偏微分方程的数值方法种类较多，最常用的方法是差分方法。
差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。
本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。
9.1椭圆型方程边值问题的差分方法9.1.1 差分方程的建立最典型的椭圆型方程是Poisson （泊松）方程G y x y x f yux u u ∈=∂∂+∂∂-≡∆-),(),,()(2222 （9.1）G 是x ,y 平面上的有界区域，其边界Γ为分段光滑的闭曲线。
当f (x ,y )≡0时，方程（9.1）称为Laplace(拉普拉斯）方程。
椭圆型方程的定解条件主要有如下三种边界条件第一边值条件 ),(y x u α=Γ （9.2） 第二边值条件),(y x nuβ=∂∂Γ （9.3） 第三边值条件 ),()(y x ku nuγ=+∂∂Γ （9.4） 这里，n 表示Γ上单位外法向，α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数，k (x,y )≥0。
满足方程（9.1）和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。
用差分方法求解偏微分方程，就是要求出精确解u (x ,y )在区域G 的一些离散节点（x i ，y i ）上的近似值u i ,j ≈(x i ，y i ）。
差分方法的基本思想是，对求解区域G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。
设G ={0<x <a , 0<y <b }为矩形区域，在x ,y 平面上用两组平行直线x =ih 1, i =0,1,…,N 1, h 1=a /N 1 y =jh 2, j =0,1,…,N 2, h 2=b /N 2将G 剖分为网格区域，见图9-1。
椭圆型方程差分法
(4)
令
2 1 1 2 1 A 1 1 2 ( N 1)( N 1)
系数矩阵A是不可约对角占优阵 A 0
解存在唯一,或直接求A的特征值。
8
习题：计算矩阵
A=
2 1 1 2 1 1 1 2 ( N 1)( N 1)
x
2． Poisson方程五点差分格式
u f u in
其中
(0, a ) (0, b )
建立目标点： a b y h k x 一方向步长： I 1 ； 一方向步长： J 1
21
得
( xi , y j )
1 i I,1 j J
xi ih, y j jk
返回
2
1） 数值计算是否必要？
T '( x) T '(0) f (u )du
x 0
x
T ( x) T '(0) f ( s )ds du
0
0
u
T '(0) x
T '(0) x
x 0 x
x
x
0
u
0
f ( s)dsdu
s
f ( s )duds
从而得到迭代法：
Mxk 1 Nxk b
xk 1 M 1Nxk M 1b Sxk Tb
(*1)
18
阻尼迭代法 (Damped Iterative Method)
k1 Sxk M1b x k1 (1)xk [M1N(1)I]xk M1b (*2) xk1 x
五点差分格式求解椭圆型偏微分方程(解线性方程组方法)
五点差分格式求解椭圆型偏微分方程（解线性方程组方法）五点差分格式是一种常用的数值方法，用于求解椭圆型偏微分方程。
该方法将偏微分方程中的二阶导数项用差分近似替代，并将偏微分方程转化为一个线性方程组。
本文将介绍五点差分格式的推导过程，并使用该方法求解一个简单的椭圆型偏微分方程。
假设我们要求解的偏微分方程为：∂²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)²个内部网格点。
我们可以将这些网格点按照其中一种顺序依次编号，从而将线性方程组表示为一个矩阵方程。
五点差分格式求解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）需要对数值解与精确解作比较，以及不同步长选取下的误差比较。
椭圆型微分方程
数学与计算科学学院实验报告
实验项目名称椭圆型方程数值解
所属课程名称微分方程数值解法
实验类型验证
实验日期
班级信计0902
学号
姓名
成绩
附录1：源程序
附录2：实验报告填写说明
1．实验项目名称：要求与实验教学大纲一致。
2．实验目的：目的要明确，要抓住重点，符合实验教学大纲要求。
3．实验原理：简要说明本实验项目所涉及的理论知识。
4．实验环境：实验用的软、硬件环境。
5．实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容。
概括整个实验过程。
对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。
对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设
计思路和设计方法，再配以相应的文字说明。
对于创新性实验，还应注明其创新点、特色。
6．实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。
7．实验结论（结果）：根据实验过程中得到的结果，做出结论。
8．实验小结：本次实验心得体会、思考和建议。
9．指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。
数学中的椭圆型偏微分方程
数学中的椭圆型偏微分方程在数学领域中，椭圆型偏微分方程是一类重要的方程类型。
它在物理学、工程学和计算机科学等各个领域都有广泛的应用。
本文将介绍椭圆型偏微分方程的定义、性质和求解方法，从而帮助读者更好地理解和应用这一方程类型。
一、椭圆型偏微分方程的定义椭圆型偏微分方程是指具有标准形式的二阶偏微分方程，其中二次项系数的行列式不为零。
一般而言，椭圆型偏微分方程可以表示为：∑[i,j=1 to n] {aij(x) ∂²u/∂xi ∂xj} + ∑[i=1 to n] bi(x) ∂u/∂xi + cu = f其中，a_ij、b_i、c、f是相关系数或函数；u是未知函数，表示问题的解；x_1，x_2，…，x_n是自变量。
二、椭圆型偏微分方程的性质1. 正定性：椭圆型偏微分方程的二次项系数矩阵是正定矩阵。
这意味着椭圆型方程的解在定义域上满足一定的正定性条件。
2. 内部渐进性：椭圆型方程的解在区域的内部是光滑且渐进的。
3. 边界条件：椭圆型方程需要通过边界条件来获得唯一解。
常见的边界条件包括：泊松方程中的迪里切特边界条件和诺依曼边界条件。
三、椭圆型偏微分方程的求解方法1. 分离变量法：分离变量法是椭圆型偏微分方程求解的一种常见方法。
通过假设解可以表示为各个自变量分量的乘积形式，然后将未知函数与其各个自变量的分量进行分离，最终得到一个由各自变量分量的常微分方程组成的代数方程。
2. 特征线法：特征线法适用于一类特殊的椭圆型偏微分方程。
通过求解特征方程，我们可以找到解的参数化表示，从而将原方程化为一个更简单的常微分方程。
3. 有限差分法：有限差分法是一种通过在空间和时间上离散化方程来数值求解椭圆型偏微分方程的方法。
通过将偏微分方程转化为差分方程，可以用迭代方法求解离散问题。
四、椭圆型偏微分方程的应用1. 热传导方程：热传导方程可以描述物体内部温度分布随时间变化的情况。
通过求解热传导方程，我们可以研究热量在不同材料中的传导行为。
五点差分格式求解椭圆型偏微分方程(解线性方程组方法)
五点差分格式求解椭圆型偏微分方程华南理工大学 数创班 张华富【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，得到图像如下：得到的拟合函数为一次函数，且其系数为：可见，该方程的五点差分格式收敛，且收敛阶约为二阶，这与理论收敛阶数相一致。
五点差分格式
矩阵方程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方程五点格式，方程以较快速度迭代收缩。
微分方程数值解实验报告
微分方程数值解法课程设计报告班级：姓名：学号：成绩：2017年 6月 21 日摘要自然界与工程技术中的很多现象，可以归结为微分方程定解问题。
其中，常微分方程求解是微分方程的重要基础内容。
但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。
，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。
同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。
关键词：微分方程数值解、MATLAB目录摘要 (2)目录 (3)第一章常微分方程数值解法的基本思想与原理 (4)1.1常微分方程数值解法的基本思路 (4)1.2用matlab编写源程序 (4)1.3常微分方程数值解法应用举例及结果 (5)第二章常系数扩散方程的经典差分格式的基本思想与原理 (6)2.1常系数扩散方程的经典差分格式的基本思路 (6)2.2 用matlab编写源程序 (7)2.3常系数扩散方程的经典差分格式的应用举例及结果 (8)第三章椭圆型方程的五点差分格式的基本思想与原理 (10)3.1椭圆型方程的五点差分格式的基本思路 (10)3.2 用matlab编写源程序 (10)3.3椭圆型方程的五点差分格式的应用举例及结果 (12)第四章总结 (12)参考文献 (12)第一章常微分方程数值解法的基本思想与原理1.1常微分方程数值解法的基本思路常微分方程数值解法(numerical methods forordinary differential equations)计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.1.2用matlab编写源程序龙格库塔法：M文件:function dx=Lorenz(t,x)%r=28,sigma=10,b=8/3dx=[-10*(x(1)-x(2));-x(1)*x(3)+28*x(1)-x(2);x(1)*x(2)-8*x(3)/3];运行程序：x0=[1,1,1];[t,y]=ode45('Lorenz',[0,100],x0);subplot(2,1,1) %两行一列的图第一个plot(t,y(:,3))xlabel('time');ylabel('z');%画z-t图像subplot(2,2,3) %两行两列的图第三个plot(y(:,1),y(:,2))xlabel('x');ylabel('y'); %画x-y图像subplot(2,2,4)plot3(y(:,1),y(:,2),y(:,3))xlabel('x');ylabel('y');zlabel('z');%画xyz图像欧拉法：h=0.010;a=16;b=4;c=49.52;x=5;y=10;z=10;Y=[];for i=1:800x1=x+h*a*(y-x);y1=y+h*(c*x-x*z-y);z1=z+h*(x*y-b*z);x=x1;y=y1;z=z1;Y(i,:)=[x y z];endplot3(Y(:,1),Y(:,2),Y(:,3));1.3常微分方程数值解法的应用举例及结果应用举例：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=--=--=)()()()()()()()()())()(()(t bz t y t x dt t dz t z t x t y t rx dt t dy t y t x a dt t dx a=10,b=8/3,0<r<+∞,当1<r<24.74时，Lorenz 方程有两个稳定的不动点c()1(-r b ，)1(-r b ，r-1)和c '（-)1(-r b ，-)1(-r b ，r-1），一个稳定的不动点0=（0，0，0），当r>24.74时，c 和c '都变成不稳定的，此时存在混沌和奇怪吸引子。
椭圆型偏微分方程的求解及其应用[含论文、综述、开题-可编辑]
![椭圆型偏微分方程的求解及其应用[含论文、综述、开题-可编辑]](https://img.taocdn.com/s3/m/8669e8a5a21614791611280a.png)
设计（20 届）椭圆型偏微分方程的求解及其应用所在学院专业班级信息与计算科学学生姓名学号指导教师职称完成日期年月摘要：本文叙述了椭圆型偏微分方程的历史背景，阐述了相关概念，如什么是偏微分方程，椭圆型偏微分方程以及几种定解问题的概念。
弹性力学中的平衡问题，位势场问题，热传导中的温度分布等实际应用问题都可用椭圆型方程的定解问题来描述。
本文还讨论了求解椭圆型偏微分方程的定解问题的几种基本方法，如分离变量法、积分变换法、差分法，最后综述了这三种方法的适用性和特点。
关键字：偏微分方程；椭圆型；分离变量法；积分变换法；差分法Solution of Elliptic Partial Differential Equation and ItsApplicationAbstract: This thesis describes the historical background of elliptic partial differential equation and the related concepts, such as what partial differential equation and elliptic partial differential equation are and several concepts of the solution of problems. The balance of elasticity, the potential field problems and the temperature distribution of heat conduction in the practical application are available to the solution of elliptic equation to describe the practical problems. This thesis also discusses several basic ways to solve the solution of problems of the elliptic partial differential equation, for instance, the method of separation of variables, integral transformation method and difference method. And at the end of this thesis, it summarizes the applicability and features of the three methods above.Key Words: partial differential equation; elliptic; the method of separation of variables; integral transformation method; difference method目录1 引言 (1)2 基本概念的介绍 (2)2.1 偏微分方程的基本概念 (2)2.1.2 定解条件和定解问题 (3)2.2 两个自变量的二阶线性偏微分方程的分类与化简 (3)2.3 典型方程 (5)3 椭圆型偏微分定解问题的几种基本解法 (6)3.1 分离变量法 (6)3.1.1 预备知识 (6)3.1.2 分离变量法求解定解问题的具体步骤 (7)3.1.3 具体应用（用分离变量法求解） (7)3.2 积分变换法 (9)3.2.1 傅里叶积分变换 (9)3.2.2 具体应用（用积分变换法求解） (11)3.3 差分法 (13)3.3.1 化微分方程为差分方程 (13)3.3.2 边值问题的差分逼近 (16)3.3.3 差分解的存在、唯一性和收敛性 (18)3.3.4 椭圆型差分方程的求解——逐次超松弛法 (19)3.4 总结 (21)4 致谢 (22)参考文献 (23)1 引言数学物理方程主要指从物理学及其他各门自然科学、技术科学中所产生的偏微分方程（有时也包括积分方程、微分积分方程等），它们反映了有关的未知变量关于时间的导数和关于空间变量的导数之间的制约关系[1]。
泊松方程-椭圆型方程的五点格式
称为半整数点，则由节点 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
偏微(13)椭圆型方程的差分方法
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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
用差分法解椭圆型偏微分方程
-(Uxx+Uyy)=(pi*pi-1)e^xsin(pi*y) 0<x<2; 0<y<1
U(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));
end
for(i=1:n)
for(j=1:m)
f(i,j)=(pi*pi-1)*exp(x(j))*sin(pi*y(i));
end
end
t=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;
end
end
t(i,j)=sqrt(t(i,j));
if(k>kmax)
break;
end
if(max(max(t))<ep)
break;
end
end
for(i=1:n+1)
for(j=1:m+1)
p(i,j)=exp(x(j))*sin(pi*y(i));
e(i,j)=abs(u(i,j)-exp(x(j))*sin(pi*y(i)));
end
End
在命令窗口中输入:
[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,40,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=2164
surf(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,瞧程序)有关;固定网格数,随着设定的迭代前后两次差值变小,误差由大比变小,中间有一个最小值,随着又增大一点,最后趋于稳定。
也许可以去研究一下那个误差最小的地方或者研究趋于稳定时的临界值。