中心差分格式.docVIP

中心差分格式 1、考虑问题 考虑二阶常微分方程边值问题: Lu =d2u Lu = d2u dx2 qu = (1) u(a) = :,u(b)=: 其中q, f为［a,b］上的连续函数，:-为常数 2、网格剖分与差分格式 将区间［a,b］分成N等分,分点为 xi 二 a ih,i 二 0,1, , N , h=(b-a)/N,于是我们得到区间l=［a,b］的网格剖分，x,为网格节点，h为步长 差分格式为： ..._ -25 7_1 一 丄 Lhui _ ~2 qiui _ fi h i =1,2, ,N -1, Uo =：从= 3、截断误差 将方程(1)在节点离散化，由泰勒公式展开得 u(x J -2u(x) u(x Jh2 u(x J -2u(x) u(x J h2 [d2u(x)1 +h2 ~d4u(x) I J dx2」12 ] dx4 一 O(h3) 所以截断误差为 R(u)二h2 R(u)二 h2 _d4u(x)l 77 ,1 . 4 ~ 12 |L dx O(h3) 4、数值例子u(x) = exq(x) =1 sin x 4、数值例子 u(x) = ex q(x) =1 sin x 5、求解 其中X 0,1 1 Lu d2u 捺3，且已知 u(x)二ex q(x) =1 sin x 可得 f (x)二 ex sin x 将向量式的差分格式用矩阵形式表示出来，得到矩阵形式为 2 qih 2 qih2 -1 —1 … 2 q2h2 -1 + + -1 … u2 *? : . ? 2 2 qNjh _un_1 -h2f, 1 | h2f2 ]h2fN」+ P 一 2 q-j 2 q-jh2 _1 系数矩阵A= -1 2 q2h2 -1 + + -1 +. ，我们求出矩阵A极其 2 2 + qN^h 一 逆便可求得u (x)的数值解。 6、参考文献 《偏微分方程数值解法》李荣华高等教育出版社 《科学计算中的有限差分法》 《MATLAB程序设计教程》刘卫国 中国水利水电出版社 附件 1、程序流程图 计算各初 始条件 构造系数矩阵d 1 r 计算矩阵d的逆D 1 F 构造格式 右端矩阵B F 计算数值解 u=D*B 输出u的数值解 2、程序代码 %zah.m %u(x)=exp(x), q(x)=1+s inx, f(x)=exp(x)*s inx 0x1 a=input(input a:);% 输入区间条件 a b=input(input b:);% 输入区间条件 b N=i nput( in put N:);% 输入等分数 N h=(b-a)/N;% 步长 for i=1:N-1 x(i)=a+i*h;%网格节点 q(i)=1+sin(x(i));%连续函数q网格节点上各值 f(i)=exp(x(i))*sin(x(i));%连续函数f网格节点上各值 end c1=li nspace(2+q(1)*h*h,2+q(N-1)*h*h,N-1); % 系数矩阵主对角线上各值 c2=linspace(-1,-1,N-2); % 系数矩阵次对角线上各值 d1=diag(c1,0); d2=diag(c2,-1); d3=diag(c2,1); d=d1+d2+d3 %构造系数矩阵 D=i nv(d); %系数矩阵的逆 b1=li nspace(h*h*f(1),h*h*f(N-1),N-1); for i=1:N-1 if i==1 b2(i)=exp(a); else if i==N-1 b2(i)=exp(b); else b2(i)=0; end end end B1=b1+b2;%构造差分格式等式右边的向量 B=B1 %差分格式等式右边向量的转置 u=D*B %输出u (x)数值解 输出结果: clear zah in put a:0 in put b:1 in put N:10 2.0110 -1.00000 2.0110 -1.0000 0 -1.0000 2.0119 0 0 -1.0000 0 0 0 0 0 0 0 0 0 0 0 -1.0000 0 0 2.0127 -1.0000 0 -1.0000 2.0136 -1.0000 0 -1.0000 2.0144 0 0 -1.0000 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 0 2.0153 -1.0000 0 -1.0000 0 0 0 0 0 -1.0000 2.0161 -1.0000 0 0 0 0 0 0 -1.0000 2.0170 0 0 0 0 0 0 0 -1.0000 2.0178 1.0011 0.0034 0.0056 0.0079 0.0102 0.0125 0.0147 0.0170 2.7