中心差分格式

中心差分格式

1、考虑问题

考虑二阶常微分方程边值问题：

d2uLu??2?qu?f （1） dxu(a)??,u(b)??

其中q，f为[a,b]上的连续函数，?,?为常数。

2、网格剖分与差分格式

将区间[a,b]分成N等分,分点为

xi?a?ih,i?0,1,???,N,

h=(b-a)/N,于是我们得到区间I=[a,b]的网格剖分，xi为网格节点，h为步长。

差分格式为：

ui?1?2ui?ui?1?qiui?fi2hi?1,2,???,N?1, u0??,uN??.Lhui??

3、截断误差

将方程（1）在节点离散化，由泰勒公式展开得

u(xi?1)?2u(xi)?u(xi?1)?d2u(x)?h2?d4u(x)?3????(h) ??22?4?h?dx?i12?dx?i所以截断误差为

h2?d4u(x)?Ri(u)?????(h3) 4?12?dx?i4、数值例子

u(x)?exq(x)?1?sinx

其中x??0,1?

5、求解

d2uLu??2?qu?f由，且已知 dxu(x)?exq(x)?1?sinx 可得

f(x)?exsinx

将向量式的差分格式用矩阵形式表示出来，得到矩阵形式为

?2?q1h2???1?????12?q2h2??????????12?qN?1h2????1?12?q2h2???1?u1??h2f1????u???2hf22????=????? ?????2uhf???N?1??N?1????2?q1h2???1系数矩阵A=?????????，我们求出矩阵A极其逆??2??12?qN?1h??便可求得u（x）的数值解。

6、参考文献

《偏微分方程数值解法》李荣华 高等教育出版社 《科学计算中的有限差分法》

《MATLAB程序设计教程》刘卫国 中国水利水电出版社

附件

1、程序流程图

开始 计算各初 始条件

结束 输出u的数值解 计算数值解 u=D*B 构造格式 右端矩阵B 计算矩阵d的逆D 构造系数矩阵d 2 、程序代码 %zah.m

%u(x)=exp(x), q(x)=1+sinx, f(x)=exp(x)*sinx 0

x(i)=a+i*h;%网格节点

q(i)=1+sin(x(i));%连续函数q网格节点上各值 f(i)=exp(x(i))*sin(x(i));%连续函数f网格节点上各值 end

c1=linspace(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=inv(d); %系数矩阵的逆

b1=linspace(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 input a:0 input b:1 input N:10 d =

2.0110 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0119 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0127 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0136 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0144 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0153 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0161 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0170 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0178

B =

1.0011 0.0034 0.0056 0.0079 0.0102 0.0125 0.0147 0.0170 2.7375 u =

1.1117 1.2345 1.3685 1.5143 1.6727 1.8451 2.0331 2.2393 2.4664

本文来源：https://www.bwwdw.com/article/9f9p.html