利用中心差分格式数值求解导数

利用中心差分格式数值求解导数

目录

一、问题描述 (2)

二、格式离散 (2)

二阶导数中心差格式离散 (2)

追赶法求解线性方程组简述 (3)

计算流程图 (5)

三、程序中主要符号和数组意义 (5)

四、计算结果与讨论 (6)

五、源程序 (9)

一、问题描述

利用中心差分格式近似导数2

2/dx y d ，数值求解 ()x dx y d 2sin 22= ()10≤≤x

1

/,0/10====x x y y 步长分别取 0001.0,001.0,01.0,

05.0=∆x

二、格式离散

将x 轴上[0,1]之间的线段按上述步长，等步长的离散为n 个小段，包括端点，共n+1个网格节点，示意图如下：

线段上边的数字表示x 轴上的坐标值，线段下边的数字表示节点编号，从0到n 编号。

二阶导数中心差格式离散

211222)2sin(x y y y dx y d x i i i ∆+-==+- 整理为线性方程形式

)2sin(2211x x y y y i i i ∆=+-+-

其中，x ∆ 为空间离散步长；i=1,2，……，n-1

包括边界条件的线性方程组如下：

边界条件

边界条件0

)

*)1(*2sin(2.........

..........)

**2sin(2..................)

*1*2sin(20

21221122100=∆-∆=+-∆∆=+-∆∆=+-=--+-n n n n i i i y x n x y y y x i x y y y x x y y y y 改写成矩阵形式：

f Ay = 其中，⎪⎪⎪⎪⎪⎭

⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧----=1012112112112101 A ，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-n n i y y y y y y 110 ，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-n n i f f f f f f 110 系数矩阵A 中仅三对角线上的数值不全为0，其余位置上的数值全为0，是

典型的对角占优的三对角矩阵，列向量f 中，)2sin(2x i x f i ∆∆=，且10==n f f ，

作为边界条件。

追赶法求解线性方程组简述

⎪⎪⎪⎪⎭

⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧----=---n n n n n i i i b a c b a c b a c b a c b A 1111110

01012112112112101

对A 做Crout 分解，即

LU A =,⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩

⎪⎪⎪⎪⎨⎧=--n n n n s r s r s r s r s L 112211

0 ，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-111111210n t t t t U 其中i i i t r s ,,为待定常数，由矩阵乘法可以得到下面的式子：

1,,3,2,1,/,,3,2,1,

,,,3,2,1,10

0000-===-=====-n i s c t n i t a b s b c t b s n

i a r i i i i i i i i i

将对角占优三对角矩阵线性方程组f Ay =等价为如下两个方程组

f L

g =，g Uy =

求解对角占优三对角矩阵线性方程组的追赶法步骤：

①输入数据i i i c b a ,,

②计算i i t s ,

③求解方程组f Lg =

n i s g a f g b f g i i i i i ,,3,2,1,

/)(/10

00 =-==-

④求解方程组g Uy = 0

,1,,2,1,1 --=-==+n n i y t g y g y i i i i n

n

⑤输出T n y y y y ),,,(10 =

计算流程图

三、程序中主要符号和数组意义

符号或数组意义

A、B、C、D、filename 用于自动更改dat文件名的字符串变量

h 离散步长

n 离散网格数，共n+1个网格节点

p 辅助变量，暂时记录网格节点上的y值

数组x，y 离散节点的x，y坐标

子程序数组a,b,c 记录系数矩阵占优对角线上的值

子程序数组f 记录线性方程组常数向量

子程序数组s,r,t,g 追赶法求解线性方程组需要用到的中间辅助向量

四、计算结果与讨论

不同步长的数值结果函数曲线与精确解的对比

从对比图中可以看到，在所取的四个计算步长下，数值计算结果与精确解都符合得相当好

Step Accuracy（Max-error）

0.05 8.152404966677018E-005

0.01 3.264604683472783E-006

0.001 3.817221055912867E-008

0.0001 9.862256566961491E-009

不同网格步长的计算精度由相应步长下所有网格节点上数值解与精确解的误差的最大值来度量，即上表中的Max-error来度量。从上表中可以看出，随着网格节点的加密，Max-error的数量级在降低，即数值解的精度提高。上述四个步长中，将线段离散成一万个网格时，数值结果的精度最高。