中心差分解两点边值问题.docx

1、中心差分解两点边值问题一题目用中心差分格式计算如下两点边值问题d x du xe sinx 1 x u x e 2x 1 sinx 1 x x x 1 ,1 x 2dx dxu(1) 0,u(2) 2已知其精确解为u(x) x(x 1)二.理论作为模型，考虑两点边值问题：Lu -d(p(x)du)晋 q(x)u f (x),a x b ( 1.1)dx dx dxu(a) ,u(b) ( 1.2)1是给定的常数。假定 p C a,b, p(x) Pmin 0,r,q, f Ca,b,1.建立差分格式 (1) 区域网格剖分 首先取N 1个节点:将区间I a,b分成N个小区间:于是得到区间I的一个

2、网格剖分。记Ii : xi 1 x xi, i 1,2,L N.xi xi 1，称h maxh为网格最大步长。用Ih表i示网格内点x1, x2, L , xN勺集合，Ih表示内点和界点 沧 a,xN b的集合。取相邻节点xi 1, xi的中点x 1i 一 21 .(Xi 1 xj(i 1,2,L , N )，称为半整数点。则由节点2xox-i x? L2 2Xn b又构成a,b的一个网格剖分，称为对偶剖分。用差商代替微商，将方程(2) 微分方程的离散，建立相应差分格式1.1)在内点Xi离散化注意对充分光滑的u，由Taylor展式有u(Xi J u(Xi Jhi hi 1hi 1 hi2霁idx

3、O(h2)(1.3)P;)U(Xi).U(Xi1) P 乳 1dx i 712i O(h3)p%. dx i.2 , 3hi r d U 3.:p 3j 1 O(h )24 dx | ?(1.4)p(X1)U(Xi1)U(Xi)i 2hi 1pl 1dx i 22 3h 1 r du】 -rylP ali O(h )24 dx(1.5)hi h 1 得2 ，寸p(X1 严 J 呛)22 dunb(PdXi由(1.5)减(1.4),并除以hi 1呂(Pli hdx dxr du_ P dx ihid21)2p(xi1)U(Xi)hU(Xi1)i 2 hi3P 吗 O(h2)12 dx3 ip 岁

4、i O(h2) dxdu X h 硬(Pi1.6)令P. 1i2满足方程:p(xi1), ri r(xj,qi2q(xj, fif(X),则由(1.3) ( 1.6)知，边值问题的解u(x)LhU(xJu(X1) u(xjhi 1u(x)12u(xi 1)hi1.7)其中亡曲1)u(Xi 1)qiU(Xi)R(u)R(u) (hi 11hi)(；dx- dx1 d3u扣i) O(h2)2 dx(1.8)为差分算子Lh的截断误差，舍去 Ri(u)，便得逼近边值问题(1.1) (1.2)的差分方程:LhU(Xi) 土 p 1U(X1) U(Xi)hi hi 1 i 2hi 1U(Xi) U(Xi1

5、)1 . 2hi(1.9)u(Xi 1) u(Xi 1) qiU(Xi)hi h 1i=1,2,N-1 ,Uo ,Un由方程(1.7) (1.9),截断误差 R(u)可表示为(1.10)当网格均匀，即R(u) LhU(xJ LhUi Lh(u(xJ ujh（i 1,2,L ,N）时差分方程（1.9）简化为（1.11）(1.12)xe 2x 1 sinx 1 x x x 1所以截断误差按|Rh（u）|o或|Rh（u）|C的阶为O（h）。在本题中，p(x) ex , r(x) 0 , q(x) sinx 1 x,f(x)因为r=0方程（1.11 ）的系数对角矩阵是三对角矩阵。我们可以用消元法或迭代

6、法求解方程组（1.1）（ 1.2）式（1.11 ）用方程组展开：P1)弭21 f2 P1 U0 f1h 2ku1 -2kp1 -2pk/Vku1 -2-pkUN3 - 2NqPN1 1P3U 2(P3h 2 h 2占(5 Pi) qh 2 211 h2P30L L001h2P32(p5 p3) Qh 2 21严L L00LLL0L L01h2PN|右(Pn3 PN 5)qN 22 21h2PN|0L L001RPn；*Pn1 PN 3) Ch2 2写成矩阵形式为:uU2MUN 2UN 1fl2收敛性分析根据（1.10）我们引进误差e u(x) Ui则误差函数eh （x）e满足下列差分方程:Le

7、i 驾 i 1,2,L,Neo eN 0的问题。由（1.12）我们知，有lhimllR(u)ll 0从而差分方程满足相容条件。 若引进记号设 Pmin C0则可将（1.9）改写为LhUi (p. i(uj) qq fiI XX2将差分解U|表成U| UI Ul,i I h ( 2.1 )其中UI满足(ui) 0,l Ih,U0 Uo,UN Un ( 2.2)XX而Ui满足LhUi fl Lh Ui ,i I h,Uo un 0 ( 2.3)先估计Uh，由Coll(Uh)xl2 (fh LhUh,Uh)|h (fh,Uh)|h (LhUh,Uh)|h ( 2.4)据差分格林公式(LhUh,Uh)

8、ih (Ph(Uh) ,(Uh) )ih (qhUh,Uh)ihX X再利用柯西不等式，有常数 c1使|(LhUh,Uh)|h | Cill(Uh) llo C2|Uh |0| (Uh) ll ( 2.5)X X将不等式(2.6)用于(2.5)右端，则 1 c1 C2ll(Uh) |0 -|fhll1 -ll(Uh) |0 -|Uh|0 ( 2.6)X Q Q X Q解差分方程（2.2，易得）从而Un UibU0 (Xiaa) U01ll(Uh) |0i|Un U0lX.(b a)|uh| , (b a)max|ui| 、(b a)(|u| |Un |)I I h这样,l|uh|1 ._(b

9、a) (b a)1(|Uo| |Un|)(2.7)1ll(Uh) |o -|fh|1x CoC1 色|3|1Co (2.8)因为b a .l|Uh|o_ |(uh) |o因此2 x利用范数|uh仏，从（2.7）推出l|Uh|1 |Uh|。|(Uh) |oX(1 )ll(Uh) |02 xb a 1 c1 c2(1 ) llfhlli - 2|uh|i2 Co Co(2.9)联结（2.1） （2.7）及（2.9）即得差分解的先验估计:|Uh|1 |uh|1 |uh|1M1|fh|1 M2(|u| |Un|)其中(2.10)1 “ b a、M1 -(1 ),co 2M2 ；(b a) (b a)1

10、1 (1 呼声绍2 Co不等式（2.10）说明差分解连续依赖于右端和边值，因此差分格式（ 1.11）关于右端及边值稳定.根据定理1.1 : 若边值问题的解u充分光滑，差分方程按|g|R满足相容条件且关于右端稳定，则差分解uh按|g|收敛到边值问题的解，且有和 |Rh（u）|R相同的收敛阶。所以差分方程的解的收敛速度为 O（h2）。三程序代码：clcclfclfsyms 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

11、; %这是 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:Ndiag_0(j-1)=(subs(p,x,(X(j+1)+X(j)/2)+(subs(p,x,(X(j-1)+X(j)/2)/(hA2

12、)+(subs(q,x ,X(j);end % 获取对角元素for j=3:N diag_2(j-2)=-(subs(p,x,(X(j-1)+X(j)/2)/(hA2)-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)/(hA2)+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*(s

13、ubs(p,x,(X(2)+X(1)/2)/(hA2);value_oLf(N-1)=value_of_f(N-1)+u_b*(subs(p,x,(X(N)+X(N+1)/2)/(hA2); A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);% 组装系数矩阵 format longU=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_v

14、alue ; % 误差deta_max=max(abs(deta);% 最大误差fprintf(最大的误差是 %fn,deta_max) plot(X(2:N),U,b*,X(2:N),precise_value,r-) % 差分解与精确解对比表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.56009

15、0.750090.960081.190071.440051.710030.0000950.0250.110000.240010.390020.560020.750020.960021.190021.440011.710010.0000240.01250.110000.240000.390000.560010.750010.960001.190001.440001.710000.000006精确解0.110000.240000.390000.560000.750000.960001.190001.440001.710000以下仅给出步长为 N=20,h=0.05的精确值和差分值图与误差图 ，其他的只要修改程序中的 步长（N ）值，即可同理得出。N=20,h=0.05（图中蓝点和紫线分另忧精确值图和差分图 ）21.31.6141.21D8D.60.4D21 11 12 13 1 4 15 1.6 1.7 1.6 1.S 2四结论：由本题可以总结出，通过用中心差分解两点边值所得的数值解能够较好地逼近方程的精确解，且区域剖分得越细，即步长越小， 数值解与精确解的误差就越小，数值解越逼近精确解。PN J 1UN 1