中心差分解两点边值问题

1、中心差分解两点边值问题一.题目用中心差分格式计算如下两点边值问题ddxex j + (sinx + l + x)w(x) = -ex (2x +1) +(sin x + l + x)x(x-1),1 <x<2 k dx)"=0/(2) = 2已知其精确解为二理论作为模型，考虑两点边值问题:汕莎w矿能）,(1.1)u(a) = a.u(h) = "" (1.2)假定 P e C也,b, p(x) > pn > 0,厂,g, / w Ca,b, a, 0 是给定的常数。1.建立差分格式（1） 区域网格剖分首先取"+1个节点:a = x

2、Q<x vV兀 <<& =b,将区间/ =小分成N个小区间：/.:齐.I <x<xpz = l,2,-AT于是得到区间/的一个网格剖分。记渥/z = max说为网格最大步长。用人表不网格内点旺,，畑的集合，匚表示内点和界点 x = a,xN =b 的集合。取相邻节点兀一,x,的中点x丄=(x,._! +x,.)(z = l,2,L ,N),称为半整数点。则由节点a = x0<x <xy <L <x )<L x <x =b55二*v*又构成Mb的一个网格剖分,称为对偶剖分。(2) 微分方程的离散，建立相应差分格式用差商代替

3、微商，将方程(1.1)在内点兀离散化. 注意对充分光滑的“，由Taylor展式有害严嚕+呼愣ZS)(1.3)心气週嗨. +0(巧r&dx 二 24 ax 二(1.4)=-1严尤"仝a+o闪) dxQ 24 7 dxy z/心2"(几)一"(兀)哙屮钞卩鈴+则)(1.5)由(1.5)减(1.4),并除以葺5得-px j"g)W心心)gj2“ digr du 、 hM -/?, di.八八 2、=(p=L -=)+ 丄1卩曰+0(力)勺 +1r d z du J 咤F+ 4检(1.6)du r du ndxJ "12 dx3/r+l -lr

4、 , d2 , du、 h. -/. r亠八、丄(p )1 +-jrl +。"厂)dx dx 12 dxdx令 P I = p(x , ),< = r(xiqt = (x.),/ = f(x, 由(1.3)( 1.6)知, ii边值问题的解如满足方程:点gjg丿+gr+恥)(1.7)其中d -R3叽兀怙(撤+9 弘-新等山曲)(1.8)为差分算子®的截断误差，舍去弘)，便得逼近边值问题(1.1) (1.2)的差分方程:LT宀p九广心5 |牛皿+ 人 + 叱-7 Igj ht； 【”(和)一) + qg) = fhi + h®(1.9)i=l，2,，N-1,%

5、=Sn =卩由方程(1.7) (1.9),截断误差弘)可表示为K (“) = S(xJ 一 Lhui = Lh (心)一的)(1.10)当网格均匀，即 h, =h（i = ,2,.,N） 时差分方程（1.9）简 化为m 扣屮厂（务“严+匕导屮仝尹讪 rJ厶(1.11)这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。截断误差为:R3 = h（£咲（半几 +【"一扣兽，）+%）4 dx clx 12 ax' 2 dx(1.12)所以截断误差按II恥）II。或II恥）II,的阶为。（心O在本题中f = exr(x

6、) = O f g(x) = sinx + l + x ,/(x) = -ev(2j+l)+(sinx+l+x)x(x-l)a=090=2因为r=0方程（1.11）的系数对角矩阵是三对角 矩阵。我们可以用消元法或迭代法求解方程组（1.1） （1.2）式（111）用方程组展开：-2 Pf '2 + 吉（卩扌 + P）+ 9 J也一 JT P1«O = /1«M1r1z、1f-V仏严+中3碍+你帆-F罠二fk-和M+占（叮+P.J +汕%-h.如2认N< 片'N=N< 一加3"2+£（必+从）+讪二久+加“/广7/厂于于” 7=&

7、gt;<h- N飞一吒 丄（化 1 + /<_ 3 ） + V-l kv-1 - /<. 3«.V-2 = Zv-1 + p- l,丄0写成矩阵形式为:也 000和作孚)+色洛作 00 0-右g瓠家P护MfzM厶 00和务+%)+如p（p3+pl）+t?lJ?说00fl2 收敛性分析根据(110)我们引进误差ei=u(xi)-ui则误差函数曲r满足下列差分方程:厶局=尺(“)®=s=°于是收敛性及收敛速度的估计问题，就归结到通过右端a (“)(截断误差)估计误差函数S的问题。由(1.12)我们知，有limll/?/r(W)ll=O 丿D从而差分方

8、程满足相容条件。若引进记号1 1 1 hi = -+&利)9 ho =-/?! , Hn = -hN则可将(1.9)改写为Lm = -5 (",)+9" =/(2.1)将差分解©表成的=Ui + UiJ e Ih其中匚满足()“ =OJeZ/p wo = uun = uNXX(2.2)lUi = / -£y,Mij=un =0(2.3)而债满足先估计Uh 9由z.一.Co II (讥II詁(齐一厶，仙川 =(fUh)lb 一(厶動M九V(2.4)据差分格林公式(Uh 上 h)h =(Ph(Uh)3)h + 冋 hVX再利用柯西不等式，有常数q使l

9、(Ift UhyUh)lfi l< q II (血).ll0 +c2 II Uh ll0ll (w/r)- II(2.5)将不等式(2.6)用于(2.5)右端，则1rC "II (渤) 110< II All | +4II (Uh). 110 +二 II Uh ll0x Coc°x c°(2.6) 解差分方程(2.2,易得)弘=一(兀_a) +如b_a从而11(劲).|。=|你_心|x J(b-d)II Uh ll0< J(b一g) max I Hi l< J(b_a)(l m0 I +1 uN I)re/*这样，II Uh ll5 J(一

10、 d) + (/?-G)“(I H() I + I UN I)(2.7)利用范数讪，从(2.7)推出II (W/,)- 110< 丄 II 厶 II !r %c°(2.8)因为llwd ll(w/r)_ ll02 *因此II uh IIj<IIm/i ll0 +11(附).ll()X"+¥)ii (h)-%(2.9)2-vs(i + ¥)丄 ii 齐|“+匚£1| 心 iij2联结(2.1) (2.7)及(2.9)即得差分解的先验估计:II Uh IIVIm/t 11| + II Uh 11|HA11-! +M2(|MOI + |W

11、A'1)(2.10)其中心(1+字，c° 2M, = J(b-a) + (b-aYl + (1 +-"2 c“不等式（2.10）说明差分解连续依赖于右端和边值，因此差分格式（1口）关于右端及边值稳定.根据定理1.1 : 若边值问题的解u充分光滑, 差分方程按叫满足相容条件且关于右端稳定，则 差分解你按*收敛到边值问题的解，且有和 | 恥）相同的收敛阶。所以差分方程的解的收敛速度为0(冋O二程序代码:syms x;a=l;b=2; p=exp(x);%区间界点%区间界点%这是P函数q=sin(x)+l+x;%这是 q 函数f=-exp(x)*(2*x+l)+(sin(

12、x)+l+x)*x*(x-l);% 这是 f函数 r=0;N=10;%这是r函数.%将区间划分的等分,这里控制！ h=(b-a)/N;%这里确定步长value_of_f=zeros(N-l,l)； % 这是 f diag_0=zeros(N-l9l);%确定 A 的对角元 diag_l=zeros(N-24)；%确定A的偏离对角的上 对角元diag_2=zeros(N-24)；%确定A的偏离对角的下 对角元X=a:h:b;u_a=O; %边界条件 u_b=2; %边界条件forj=2:N diag_0(j.l)=(subs(p,xJ(X(j+l)+X(j)/2)+(su bs(pJxJ(X(j

13、-l)+X(j)/2)/(hA2)+(subs(qJxJ x(j)；end%获取对角元素for j=3:Ndiag_20-2)=-(subs(p,x,(X(j-l)+X(j)/2)/(hA 2)-subs(r,x,X(j)/(2*h);end%获取A的第三条对角forj=2:N-ldiag_10-l)=-(subs(p,x,(X(j+l)+X(j)/2)/(h 人 2)+subs(r,x,X(j)/(2*h);end%获取A的第二条对角for j=2:N;value_of_f(j-l)=subs(fJxJX(j);end%获取F值value_of_f(l)=value_of_f(l)+u_a(

14、subs(pjxj( X +X(l)/2)/(22);value_of_f(N-l)=value_of_f(N-l)+u_b(subs(pJ x,(X(N)+X(N+l)/2)/(22);A=diag(diag_0)+diag(diag_l,l)+diag(diag_2rl); %组装系数矩阵%差分解%精确解%误差format long U=inv(A)*value_of_f %fprintf(,%11.5fU) fprintfC1);dx=X(2:N);precise_value=dx. * (dx 1) %fprintf(f % 11.5fprecise_value) deta=U-pre

15、cise_value1 ;deta_max=max(abs(deta); % 最大误差fprintf(f 最大的误差是 %fndeta_max)plot(X(2:N)，LVb*'，X(2:N)，precise_valii£，t)%差分解与精确解对比表%误差figure(); plot(X(2:N),deta) 图结果:值2.12.22.32.42.52.62.72.82.9最大误差0.10.110150.240260.390330.560370.750380.960381.190311.440221.710120.0003800.050.110030.240060.390080.560090.750090.960081.190071.440051.710030.00()0950.0250.110000.240010.390020.560020.750020.960021.190021.440011.710010.0000240.01250.110000.240000.390000.560010.750010.960001.190001.440001.710000.000006精确解0