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

利用中心差分格式数值求解导数 目录 TOC \o "1-5" \h \z 一、问题描述 2 HYPERLINK \l ”bookmark14” \o "Current Document” 二、格式离散 2 \o "Current Document” 二阶导数中心差格式离散 2 HYPERLINK \l ”bookmark34” \o "Current Document" 追赶法求解线性方程组简述 3 \o "Current Document” 计算流程图 5 Current Document" 三、程序中主要符号和数组意义 5 \o "Current Document" 四、计算结果与讨论 6 \o ”Current Document” 五、源程序 9 问题描述 利用中心差分格式近似导数 d2y∕dχ2 ，数值求解 d2y d2ydx2 =Sin 2x 0<x<1 y∕χ=o = 0， y/X^= 1 步长分别取 .05, 0.01, 0。001， 0。0001 】、格式离散 将X轴上［0,1］之间的线段按上述步长，等步长的离散为 n个小段，包括端点, 共n+1个网格节点，示意图如下: n-2 n-1 n n-2 n-1 n 线段上边的数字表示 X轴上的坐标值，线段下边的数字表示节点编号，从 0到n 编号。 二阶导数中心差格式离散 Sin(2x)二yi 1dx Sin(2x)二 yi 1 dx x2 整理为线性方程形式 2 W= —2yi f x sin(2x） 其中， X为空间离散步长；i=1，2 ,……，n-1 包括边界条件的线性方程组如下： y° = 0 边界条件 2 y° —2yι y2 — X Sin(2＊1＊ :x） yij - 2yi yi 1 = x2 sin（2* i * :x） 2 y∏^ -2y∏j y^ X sin(2＊ （n—1）* ：x) y∏ =0 边界条件 改写成矩阵形式： Ay = f 其中,A = * 1 1 0 ’ -2 1 1 —2 1 + + + + + + ?，y = * y。 * ≡ yi ≡ ,，f 'f°、 f1 ■- f ? 1 i ≡ 1 —2 1 ■■ ■ y∏-1 f ∏ _1 1 —2 1 、y∏ 、、f∏ 「 0 1“ 系数矩阵A中仅三对角线上的数值不全为 0，其余位置上的数值全为 0，是 典型的对角占优的三对角矩阵， 列向量f中，fi = x2 Si∏(2i?〉x),且f0 = f∏ =1， 作为边界条件. 追赶法求解线性方程组简述 彳0 1 -2 1 1 —2 1 ?0 C0aIbIGIA1 -2 11 -2 10 1a∏4 b ?0 C0 aIbIGI A 1 -2 1 1 -2 1 0 1 a∏4 b∏ -A a∏ bn 对A做CroUt分解，即 r1 to 、 「1 S1 1 t1 「2 S2 1 t2 A = LU ,L + + ＞,U = ? 「n_d Sn_J 1 tn」 rn Sn 1 其中S ， Γi， I为待定常数,由矩阵乘法可以得到下面的式子: ai， i = 1,2，3， ,n SO= b0 ， tO =汁 bo Si= bi - ajtj_i， i = 1，2,3, ，n I = & / Si, i = 1，2，3, ，n -1 将对角占优三对角矩阵线性方程组 Ay = f等价为如下两个方程组 Lg = f，UyHg 求解对角占优三对角矩阵线性方程组的追赶法步骤： 输入数据ai， b , Ci 计算Si , ti 求解方程组Lg= f g。= f°∕b° gi = (fi —aigi^）∕Si， i =1，2,3， ，n 求解方程组Uy = g yn =gn y^ gi —ti? 1， i = n -1， n —2, ，1，0 输出y =(y°, ％,…y)τ 计算流程图 开始 对赋值 追赶法求解线性方程 组，得到X f 输出结果，后处理 ! E 结束 三、程序中主要符号和数组意义 符号或数组 意义 A B、C、D file name 「用于自动更改dat文件名的子符串变量 h 离散步长 n :离散网格数，共n+1个网格节点 P 辅助变量，暂时记录网格节点上的 y值 数组X, y 离散节点的X，y坐标 子程序数组a，b,c ：记录系数矩阵占优对角线上的值 子程序数组f 「记录线性方程组常数向量 子程序数组s,r,t,g 追赶法求解线性方程组需要用到的中间辅助向量 四、计算结果与讨论 不同步长的数值结果函数曲线与精确解的对比 X X SteP=O。01 0。8 numerical exact 0。6 0。4 0.2 i 1 i 1 I 0 0.2 0。4 0.6 0。8 1 X