中心差两边值问题.docx

概念定义及理论基础 在数值计算中，常用差分近似微分 最简单的差分格式有向前、向后和中心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 - 1 i于是得到区间I的一个网格剖分.记h = x - x ,称h = maxh为最大网格步长i i i - 1 ii用1表示网格内点x ,x…,x 的集合，1表示内点和界点x = a,x = b的集合 h 1 2 N - 1 h 0 N取相邻节点x ,x的中点x = (x + x )(i = 1,2,…,N),称为半整数点，则由i- 1 i .丄 2 i- 1 ii-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 i i一1 = p + _i-dxJ 1 24h• 1i-h+h由(2.5)减(2.4)并除以i 2 i+i,得22h + hi i+1_ 2 h + hi i+1u( x ) — u ( x ) p( X ) i+1 1i+12(「du■ pdXt+F-hh - h+ ―i+1 i12d 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 )i i —1 id3uP~dX+ O(h3).ii -2d 3Udx3 Jih2 r d 3U胡 p I + O(h3)24 L dx3 J 丄i+2+ O(h3) (2.4)d^l + O(h3) (2.5)dx3 Jiu(x ) -u(x )—p( X ) i i-.1i—2d3u li+O(h2)h —h+ ―i+1 i4i'du'v dx 丿h —h+ ―i+1 i12id 3udx3+O(h2).i(2.6)r(x ),q = q(x ),f = f (x ),则由(2.3), (2.6)知，边值问题的解i i i i ih + hi i +1p( x )• 1i +2u(x )—u(x )i+1 1hi+1u(x ) —u(x ) —p( x ) i + i—1 h2i[u(x ) — u(x )]+ qu(x )i+1 i—1 i 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 + hi i+iu — up -^+1 ii+1 h2 i+1u —u—P —i i—1i —1 h2iri h + hi i+1[u—uquiiii+1i—1(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)i i+ 1就得出解u(x)在x的近似值uii由方程(2.7)，(2.9)，截断误差R (u)可表示为iR(u)= Lu(x)- Lu = L (u(x)- u). (2.11)i h i h i h i i以R (u)表示由R (u)定义的网函数，则由(2.8)可知截断误差按||・||或||・||的阶都是0(h)h i c 0当网格均匀，即件=h(i = 1,2,…,N)时，叫""或|Rh""I的阶提高为O(h2)此Z c 0时差分方程(2.9)简化为Luhi1h2.1 i+1.1.1i.1 i —1i+i +i—i—2222p u — (p + p )u + p u(2.12)u — ur —i+1 — + q u = f .i 2h i 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。