中心差两边值问题

概念定义及理论基础 在数值计算中，常用差分近似微分。 最简单的差分格式有向前、向后和中心3种向前差分：f'(n)=f(n+1)-f(n)向后差分： f'(n)=f(n)-f(n-1)中心差分： f'(n)=f(n+1)-f(n-1)/2考虑一般的两点边值问题：r d ( duLu = - I p <dx v dx 丿du+ r + qu = dx(2.1)(2.2) p e Ci a,b, p(x) > p > 0, r, q, f e C a,b,a, P min我们采用直接差分法来构造逼近差分格式，步骤如下： 进行网格剖分取 N +1 N+1 个节点a二x < x < < x < < x =b,0 1 i N将区间I=a,b分成N个小区间:I : x #x x ,i = 1,2,，N.i i - 1i于是得到区间I的一个网格剖分.记h = x - x ,称h = maxh为最大网格步长。i i i - 1ii用1表示网格内点x ,x,x 的集合，1表示内点和界点x = a,x = b的集合。 h1 2N - 1h0N取相邻节点x,x的中点x =(x+ x )(i = 1,2,N),称为半整数点，则由i- 1i.丄 2i- 1ii-2x < < x < < x< x =b节点:a=x < x <0 12又做成a, b的一个网络剖分，称为对偶剖分 用差商代替微商方程进行离散化 对于充分光滑的 u ，由 Taylor 展式有u(x ) -u(x )d 2udx2+ O(h 2),(2.3)p(x )u (x ) 一 u (x ) I du l h 2 ii一1 = p+ _i-dxJ 124h 1i-h+h由(2.5)减(2.4)并除以i 2 i+i,得22h + hii+1_ 2 h + hii+1u( x ) u ( x ) p( X )i+11i+12(du pdXt+F-hh - h+ i+1i12d 2dx2d ( du' I P dx v dx 丿令 p = p(x), r =.1.1 ii-i-22u(x) 满足方程：L u( x ) _ hi其中R (u ) _ (h h )ii 1 id3uPdX+ O(h3).ii -2d 3Udx3 Jih2 r d 3U胡 p I + O(h3)24 L dx3 J 丄i+2+ O(h3)(2.4)dl + O(h3)(2.5)dx3 Jiu(x ) -u(x )p( X ) ii-.1i2d3u li+O(h2)h h+ i+1i4i'du'v dx 丿h h+ i+1i12id 3udx3+O(h2).i(2.6)r(x ),q = q(x ),f = f (x ),则由(2.3), (2.6)知，边值问题的解i ii iih + hii +1p( x) 1i +2u(x )u(x )i+11hi+1u(x ) u(x ) p( x ) i + i1h2iu(x ) u(x )+ qu(x )i+1i1i ii i+1_f R(u),iid 2 f du 'pdx2 v dx 丿1+12id 3udx3d 2 u rdx2+ O(h 2).丿i(2.7)(2.8)为差分算子Lh的截断误差。舍去R.(u)，便得逼近边值问题2D, (2.2)的差分方程:Lu =hih + hii+iu up-+1ii+1 h2i+1u uPii1i 1h2iri h + hii+1uuquiiii+1i1(2.9)u = a, u = B.0N(2.10)方程(2.9), (2.10)是N- 1阶的线性代数方程组。若节点次序由左到右排列，则系数矩阵A是三对角矩阵。由于r./0，矩阵A不对称。当r ° 0即(2.1)对称时，若网格布均匀，矩 阵A亦不对称,但可以对称化，这只要在(2.9)同时乘以(h + h )即可，求解(2.9)，(2.10)ii+ 1就得出解u(x)在x的近似值u。ii由方程(2.7)，(2.9)，截断误差R (u)可表示为iR(u)= Lu(x)- Lu = L (u(x)- u).(2.11)ih ih i h i i以R (u)表示由R (u)定义的网函数，则由(2.8)可知截断误差按|或|的阶都是0(h)。hic0当网格均匀，即件=h(i = 1,2,N)时，叫""或|Rh""I的阶提高为O(h2)。此Zc0时差分方程(2.9)简化为Luhi1h2.1i+1.1.1i.1 i 1i+i +ii2222p u (p + p)u + p u(2.12)u ur i+1 + q u = f .i 2hi i i二、实践题目用中心差分格式计算如下两点边值问题dx v dxexdu | + x2u (x) = ex (x2 2ex ),2 < x < 3u(2) = e2,u(3)= e3已知其精确解为：u(x) = ex三程序编写:程序代码：运行环境MATALB 7.11.0(R2010b)clcsyms x;a=2;%区间界点b=3;%区间界点p=exp(x);%这是 p 函数q=xA2;%这是q函数f=exp(x)*(xA2-2*exp(x);%这是 f 函数 r=0;%这是 r 函数.N=1000;%将区间划分的等分,这里控制步长！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=exp(2);%边界条件u_b=exp(3);%边界条件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)+(subs(q,x,X(j); end%获取对角元素for j=3:Ndiag_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-1diag_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*(subs(p,x,(X(2)+X(1)y2)/(hA2); value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,x,(X(N)+X(N+1)y2)/(hA2); A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);% 组 装系数矩阵 format longU=inv(A)*value_of_f;%差分解precise_value=exp(X(2:N)'%精确解deta=U-exp(X(2:N)'%误差deta_max=max(U-exp(X(2:N)')%最大误差plot(X(2:N),U,'y*',X(2:N),precise_value,'r-') %差分解与精确解对比表 figure();plot(X(2:N),deta)%误差图 ：结果、X的值步长h2.12.22.32.42.52.62.72.82.9最大误差0.18.168419.028979.9793511.0290712.1886413.4696614.8849416.4486518.176420.00614670.058.166739.026009.9754811.0246512.1840313.4652214.8810416.4456518.174710.00153860.0258.166319.025269.9745111.0235512.182813.4641114.8800616.4449018.174290.00038470.01258.166209.025089.9742611.0232712.1825913.4638314.8798116.4447118.174180.0000962精确解8.166179.025019.9741811.0231812.1824913.4637414.8797316.4446518.174150以下给出以上各步长的精确值和差分值图与误差图：N=10,h=0.1N=20,h=0.05N=40,h=0.025误差N=80,h=0.025精 确 值 与 差 分 值 N=80， h=0.025四：总结 误差在数值计算中是不可避免的，误差的传播和积累直接影响到计算结果的 精度。在研究算法的同时，必须注重误差分析，使建立起来的算法科学有效。 通常用函数的泰勒展开对误差进行估计 。从这次的课程设计中我学到，在一定的误差范围内，我们完全可以用数值解 来代替精确解（解析解），因为大多时候精确解的求解过程是复杂，并且繁琐的。 很少有复杂的方程得到其精确解，而从以上的中心差分看，误差是如此的小，即 便分割的不是很细。