中心差分解两点边值问题

一．题目

用中心差分格式计算如下两点边值问题

dxduxsinx1xuxe2x1sinx1xxx1,1x2e dxdxu(1)0,u(2)2已知其精确解为

u(x)x(x1)

二．理论

作为模型，考虑两点边值问题：

Luddudu(p(x))rq(x)uf(x),axb…………（1.1） dxdxdxu(a),u(b)…………（1.2）

1假定pC[a,b],p(x)pmin0,r,q,fC[a,b],,是给定的常数。

1. 建立差分格式 （1）.区域网格剖分 首先取N1个节点：

ax0x1将区间I[a,b]分成N个小区间：

xixNb,

Ii:xi1xxi,i1,2,N.

i于是得到区间I的一个网格剖分。记hixixi1，称hmaxhi为网格最大步长。用Ih表示网格内点x1，x2，

，xN1的集合，Ih表示内点和界点x0a,xNb的集合。

取相邻节点xi1,xi的中点xi121(xi1xi)(i1,2,2xN12,N)，称为半整数点。则由节点

ax0x1x322xi12xNb

又构成[a,b]的一个网格剖分，称为对偶剖分。 （2）.微分方程的离散，建立相应差分格式

用差商代替微商，将方程（1.1）在内点xi离散化.注意对充分光滑的u，由Taylor展式有

u(xi1)u(xi1)duhi1hid2u[]i[2]iO(h2)………（1.3）

hihi1dx2dx.

.

u(xi)u(xi1)hi2dud3up(x1)[p]1[p3]1O(h3)ihidxi224dxi22[phdudu]1i[p3]iO(h3)（1.4）dxi224dx23

u(xi1)u(xi)hi12d3udup(x1)[p]1[p3]iO(h3)

ihi1dxi224dx2 ………（1.5）

由（1.5）减（1.4），并除以

hihi1，得 2u(xi1)u(xi)u(xi)u(xi1)2[p(x1)p(x1)]iihihi1hhi1i22…………（1.6）

hi1hid3u2dudu([p]1[p]1)[p3]iO(h2)hihi1dxi2dxi212dxhi1hid2hi1hid3uddudu[(p)]i[2(p)]i[p3]iO(h2)dxdx4dxdx12dx令pi12p(xi12),rir(xi),qiq(xi),fif(xi),则由（1.3）（1.6）知，边值问题的解u(x)满足方程：

Lhu(xi)u(xi1)u(xi)u(xi)u(xi1)2[p1p1]ihihi1i2hi1hi2ri[u(xi1)u(xi1)]qiu(xi)fiRi(u)hihi1 其中

…………（1.7）

1d2du1d3u1d2uRi(u)(hi1hi)([2(p)]i[p3]i[r2]i)O(h2)

4dxdx12dx2dx …………（1.8） 为差分算子Lh的截断误差，舍去Ri(u)，便得逼近边值问题（1.1）（1.2）的差分方程：

Lhu(xi)

u(xi1)u(xi)u(xi)u(xi1)2[p1p1]iihihi1hi1hi22ri[u(xi1)u(xi1)]qiu(xi)fihihi1i=1,2,…,N-1，

…………（1.9）

u0,uN

由方程（1.7）（1.9），截断误差Ri(u)可表示为

.

.

Ri(u)Lhu(xi)LhuiLh(u(xi)ui) …………（1.10）

当网格均匀，即hih(i1,2, Lhui,N)时差分方程（1.9）简化为

ui1ui11[pu(pp)upu]rqiuifi 1i111i1i1i2iiiih2h2222 …………（1.11）

这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。 截断误差为：

1d2du1d3u1d2uRi(u)h([2(p)]i[p3]i[r2]i)O(h2)

4dxdx12dx2dx …………（1.12） 所以截断误差按||Rh(u)||0或||Rh(u)||c的阶为O(h)。

在本题中，p(x)e，r(x)0，q(x)sinx1x，f(x)e2x1sinx1xxx1

x2x0，2

因为r=0方程（1.11）的系数对角矩阵是三对角矩阵。我们可以用消元法或迭代法求解方程

组（1.1）（1.2）

式（1.11）用方程组展开：

111pu[(pp)q]up1u0f1323111222h2hh2221pu[1(pp)q]u1puf11kkk2k1k122k1k1kk hhh22221pu[1(pp)q]u1puf1N13N1N13N2N1222NNNNhhh2222

111pu[(pp)q]ufp32311112221h2h2h22[1(pp)q]u1puf1p13N1N1N122N3N22N1NNhhh2222

.

.

写成矩阵形式为：

11(pp)qp0002311h23222hu1112p3(pp)qp001532225h2h22h2u2uN2111002p52(p3p5)qN22p3NNNhhhN2222uN1110002p3(pp)qN12NN1N3hh2221f1h2p12f2fN21fpN1h2N12

2.收敛性分析

根据（1.10）我们引进误差

eiu(xi)ui

则误差函数eh(xi)ei满足下列差分方程：

LheiRi(u)e0eN0i1,2,,N1

于是收敛性及收敛速度的估计问题，就归结到通过右端Ri(u)（截断误差）估计误差函数eh的问题。

由（1.12）我们知，有

lim||Rh(u)||0

h0从而差分方程满足相容条件。 若引进记号

(vi)vxi,xvivi1vvivv，(vi)xvi,xi1，(vi)vi1i，

xi,xhihi1hi111hi(hihi1)，h0h1，hNhN

222设pminc0

则可将（1.9）改写为

.

.

Lhui(p将差分解ui表成

i12(ui))qiuifi

xxuiuiui,iIh …………（2.1）

其中ui满足

~(ui)0,iIh,u0u0,uNuN …………（2.2）

xx

而ui满足

~LhuifiLhui,iIh,u0uN0 …………（2.3）

先估计uh，由

~~~~c0||(uh)||(fhLhuh,uh)Ih(fh,uh)Ih(Lhuh,uh)Ih …………（2.4）

x~20~~~

据差分格林公式

~~~(Lhuh,uh)Ih(ph(uh),(uh))Ih(qhuh,uh)Ih

xx再利用柯西不等式，有常数c1使

|(Lhuh,uh)Ih|[c1||(uh)||0c2||uh||0]||(uh)|| …………（2.5）

xx~~

将不等式（2.6）用于（2.5）右端，则

c1c21||(uh)||0||fh||1||(uh)||0||uh||0 …………（2.6）

xxc0c0c0~解差分方程（2.2，易得）

ui从而

uNu0(xia)u0 ba1|uNu0|

(ba)||(uh)||0x||uh||0(ba)max|ui|(ba)(|u0||uN|) iIh.

.

这样，

||uh||1(ba)(ba)1(|u0||uN|) …………（2.7）

利用范数||uh||1，从（2.7）推出

c1c21||(uh)||0||fh||1||uh||1 …………（2.8）

xc0c0因为

ba~||uh||0||(uh)||0

x2~因此

||uh||1||uh||0||(uh)||0x~ba(1)||(uh)||0 …………（2.9）

x2c1c2ba1(1)[||fh||1||uh||1]2c0c0~~~

联结（2.1）（2.7）及（2.9）即得差分解的先验估计：

||uh||1||uh||1||uh||1M1||fh||1M2(|u0||uN|)其中

~ …………（2.10）

M11ba(1),c02bac1c2)]2c0

M2(ba)(ba)1[1(1不等式（2.10）说明差分解连续依赖于右端和边值，因此差分格式（1.11）关于右端及边值

稳定.

根据定理1.1 ： 若边值问题的解u充分光滑，差分方程按||||R满足相容条件且关于右端稳定，则差分解uh按||||收敛到边值问题的解，且有和||Rh(u)||R相同的收敛阶。所以差分方程的解的收敛速度为O(h)。

.

2.

三．程序代码： clc clf clf

syms x;

a=1; %区间界点 b=2; %区间界点

p=exp(x); %这是p函数

q=sin(x)+1+x; %这是q函数

f=-exp(x)*(2*x+1)+(sin(x)+1+x)*x*(x-1);%这是f函数 r=0; %这是r函数.

N=10; %将区间划分的等分,这里控制！ h=(b-a)/N; %这里确定步长 value_of_f=zeros(N-1,1);%这是f

diag_0=zeros(N-1,1);%确定A的对角元

diag_1=zeros(N-2,1);%确定A的偏离对角的上对角元 diag_2=zeros(N-2,1);%确定A的偏离对角的下对角元 X=a:h:b;

u_a=0; %边界条件 u_b=2; %边界条件 for j=2:N

diag_0(j-1)=((subs(p,{x},{(X(j+1)+X(j))/2}))+(subs(p,{x},{(X(j-1)+X(j))/2})))/(h^2)+(subs(q,{x},{X(j)}));

end %获取对角元素 for j=3:N

diag_2(j-2)=-((subs(p,{x},{(X(j-1)+X(j))/2})))/(h^2)-subs(r,{x},{X(j)})/(2*h); end %获取A的第三条对角 for j=2:N-1

diag_1(j-1)=-((subs(p,{x},{(X(j+1)+X(j))/2})))/(h^2)+subs(r,{x},{X(j)})/(2*h); end %获取A的第二条对角 for j=2:N;

value_of_f(j-1)=subs(f,{x},{X(j)}); end %获取F值

value_of_f(1)=value_of_f(1)+u_a*(subs(p,{x},{(X(2)+X(1))/2}))/(h^2);

value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,{x},{(X(N)+X(N+1))/2}))/(h^2); A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%组装系数矩阵 format long

U=inv(A)*value_of_f %差分解 %fprintf('%11.5f',U) fprintf('\\n'); dx=X(2:N);

precise_value=dx.*(dx-1) %精确解 %fprintf('%11.5f',precise_value)

.

.

deta=U-precise_value' ; %误差 deta_max=max(abs(deta));%最大误差 fprintf('最大的误差是%f\\n',deta_max)

plot(X(2:N),U,'b*',X(2:N),precise_value,'r--') %差分解与精确解对比表 figure();

plot(X(2:N),deta) %误差图 结果： X的值 步长h 0.1 0.05 0.025 0.0125 精确解 2.1 0.11015 0.11003 0.11000 0.11000 0.11000 2.2 0.24026 0.24006 0.24001 0.24000 0.24000 2.3 0.39033 0.39008 0.39002 0.39000 0.39000 2.4 0.56037 0.56009 0.56002 0.56001 0.56000 2.5 0.75038 0.75009 0.75002 0.75001 0.75000 2.6 0.96038 0.96008 0.96002 0.96000 0.96000 2.7 1.19031 1.19007 1.19002 1.19000 1.19000 2.8 1.44022 1.44005 1.44001 1.44000 1.44000 2.9 1.71012 1.71003 1.71001 1.71000 1.71000 0.000380 0.000095 0.000024 0.000006 0 最大误差

以下仅给出步长为N=20,h=0.05的精确值和差分值图与误差图，其他的只要修改程序中的步长（N）值，即可同理得出。

N=20,h=0.05(图中蓝点和紫线分别是精确值图和差分图)

.

.

四．结论：

由本题可以总结出，通过用中心差分解两点边值所得的数值解能够较好地逼近方程的精确解，且区域剖分得越细，即步长越小，数值解与精确解的误差就越小，数值解越逼近精确解。

.