新乡市政协委员张亚玲,石家庄张亚玲,新乡张亚玲,张亚玲,数值分析实验指导书2013

实验 1.1(体会数值计算精度与步长关系的实验)实验目的:数值计算中误差是不可避免的,要求通过本实验初步认识数值分析中两个重 要概念:截断误差和舍入误差,并认真体会...

数值分析课程 实验指导书 主编：张亚玲 张卫国 西安科技大学 计算机学院 目 第一部分 第二部分 第三部分 录 数值实验报告格式 ..........................................- 2 数值实验报告范例 ..........................................- 4 内容回顾及数值实验 ......................................- 8 - 第 1 章 绪论 ...................................................................- 8 第 2 章 非线性方程求解 ............................................ - 11 第 3 章 线性方程组的解法 ........................................- 17 第 4 章 插值法 .............................................................- 20 第 5 章 函数逼近与曲线拟合 ....................................- 23 第 6 章 数值积分 .........................................................- 25 第 7 章 矩阵特征值与特征向量的计算 ....................- 31 第 8 章 常微分方程数值解 ........................................- 35 第四部分 MATLAB 入门 ................................................- 38 - 数值分析实验指导书 前 言 该实验指导书是根据《数值分析》课程编写的配套数值实验教材。

《数值分析》是西 安科技大学本科生与硕士研究生的必修课程，学习本课程的最终目的，是用计算机解决 科学和工程实际中的数值计算问题，因此熟练地在计算机上实现算法是必备的基本技能。

数值实验是数值分析课程中不可缺少的部分，利用计算机进行数值实验，以消化巩固所 学的内容，增加对算法的可靠性、收敛性、稳定性及效率的感性认识，体会和重视算法 在计算机上实验时可能出现的问题。学生通过选择算法、编写程序、分析数值结果、写 数值实验报告等环节的综合训练，逐步掌握数值实验的方法和技巧，获得各方面的数值 计算经验，培养学生运用所学算法解决实际问题和进行理论分析的能力。本书的编写过 程中首先给出数值实验的报告形式，并给出了数值实验报告范例供学生参考。其次是对 每章节内容的复习和回顾，让学生更好的掌握所学的内容，从而给出需要学生编写的数 值实验，并给出了实验编程的思路。最后介绍了 MATLAB 的入门知识。

该实验指导书由张卫国老师指导，张亚玲编写完成，在此对张卫国老师表示感谢！ -1- 数值分析实验指导书 第一部分 数值实验报告格式 一个完整的实验，应包括数据准备、理论基础、实验内容及方法，最终对实验结果 进行分析，以达到对理论知识的感性认识，进一步加深对相关算法的理解，数值实验以 实验报告形式完成，实验报告格式如下： 一、实验名称 实验者可根据报告形式需要适当写出。 二、实验目的及要求 首先要求做实验者明确，为什么要做某个实验，实验目的是什么，做完该实验应达 到什么结果，在实验过程中的注意事项，实验方法对结果的影响也可以以实验目的的形 式列出。 三、算法描述（实验原理与基础理论） 数值实验本身就是为了加深对基础理论及方法的理解而设置的， 所以要求将实验涉及 到的理论基础，算法原理详尽列出。 四、实验内容 实验内容主要包括实验的实施方案、步骤、实验数据准备、实验的算法以及可能用 到的仪器设备。 五、程序流程图 画出程序实现过程的流程图，以便更好的对程序执行的过程有清楚的认识，在程序 调试过程中更容易发现问题。 -2- 数值分析实验指导书 六、实验结果 实验结果应包括实验的原始数据、中间结果及实验的最终结果，复杂的结果可以用 表格形式列出，较为简单的结果可以与实验结果分析合并出现。 七、实验结果分析 （对算法的理解与分析，包括改进与建议） 实验结果分析是数值实验的重要环节，只有对实验结果认真分析，才能对实验目的、 实验方法进一步理解，对实验的重要性充分认识，明确数值分析的实用范围及其优缺点。 -3- 数值分析实验指导书 第二部分 数值实验报告范例 为了更好地做好数值实验并写出规范的数值实验报告，下面给出一简单范例供读者 参考。 数值实验报告 一、实验名称 误差传播与算法稳定性 二、实验目的 1．理解数值计算稳定性的概念。

2．了解数值计算方法的必要性。

3．体会数值计算的收敛性与收敛速度。

三、实验内容 1 计算 四、算法描述 In ? ? 0 xn dx x ? 10 ， n ? 1, 2, ,10 由 则 In ? ? 1 0 n ?1 1 x xn dx I n ?1 ? ? dx 0 x ? 10 x ? 10 ，知 I n ? 10I n?1 ? ? 1 0 1 x n ? 10x n?1 1 dx ? ? x n?1dx ? 0 x ? 10 n -4- 数值分析实验指导书 可得递推关系 1 ? 10 I n ?1 I ? 1． n n I n ?1 ? ， n ? 1,2,?,10 2． 1 1 ( ? In ) 10 n ， n ? 10,9,?,1 下面分别以 1，2 递推关系求解 方案 1 1 ? 10 I n ?1 In ? n ， n ? 1,2,?,10 当 n ? 0 时， I0 ? ? 1 0 11 1 ? dx ? x ? 10 ㏑ 10 ㏑ 1.1 , 递推公式为 1 ? ?I n ? ? 10I n ?1 , n ? 1,2,?,10 n ? ? ?I 0 ? ln 1.1 I n ?1 ? 1 1 ( ? In ) 10 n ， n ? 10,9,?,1 （1） 方案 2 当 0 ? x ? 1 时， 则 1 n xn 1 x ? ? xn 11 x ? 10 10 ? 即 1 0 1 1 1 n xn 1 n x dx ? ? dx ? ? x dx 0 0 10 11 x ? 10 1 1 ? In ? 11(n ? 1) 10(n ? 1) 1 1 1 21 I 10 ? [ ? ]? 2 11(10 ? 1) 10(10 ? 1) 220(10 ? 1) 取递推初值 -5- 数值分析实验指导书 1 1 ? I n ?1 ? ( ? I n ), n ? 10,9,?,1 ? 10 n ? ? 21 ?I 10 ? ? 220(10 ? 1) 递推公式为 ? 取递推公式（1）中的初值 （2） I 0 ? ln1.1 ? 0.095310 ，得 n ? 1,2,?,10 1 ? ?I n ? ? 10I n ?1 , n ? ? ?I 0 ? 0.095310 取递推公式（2）中的初值 I10 ? 0.008678 ，得 1 1 ? ?I n ?1 ? ( ? I n ), n ? 10,9,?,1 10 n ? ? ?I 10 ? 0.008678 五、程序流程图 由于实验方案明显、简单，实现步骤及流程图省略。

六、实验结果 计算结果如下：

n 0 1 2 3 4 5 6 7 ~ I n (1) 0.095310 0.046900 0.031000 0.023333 0.016667 0.033333 -0.166667 1.809524 -6- * In (2) 0.095310 0.046898 0.031018 0.023153 0.018465 0.015353 0.013138 0.011481 数值分析实验指导书 8 9 10 七、实验结果分析 0.010188 0.009232 0.008678 由递推公式（1）知当 I 0 ? ln1.1 时， I n 应当为精确解，递推公式的每一步都没有误 差的取舍， 但计算结果 I 5 ? 0.033333 ＞ 0.016667 ? I 4 ，I 6 出现负值。

由此看出， 当n 较大时，用递推公式（1）中的 I n 近似 ~ ~ ~ ~ I n 是不正确的。 ~ ~ I 主要原因是初值 0 ? 0.095310 不是精确值，设有误差 e( I 0 ) ，由递推公式（1）知 ~ ~ e(I n ) ? ?10e(I n?1 ) 则有 ~ ~ ~ ~ e(I n ) ? ?10e(I n?1 ) ? ?100e(I n?2 ) ? (?10) n e(I 0 ) n 误差 e( I n ) 随 n 的增大而迅速增加，增加到 e( I 0 ) 的 (?10) 倍。由此可见，递推公式计 ~ ~ 算的误差不仅取决于初值的误差，公式的精确性，还依赖于误差的传递即递推计算的稳 定性。 I ? 0.008678， I n 为估计值，并不精确，有 由递推公式（2）知 10 * e( I n ?1 ) ? ? e( I 10 ) ? 1 1210 ， 而由 1 * e( I n ) 10 得 * e( I 0 ) ? (? 1 n * ) e( I n ) 10 误差 * e( I 0 ) 随递推公式逐步缩小。

综上所述， 在递推计算中， 数值计算方法是非常重要的， 误差估计、误差传播及递推计算的稳定性都会直接影响递推结果。 -7- 数值分析实验指导书 第三部分 内容回顾及数值实验 第 1 章 绪论 一、主要内容回顾 1、算法：电子计算机实质上只会做加、减、乘、除等算术运算和一些逻辑运算，由 这些基本运算及运算顺序规定构成的解题步骤，称为算法．它可以用框图、算法语言、 数学语言或自然语言来描述。用计算机算法语言描述的算法称为计算机程序。

（如 c—语 言程序，c++语言程序，Matlab 语言程序等） 。

2、最有效的算法：应该运算量少，应用范围广，需用存储单元少，逻辑结构简单， 便于编写计算机程序，而且计算结果可靠。

3、 算法的稳定性：

一个算法如果输入数据有误差， 而在计算过程中舍入误差不增长， 则称此算法是数值稳定的，否则称此算法为不稳定的。换句话说：若误差传播是可控制 的,则称此算法是数值稳定的,否则称此算法为不稳定的。

4、控制误差传播的几个原则：

1）防止相近的两数相减； 2）防止大数吃小数； 3）防止接近零的数做除数； 4）要控制舍入误差的累积和传播； 5）简化计算步骤,减小运算次数,避免误差积累。

二、数值实验（以下实验都需利用 Matlab 软件来完成） 1． 误差传播与算法稳定性 实验目的：体会稳定性在选择算法中的地位。误差扩张的算法是不稳定的，是我们 所不期望的；误差衰竭的算法是稳定的，是我们努力寻求的，这是贯穿本课程的目标。 -8- 数值分析实验指导书 实验内容：计算 En ? ? x n e x ?1dx, 0 1 n ? 1,2,? 算法一： E1 ? 1 e ， En ? 1 ? nEn?1 , n ? 2,3,? 算法二：

实验要求： EN ? 0 ， E n ?1 ? 1 ? En , n n ? N - 1, N - 2,?,3,2 （1） 分别用算法一、算法二采用 6 位有效数字计算 精确的结果。 En ，请判断哪种算法能给出更 （2） 请从理论上证明你实验得出的结果，解释实验的结果。设算法一中 E1 的计算 误差为 e 1 ， 由 E1 递推计算到 前递推计算到 En 的误差为 e n ；算法二中 E N 的计算误差为 ? N ，由 E N 向 En （ n ? N ）的误差为 ? n 。如果在上述两种算法中都假定后面的计算不再 引入其他误差，试给出 en 与 e1 的关系和 ? n 与 ? N 关系。 ?N 通 （3） 算法一中通常 e1 会很小，当 n 增大时， en 的变化趋势如何？算法二中 常相对比较大，当 n 减小时，误差 ? n 又是如何传播的？即比较两个算法，当某一步产生 误差后，该误差对后面的影响是衰减还是扩张的。

（4） 通过理论分析与计算实验，针对两个算法的稳定性，给出你的结论。

2． 不同方案收敛速度的比较 实验目的：通过实验体会数值计算中算法选择的重要地位。

实验内容：三种求㏑ 2 的算法比较。 ln 2 ? 1 ? 方案一：利用级数` ? 1 1 1 (?1) k ?1 ? ? ? ? ? .? 2 3 4 k ， k ?1 设 Sn ? ? -9- (?1) k ?1 k ， k ?1 n 数值分析实验指导书 则 ln 2 ? S n 方案二：对上述 ? Sn ? ? (?1) k ?1 k ， k ?1 n (s n ? sn?1 ) 2 sn ? S n ? (n ? 3,4,?) sn ? 2sn?1 ? sn?2 按 生成新数列 s n ，则 ln 2 ? s n ? 1 1 1 1 1 ? ? ? ? ? ? ? 2 3 4 k 3? 2 4?2 k ?1 k 2 方案三：利用级数 1 ? 2 2 ? 2 ? ? 设 实验要求： Sn ? ? 1 k k ?1 k 2 ，则 ln 2 ? S n n 分别用三种方案求出 ln 2 的近似值，要求 ?? 1 ? 10 ?5 2 ，观察比较三种计算方案的收 敛速度。在 MATLAB 命令窗口输入 log(2)求解，并与三种方案计算的结果进行比较。 - 10 - 数值分析实验指导书 第 2 章 非线性方程求解 一、主要要求 编写二分法、Newton 迭代法和快速弦截法通用子程序。

二、主要结内容回顾 1、二分法的基本思想：

二分法就是将方程的有根区间对分,然后再选择比原区间缩小一半的有根区间,如此 继续下去,直到得到满足精度要求的根为止的一种简单的区间方法。

基本法原理：给定方程 f(x)=0，设 f(x)在区间 [a，b] 连续,且 f(a)×f(b)<0，则 方程 f(x)在(a，b)内至少有一根，为便于讨论,不妨设方程 f(x)=0 在(a，b)内只有一实 根 x *采取使有根区间逐步缩小，从而得到满足精度要求的实根 x * 的近似值 x k 。

取[a,b]区间二等分的中点 x0 ? a?b 2 ，若 f(x0)=0,则 x0 是 f(x)=0 的实根； 若 f(a)f(x0)<0 成立,则 x * 必在区间(a, x0)内,取 a1=a，b1= x0；否则 x *必在 区间(x0,b)内，取 a1= x0，b1=b, 这样,得到新区间(a1,b1),其长度为[a,b]的一半, 如 此继续下去，进行 k 次等分后,得到一系列有根区间：

[a, b] ? [a1 , b1 ] ? ? ? [ak , bk ] ，其 中每一个区间长度都是前一个区间长度的一半，从而[ak ,bk ]的长度为 bk ? a k ? b?a 2 k 如此 继续下去,则有这些区间将收敛于一点，该点即为所求的根.当做到第 k 步时，有: | x * ? xk |? xk ? bk ? a k b ? a ? k ?1 ? ? 2 2 （ε a k ? bk 2 为给定的精度） 此时 即为所求方程的近似解。 2、二分法算法：

1）输入有根区间的端点 a,b 及预先给定的精度ε ； 2）计算(a+b)/2 并将其赋值给变量 c，记为(a+b)/2→c； 3）若 f(a)f(c)<0，则 c→b，转向 4） ，否则 c→a，转向 4） ； - 11 - 数值分析实验指导书 4）若 b-a<ε ，则输出满足精度要求的跟 c，否则转向 2） 3、二分法程序框图（见图 2） 4、不动点迭代法：

不动点迭代法又称简单迭代法，基本思想是构造不动点方程，以求得近似根。即由 方程 f(x)=0 变换为等价方程 x=?(x) ，这样原方程的根必满足：

x*= ?(x*) ，即 ?(x) 作用在 x*上，其值不发生变化，因此我们也称 x*为 ?(x) 的 不动点，要求方程 f(x)=0 的根就转化为求 ?(x) 的不动点了。

具体迭代如下：

先取一个估计值 x0 来试探，若 ?(x0) =x0，则 x*=x0(可能性很小) 一般 ?(x0) ≠x0，记 x1=?(x0) ，若 x1=?(x1) , 则 x*=x1 若 x1≠?(x1) ， 记 x2=?(x1) ，再用 x2 继续试探?? xk+1=?(xk) ，(k=0,1,2,?) 如此反复计算，即形成一迭代公式 如果{xk}收敛，则称迭代公式是收敛的；否则称迭代公式是发散的。

如果{xk}收敛于 x*，而 ?(x)是连续函数时，那么 x*即是 ?(x)的不动点。也即 x*就 是方程的根。 - 12 - 数值分析实验指导书 开始 输 入 a,b,e y1= f(a),y2=f(b ) 是 输出失败信 息 否 y1* y2>0 n=1 c=(a+b)/2，y=f(c) 否 a=c，y1=y y1*y<0 是 b=c y2=y 图2 否 n=n+ 1 输出x sto p b-a<e 是 5、Newton 迭代法的思想： xk ?1 ? xk ? 给定方程 f(x)=0，迭代公式为：

就称为牛顿迭代法。

6、Newton 迭代法的算法：

1）给出初始近似根 x0 及精度ε ； f ( xk ) f '

( xk ) ，应用该公式来解方程的方法 x0 ? 2）计算 f ( x0 ) f '

( x0 ) →x1； 3）若|x1-x0|<ε 或迭代次数达到预定指定的次数 N ，则转向 4） ；否则 x1→x0，转 向 2） ； 4）输出满足精度的根 x1，结束。

（程序框图略） 7、收敛性判定定理： - 13 - 数值分析实验指导书 定理 1：假定函数 ? (x) 满足下列条件：

1、对任意 x∈[a,b] 有 ?(x) ∈[a,b] 2、存在正数 L<1，使对任意 x1，x2∈[a,b] 有 |?(x1) - ?(x2)|≤L|x1-x2| (0≤L<1) 则迭代过程 xk+1= ? (xk)对于任意初值 x0∈[a,b]均收敛于方程 x= ?(x)的根 x* ， 且有如下的误差估计式： Lk xk ? x ? x1 ? x0 1? L * 局部收敛性定义：

如果存在不动点 x*的某个邻域 U(x*,δ )，使得对于任意初值 x0 ∈ U(x*,δ )，迭代公式 xk+1= ?(xk) （k=0,1,2......)产生的数列{xk}均收敛与 x*， 则称迭代公式 xk+1= ?(xk)是局部收敛的。

定理 2：设 ?(x)在 x= ?(x)的根 x*邻近有连续的一阶导数， 且| ?’(x*)|<1， 则迭代公式 xk+1=?(xk)具有局部收敛性。

8、收敛速度定义 | ek ?1 | ?C p 定义 当 k→∞时，有 | ek | （C≠0,且为常数） 则称迭代过程是 p 阶收敛的。

特别地，当 p=1，0<C<1 时，称为线性收敛； p=2 称为平方收敛 其中：ek=(x*-xk)为迭代误差 。

定理 3：对于迭代过程 xk+1= ?(xk) ，如果 ?(p)(x) 在所求根 x*的邻近连续，并且 ?’(x*)= ?’’(x*) =...= ?(p-1)(x*) =0，?(p)(x*)≠0，则该迭代过程在点 x* 邻近 是 P 阶收敛的。

注：牛顿迭代公式在根 x* 的邻近至少是平方收敛的。

三、数值实验（以下实验都需利用 Matlab 软件来完成） - 14 - 数值分析实验指导书 实验 2.1（求非线性方程根实验（一）） 实验目的：会求非线性方程的根 实验内容：1、用 Matlab 软件做出下列方程的图形，观察根所在的区间： 3x 2 ? e x ? 0 2、用二分法求上述方程的根，并分析其误差。

实验 2.1（求非线性方程根实验（二）） 实验目的：会求非线性方程的根 实验内容：设有方程：

3x ? e ? 0 2 x 1、分别用 Newton 迭代法和快速弦截法求上述方程的根； 2、比较两种迭代法的收敛速度。

具体操作见下例。

（以下实验都需利用 Matlab 软件来完成） 例 1：用二分法求方程 x2-x-1=0 的正根，要求误差小于 0.05 解：1）首先根据程序框图编制二分法的函数文件：

erfen.m 2）再编写函数文件 fc.m funtion y=fc(x) y=x^2-x-1;

3）最后在 Matlab 命令窗口输入调用命令：erfen(‘fc’,0,10,0.05) 4）输出结果为 n= 56 ans= 1.6180 例 2：求解非线性方程 5x2sinx-e-x=0，观察知有多解，如何求出所有解？ 解：1）编写 M—文件:newton.m 2）编写函数文件：fc1.m 和 df.m %fc1.m function y=fc1(x) y=5*x^2*sin(x)-exp(-x);

%df.m function y=df(x) y=10*x*sin(x)+ 5*x^2*cosx+exp(-x);

3）用命令 fplot(‘[5*x^2*sin(x)-exp(-x),0]’,[0,10])）作图，%在指定的范围内画出 函数的图形。

（注意：[5*x^2*sin(x)-exp(-x),0]中的‘[…，0]’是作 y=0 直线，即 x 轴） （见图 3） 从图中可看出方程在[0，10]区间有 4 个解，分别在 0、3、6、9 附近，所以 4）用命令 >>

>>

newton(x0) - 15 - 数值分析实验指导书 即可求出 求方程 fc1=0 在 x0 附近的近似解 得出结果 ans = 0.5018 3.1407 6.2832 9.4248 图3 - 16 - 数值分析实验指导书 第 3 章 线性方程组的解法 一、主要要求 编写方程组求解的通用程序：Jacobi 迭代、Seidel 迭代以及 SOR 迭代程序 二、主要结果回顾 1、迭代法：

? 它的基本思想是将线性方程组 AX=b 化为 ：

X=BX+f ? 再由此构造向量序列{X(k)}：

X(k+1)=BX (k)+f ? 若{X(k)}收敛至某个向量 x *,则可得向量 X *就是所求方程组 AX=b 的准确解. ? 线性方程组的迭代法主要有 Jocobi 迭代法、Seidel 迭代法和超松弛(Sor)迭代法. 2、收敛性判定定理：

定理 1：对任意初始向量 X(0)及任意向量 f，由此产生的迭代向量序列{X(k)}收敛的 充要条件是 ? ( B) ? 1. 。

定理 2：若||B||<1,则迭代公式 X(k)=BX(k-1)+f 收敛. ( k ?1) ? 3、Jacobi 迭代公式：

xi n a bi ij ?? x (k ) aii j ?1 aii j ?i i ? 1,2,?, n Jacobi 迭代的矩阵格式: X ( k ?1) ? BJ X ( k ) ? f J 其中：BJ = -D-1(U+L)，fJ=D-1b。

4、Seidel 迭代公式：

xi( k ?1) ? bi i ?1 aij ( k ?1) ? ? xj ? aii j ?1 aii j ?i ?1 ii ?a n aij x (jk ) i ? 1,2,?, n Seidel 迭代的矩阵格式: X ( k ?1) ? BS X ( k ) ? f S 其中：Bs= -(D+L)-1U，fs=(D+L)-1b。

5、Jacobi 迭代、Seidel 迭代的算法：

Jacobi 迭代算法：

1)给出矩阵 A，b，x0 2)计算 B= D-1(-U-L), f=D-1b 3) y=BX 0+f 4)若 ||y-x0||<=1.0e-6，则转到 5） ，否则，转到 3） 5）输出 y 和 n。 - 17 - 数值分析实验指导书 Seidel 迭代算法：

1)给出矩阵 A，b，x0 -1 -1 2)计算 B= (D-L) U，f=(D-L) b，则有：

3) y=BX 0+f 4)若 ||y-x0||<=1.0e-6，则转到 5） ，否则，转到 3） 5）输出 y 和 n。

6、程序框图：

（见图 4） 开始 输入a,b,x0 由矩阵a,生成D,U,L,B,f y=B*x0+f n=1 N ||y-x0||>=1.0e-6 Y x0=y y=B*x0+f N n=n+1 输出x 结束 图4 仿 Jacobi 迭代法算法和 Seidel 迭代算法可给出超松驰迭代的算法（略） 三、数值计算实验（以下实验都需用 Matlab 软件来完成） 实验 3.1（求解线性方程组实验） 实验目的：掌握各种迭代法，比较各种迭代法的渐进收敛速度. 实验内容：令 ? 25 ? 41 10 ? 6 ? ? ? ? ? 41 68 ? 17 10 ? A?? 10 ? 17 5 ? 3? ? ? ? ? 6 10 ? 3 - 18 ? 2 ? ? ? 32? ? ? ? 23? b?? ? 33 ? ? ? 31? ? ? 数值分析实验指导书 1、计算 A 的条件数 cond(A)； （可选用任何一种范数） 2、上述方程组是否为病态方程组？若是，如何改变方程组的病态性？ 3、分别用 Jacobi 迭代、Seidel 迭代与超松驰迭代求方程组 AX=b 的近似解； 4、比较 Jacobi 迭代、Seidel 迭代与超松驰迭代的渐进收敛速度。

注：上述实验分两次完成。

四、具体操作见下例（以下实验都需用 Matlab 软件来完成） 例：用 Jacobi 方法求解下列方程组，设 X(0)=0,精度为 10-6 ? 10x1 ? ? ? x1 ?? 2 x 1 ? ? x2 ? 10x 2 ? 2 x3 ? 10x3 ?9 ?7 ?6 解：1）先根据 Jacobi 算法编写 M—文件：Jacobi(a,b,X0) 2）在 Matlab 命令窗口输入命令：

>>a=[10 -1 0;-1 10 -2;0 -1 10];

>>b=[9;7;6];

>>Jacobi(a,b,[0;0;0]) 3）输出结果：y = 0.9938 0.9381 0.6938 n=9 注：

n=9 为迭代次数。 - 19 - 数值分析实验指导书 第 4 章 插值法 一、主要要求 编写拉格朗日插值法、分段线性插值法、*三次样条插值通用程序(Matlab 自带)。

拉格朗日插值多项式：

二、主要结果回顾 插值法是函数逼近的重要方法之一，有着广泛的应用 。在生产和实验中，函数 f(x)或者其表达式复杂不便于计算或者无表达式而只有函数在给定点的函数值(或其导 数值) ，此时我们希望建立一个简单的而便于计算的函数 ?(x)，使其近似的代替 f(x)， 这就是所谓的插值法 .有很多种插值法,其中以拉格朗日 (Lagrange)插值和牛顿(Newton) 插值为代表的多项式插值最有特点,常用的插值还有 Hermit 插值,分段插值和样条插值. 1、插值：求近似函数的方法:由实验或测量的方法得到所求函数 y=f(x) 在互异点 x0 , x1, ... , xn 处的值 y0 , y1 , … , yn , 构造一个简单函数 ?(x) 作为函数 y=f(x) 的近似表达式 ：

y= f(x) ? ?(x) 使 ?(x0)=y0 , ?(x1)=y1 , ?, ?(xn)=yn , (*) 这类问题称为插值问题。

f(x) 称为被插值函数， ?(x) 称为插值函数， x0 , x1, ... , xn 称 为插值节点。(*)式称为插值条件。

2、Lagrange 插值公式 Pn ( x) ? ? j ?0 n j ?0 n ( x ? x0)(x ? x1)...(x ? x j ?1)(x ? x j ?1)...(x ? xn) ( x j ? x0)( x j ? x1)...(x j ? x j ?1)( x j ? x j ?1)...(x j ? xn) n i ?0 i? j y j ? ? (? x ?x j x ? xi )y i j 其中 x0 , x1,... , xn 为插值节点，yj 为插值节点 xj 处的函数值（j=1,2,…n） 。

3、Lagrange 插值的截断误差(插值余项) 定理: 设 Ln(x)是过点 x0 ，x1 ，x2 ，…xn 的 n 次插值多项式， f (n)(x)在[a，b]上连 续，f (n+1)(x)在[a，b]上存在，其中[a，b]是包含点 x0 ，x1 ，x2 ，…,xn 的任一区间， 则对任意给定的 x?[a，b]，总存在一点 ??（a，b） （依赖于 x）使 R n ( x ) ? f ( x ) ? Ln ( x ) ? f ( n ?1) (? ) (n ? 1)! ? n ?1 ( x) 其中：

?n?1 ( x) ? ( x ? x0 )(x ? x1 )...(x ? xn ) 。 - 20 - 数值分析实验指导书 4、拉格朗日插值法程序框图：

（见图 5） 开始 输入xi ,yi i =0,1,2,...,n L=0,k=0 T=1 T=T(x-xi)/(xk-xi) i=0,1,2,...,k-1, k+1,...,n L=L+ykT k=n Y 结束 输出L （图 5） N k=k+1 三、数值计算实验（以下实验都需利用 Matlab 软件来完成） 实验 4.1（观察龙格(Runge)现象实验） 实验目的：观察拉格朗日插值的龙格(Runge)现象. 实验内容：

对于函数 f ( x ) ? 5 进行拉格朗日插值， 取不同的节点数， 在区间[-5,5] a ? x2 2 上取等距间隔的节点为插值点，把 f(x)和 n 次插值多项式的曲线画在同一张图上进行比 较。

（a 可以取任意值）具体步骤：

1、 a=1 时，1）取 n=4，作出 f(x)和插值多项式的曲线图； 2）取 n=10，作出 f(x)和插值多项式的曲线图； 2、a=0.5 时，1）取 n=4，作出 f(x)和插值多项式的曲线图； 2）取 n=10，作出 f(x)和插值多项式的曲线图； 3、分析上述曲线图，你可以得出什么结论？ - 21 - 数值分析实验指导书 四、具体操作 例 1、 给出 f(x)=lnx 的数值表，用 Lagrange 插值计算 ln(0.54)的近似值。

x Lnx 0.4 -0.916291 0.5 -0.693147 0.6 -0.510826 0.7 -0.357765 0.8 -0.223144 解：

1）首先根据程序框图拉格朗日插值法的函数文件：

lagrange(x0,y0,x) 2）在 Matlab 命令窗口输入调用命令：

>>x0=[0.4: 0.1: 0.8];

>>y0=[-0.916291 -0.693147 -0.510826 -0.356675 -0.223144];

>>lagrange(x0,y0,0.54) 3）输出结果为: -0.616143 (精确解: -0.616186) 实验 4.2 分段插值实验 实验目的：分段插值计算，会使用一维插值函数 ,寻找最佳的插值方法。

实验内容：在 1—12 的 11 小时内，每隔 1 小时测量一次温度，测得的温度依次为：5， 8，9，15，25，29，31，30，22，25，27，24。

1）试估计每隔 1/10 小时的温度值。

2）分别用分段线性插值，立方插值，三次样条插值和最邻近插值估计其值。 - 22 - 数值分析实验指导书 第 5 章 函数逼近与曲线拟合 一、主要要求 编写用最小二乘法求拟合数据的多项式 二、主要结果回顾 1.连续函数的最佳平方逼近 在 ? ? Span{1, x, x2 , , xn } 上，法方程为 H n a ? d ， 12 ? 1 ? 12 13 其中 H n ? ? ? ? ?1 (n ? 1) 1 (n ? 2) 均方误差： 1 (n ? 1) ? 1 1 (n ? 2) ? ? ， dk ? ( f , ?k ) ? f ( x)?k dx ?0 ? ? 1 (2n ? 1) ? 2 ? 2 ? ( f , f ) ? ( P* , f ) ? f ? max f ? P* 0? x ?1 ? ? ai*di 2 i ?1 n 最大误差：

? ? 2.离散函数的最佳平方逼近（曲线的最小二乘拟合） ：

法方程 ? (? ,? )a j ?0 j k m i ?0 n j ? ( f , ?k ) (? j , ?k ) ? ? ?i? j ( xi )?k ( xi ) 其中 ( f , ?k ) ? ? ?i f ( xi )?k ( xi ) i ?0 m 三、数值计算实验（以下实验都需利用 Matlab 软件来完成） 实验 5.1 （拟合实验） 实验目的：学会用最小二乘法求拟合数据的多项式，并应用算法于实际问题。

实验内容：给定数据点 ( xi , yi ) 如下： 0．6 0．7 - 23 - xi 0 0．5 0．8 0．9 1．0 数值分析实验指导书 yi 1 1．75 1．96 2．19 2．44 2．71 3．00 实验要求：

（ 1 ） 编写程序用最小二乘法求拟合数据的多项式，并求平方误差，作出离散函数 ( xi , yi ) 和拟合函数的图形。 （2） 用 MATLAB 的内部函数 polyfit 求解上面最小二乘法曲线拟合多项式的系数及平 方误差，并用 MATLAB 的内部函数 plot 作出其图形，并与（1）的结果进行比较。

5. 以正交多项式为基底，作最小二乘多项式拟合。

实验目的：探索改善最小二乘多项式的数值不稳定性的可能。

实验内容：

（1）在[-1，1]区间上取 n=20 个等距节点，计算出以相应点上 e 的值作为数据样本， 编写程序以 1， 其中， x P 1 ( x), P 2 ( x),?, P l ( x) 为基函数作出 l ? 3,5,7,9 次的最小二乘多项式， Pl ( x) 是勒让德多项式。画出 ln(cond( A)) 随着 n 变化的曲线，其中 A 是确定最 n 小二乘多项式系数的矩阵。计算出不同阶最小二乘多项式给出的最小偏差 ? (l ) ? ? ( y( xi ) ? yi ) 2 i ?1 （2） 用 MATLAB 的内部函数 polyfit 求解上面最小二乘法曲线拟合多项式的系数，并 用 MATLAB 的内部函数 plot 作出其图形。 - 24 - 数值分析实验指导书 第 6 章 数值积分 一、主要要求 编写定步长复合梯形公式、定步长复合 Simpson 公式、变步长梯形法及龙贝格算法*通 用子程序； 二、主要结果回顾 1、梯形公式：对于连续函数 f (x)，有积分中值定理 ? 如果取 f (? ) ? b a f ( x)dx ? (b ? a) f (? )(? ? (a, b)), 其中 f (ξ )是被积函数 f (x)在积分区间上的平均值。

因此， 如果我们能给出求平均值 f (ξ ) 的一种近似方法，相应地就可以得到计算定积分的一种数值方法。 f (a) ? f (b) ，即得下列梯形公式：

2 b (b ? a ) ?a f ( x)dx ? 2 [ f (a) ? f (b) ] , 二、辛普生(Simpson)公式：先用某个简单函数 近似逼近 f (x)，然后用 ? ( x) 在[a,b]区 间的积分值近似表示 f (x)在[a,b]区间上的定积分，即取 ? b a f ( x)dx ? ? ? ( x)dx a b 我们可以取 ? ( x) 为前面介绍的 f (x)的插值多项式 Ln(x)或拟合多项式 P(x)进行近似计 算。如取 ? ( x) 为插值多项式 Ln(x)，则相应得到的数值积分公式称为插值型求积公式。

如： ? b a f ( x)dx ? ? Ln ( x)dx a b 若考虑过点 A,B,C 三点的抛物线段 L 代替曲线段 y=f (x)（见图 7） - 25 - 数值分析实验指导书 y A ·· a a?b 2 C · b B o ? ( x) ? 取c ? x （图 7） L 就是以 a,b,c 为节点构造二次插值多项式： ( x ? c)(x ? b) ( x ? a)(x ? b) ( x ? a)(x ? c) f (a) ? f (c ) ? f (b) (a ? c)(a ? b) (c ? a)(c ? b) (b ? a)(b ? c) a?b 即得下列抛物线（辛普生）公式：

2 b b b?a a?b f ( x ) dx ? ?a ?a ? ( x)dx ? 6 [ f (a) ? 4 f ( 2 ) ? f (b)], 3、误差分析：

梯形公式的截断误差(余项)为: E1 ? ? (b ? a) 3 f '

'

(? ) (? ? (a, b)) 12 (b ? a) 5 f '

'

'

(? ) ( ? ? (a, b)) 2280 辛普生(Simpson)公式的截断误差(余项)：

E 2 ? ? 注：

1)抛物线求积公式的计算精度高于梯形求积公式的计算精度。

2)当积分区间较大时，利用上述公式计算积分产生的误差也较大。

所以在建立积分公式时，应选择小区间和抛物线求积公式。

? 从余项的讨论看到，积分区间越小，也可使求积公式的截断误差变小。因此，我 们经常把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式， 如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整 - 26 - 数值分析实验指导书 个区间上的求积公式，这就是复合求积公式的基本思想。

? 常用的复合求积公式是复合梯形公式和复合辛普生(Simpson)公式 复合梯形公式：

把区间[a,b] n 等分，取节点 xi=a+ih，i=0,1,...n, h=(b-a)/n，对每个小区间[xi,xi+1]用 梯形求积公式，再累加起来得： n ?1 h I ? [ f (a) ? f (b) ? 2? f ( xi )] ? T n 2 i ?1 5、复合辛普生(Simpson)公式：

把区间[a,b] n 等分，取节点 xi=a+ih， ，i=0,1,...n, h=(b-a)/n，对每个小区间[xi,xi+1]用 抛物线求积公式，再累加起来得： In ? （其中：

xi ? h n?1 h [ f ( xk ) ? 4 f ( xk ? ) ? f ( xk ?1 )] ? 6 k ?0 2 h 为节点 xi 与 xi+1 的中点，i=0,1,...n） 2 6、变步长复合积分法 基本思想：是在求积过程中，通过对计算结果精度的不断估计，逐步改变步长（逐 次分半） ，直到满足精度要求为止。即按照给定的精度实现步长的自动选取。

下面以梯形公式为例来说明这种方法。

利用梯形公式可得积分近似值：

记为 T0，即：

T0 ? ? b a f ( x)dx ? (b ? a ) [ f (a) ? f (b)], 2 b?a [ f (a ) ? f (b)], 2 (b ? a ) 2 此时步长 h0=b-a，再将[a,b]二等分，即取步长 h1 ? 使用 n=2 时的梯形公式可得积分近似值：

(b ? a) a?b T1 ? [ f (a) ? 2 f ( ) ? f (b)] 4 2 1 b?a b?a T0 ? f (a ? ), 2 2 2 (b ? a ) 再将[a,b]四等分，即取步长 h2 ? 22 ? 使用 n=4 时的梯形公式可得积分近似值： 3 (b ? a) i(b ? a) T2 ? { f (a) ? 2? f (a ? ) ? f (b)} 3 2 2 2 i ?1` 1 (b ? a) 2 2i ? 1 ? T1 ? f [a ? 2 (b ? a)], ? 2 2 2 2 i ?1 - 27 - 数值分析实验指导书 再将[a,b]八等分，即取步长 h3 ? (b ? a ) 23 使用 n=8 时的梯形公式可得积分近似值： 7 (b ? a) i(b ? a) T3 ? { f ( a ) ? 2 f ( a ? ) ? f (b)} ? 24 23 i ?1` 1 (b ? a) 2 2i ? 1 ? T2 ? f [a ? 3 (b ? a)], ? 3 2 2 2 i ?1 2 依次类推可得变步长梯形公式： 1 (b ? a) 2 2i ? 1 Tk ? Tk ?1 ? f [a ? k (b ? a)], k ? 1,2,... ? k 2 2 2 i ?1 即： k ?1 (b ? a ) [ f (a ) ? f (b)], 2 1 (b ? a ) b?a T1 ? T0 ? f (a ? ), 2 2 2 T0 ? 1 (b ? a) 2 2i ? 1 T2 ? T1 ? f [a ? 2 (b ? a)], ? 2 2 2 2 i ?1 1 (b ? a) 2 2i ? 1 T3 ? T2 ? f [a ? 3 (b ? a)], ? 3 2 2 2 i ?1 ?????????????????? 2 1 (b ? a) 2 2i ? 1 Tk ? Tk ?1 ? f [a ? k (b ? a)], k ? 1,2,... (*) ? k 2 2 2 i ?1 若预定精度为ε ，以公式(*)计算积分的近似值，直到|Tk+1-Tk|<ε 时终止计算,并以当前 值 Tk+1 为所求近似值. 对 Simpson 公式也可类似构造出相应的变步长积分公式。

（略） 7、算法：

a) 复合梯形公式：

I ? n ?1 h [ f (a) ? f (b) ? 2? f ( xi )] ? T n 2 i ?1 k ?1 b) 复合辛普生(Simpson)公式：

I n ? 算法： h n?1 h [ f ( xk ) ? 4 f ( xk ? ) ? f ( xk ?1 )] ? 6 k ?0 2 - 28 - 数值分析实验指导书 1）令 h ? b?a h , I 1 ? f (a ? ), I 2 ? 0 n 2 h ), I 2 ? I 2 ? f (a ? kh ) 2 2）对 k=1,2,...,n-1，计算 I 1 ? I 1 ? f (a ? kh ? 3） I ? h ( f (a) ? 4 I 1 ? 2 I 2 ? f (b)) 。

6 c）变步长梯形公式公式：

T2 n ? 算法：

1）输入精度ε ，n=1， T1 ? 2）s=0； h 1 Tn ? n 2 2 ? f (x 1 k? 2 ) ，这里 hn ? b?a 。

n b?a [ f (a) ? f (b)] ； 2 h ； 2 3）对 k=0,1,2,...,n-1 计算：

s ? s ? f ( x), x ? a ? k * h ? 4） T2 ? 1 (T1 ? h * s ) ，n=2n； 2 5）如果|T2-T1|<ε ，则结束；否则执行 6） ； 6） h ? h ，T1=T2，转 2） 。

2 三、数值计算实验（以下实验都需利用 Matlab 软件来完成） 实验 6.1 Newton-cotes 型求积公式 实验目的：学会 Newton-cotes 型求积公式,并应用该算法于实际问题. 实验内容：求定积分 ? ? e x cos xdx 0 实验要求：选择等分份数 n ,用复化 Simpson 求积公式求上述定积分的误差不超过 10?8 的近似值,用 MATLAB 中的内部函数 int 求此定积分的准确值,与利用复化 Simpson 求 积公式计算的近似值进行比较。

实验 6.2 Romberg 算法 实验目的：学会数值求积的 Romberg 算法,并应用该算法于实际问题. 实验内容：求定积分 实验要求：

（1）要求程序不断加密对积分区间的等分,自动地控制 Romberg 算法中的加速收敛 - 29 - ? 1 x dx 0.5 数值分析实验指导书 过程,直到定积分近似值的误差不超过 10 为止,输出求得的定积分近似值。

（2）可用 MATLAB 中的内部函数 int 求得此定积分的准确值与 Romberg 算法计算的 近似值进行比较。

实验 6.3 Gauss 型求积公式 实验目的：学会 Gauss 型求积公式,并应用该算法于实际问题. ?6 实验内容：求定积分 实验要求： ? dx 2 ?4 1 ? x 4 （1） 把 Gauss 点的表格存入计算机,以 Gauss-Legendre 求积公式作为本实验的例子, 要求程序可以根据不同的阶数 n ,自动地用 n 阶 Gauss-Legendre 求积公式计算上述定积 分的近似值.体会 Gauss 型求积公式是具有尽可能高的代数精度的数值求积公式。

（2）可用 MATLAB 中的内部函数 int 求得此定积分的准确值与 Gauss 型求积公式求 得的值进行比较。 - 30 - 数值分析实验指导书 第 7 章 矩阵特征值与特征向量的计算 一、主要要求 编写幂法和反幂法的程序。

二、主要结果回顾 1.幂法 幂法是求实方阵 A 按模最大的特征值及相应的特征向量的一种迭代方法。在学习时， 同学们应注意对前两种情况的公式的理解和记忆。

基本思想：任取非零初始向量 X0，作迭代序列 Xk+1=AXk，k=0,1,?. 再根据 k 增大时， Xk 各分量的变化规律，求出方阵 A 的模最大的特征值及相应的特征向量。

幂法的计算公式 设矩阵 A 的 n 个特征值按模的大小排列为： ?1 ? ?2 ? ?3 ? ?n ，其相应的特征 向量为 e1, e2,?, en，且线性无关。

任取初始非零初始向量 X0，作迭代序列 Xk+1=AXk， k=0,1,?. X0=a1e1+a2e2+?+anen .根据特征值λ 的不同情况得下列三种计算公式 （1）λ 1 为实根，且│λ 1│>│λ 2│。当 a1 不为 0，k 充分大时，则有 xi?k ?1? Xk ?1 ? ?k ? , e1 ? xi Xk ? （2）λ 1 为实根，且λ 1=-λ 2，│λ 2│>│λ 3│。当 a1 ，a2 不为 0，k 充分大时，则 有 1 ? xi? k ? 2 ? ? 2 X k ?1 ? ?1 X k ?? ? ? ? ?k ? ? , e1 ? ? ? 1 ? X k ?1 ? ?1 X k ? ? xi ? ? , ? 1 ? xi? k ? 2 ? ? 2 X k ?1 ? ?2 X k ?? ? ?? ? ? , e2 ? 2 ? ? k ? ? ? X k ?1 ? ?2 X k ? ? xi ? ? 在实际应用幂法时，可根据迭代向量个分量的变化情况判断属于那种情况。若迭代 向量各分量单调变化，且有关系式 Xk+1≈cXk，则属于第 1 种情况；若迭代向量各分量不是 单调变化，但有关系式 Xk+2≈cXk，则属于第 2 种情况； 在应用乘幂法计算特征值和特征向量时，为了防止溢出，也可采用迭代公式（6.6） - 31 - 数值分析实验指导书 进行迭代计算。当 A 为对称矩阵时，可采用计算内积的方法加速迭代的收敛，即计算特 xi (Xk , Xk ) 征值λ 1 时可用公式 ?1 ? ，其精度与上式（1）中 ?1 ? k ( X ?1k , X k ) xi? ? 幂法的计算步骤： T ? k ?1? 相当。 ①任取初始非零初始向量 X0（一般取 X0=（1，1，1） 或 X0=（1，0，0） ） 。

②作迭代序列 Xk+1= AXk，k=0,1,?（也可以列表计算） 。

③根据迭代向量个分量的变化情况判断属于那种情况，选择所属公式。

④代入公式计算出方阵 A 按模最大的特征值及相应的特征向量。

2.反幂法 反幂法是求实方阵 A 按模最小的特征值及相应的特征向量的一种迭代方法。

基本思想：设非奇异矩阵 A 的 n 个特征值为：λ 1≥λ 2≥?≥λ n，其相应的特征向量 为 e1, e2,?, en，则的特征值为 -1 T 1 ?1 ? 1 ?2 ? ??? ? 1 ?n 其相应的特征向量仍为 e1, e2,?, en。 -1 则求 A 按模最大的特征值的倒数则为矩阵 A 按模最小的特征值。

再利用幂法求 A 按模最 大的特征值。

反幂法的计算公式 任取初始非零初始向量 X0，作迭代序列 Xk+1= A Xk，k=0,1,?.它等价于 AXk+1=Xk， k=0,1,?.求得 Xk+1 .当│λ n-1 -1 │>│λ n│，a≠0,k 充分大时，则有 xi?k ? X k ?1 ?n ? ?k ?1? , en ? xi X k ?1 ? 反幂法的计算步骤：

（和幂法的计算步骤基本相同） ①任取初始非零初始向量 X0（一般取 X0=（1，1，1） 或 X0=（1，0，0） ） 。

②作迭代序列 Xk+1= A Xk，k=0,1,?（也可以列表计算） 。在实际计算中，为了减少运 算量，先将矩阵 A 作三角分解 A=LR，然后再求解方程组 -1 T T ? ?LYk ? X k k ? 0,1,? ? ?. ? ? RX ? Y k ? k ?1 ③根据迭代向量个分量的变化情况判断属于那种情况，选择所属公式。

④代入公式计算出特征值及相应的特征向量。 - 32 - 数值分析实验指导书 利用反幂法还可以求得比已知特征值 ? 更准确的特征值和特征向量。

求在 ? 附近的特征值。设与 ? 最接近的特征值为 ?i 用反幂法求出矩阵 A ? ? I 的按 模最小的特征值和相应的特征向量为 ?i ? ?i ? ? ? ~ 于是得 A 在 xi? xi? k ? k ?1? , ei ? ? 附近的特征值和相应的特征向量为 X k ?1 X k ?1 ? ?i ? ? ? ? i ? ? ? ~ ~ xi?k ? X k ?1 ? k ?1? , ei ? xi X k ?1 ? 三、数值计算实验（以下实验都需利用 Matlab 软件来完成） 实验 7.1 矩阵特征值的求法 实验目的：通过实验进一步熟悉掌握求矩阵特征值的各种方法的使用与比较。 ? 1 2 1 2? ? ? ? ? 2 1 0? ? 2 2 ?1 1? ? ? A?? B ? ? ? 6 0 2? 1 ?1 1 1? ? ? ? 1 1 3? ?2 1 1 1? ? ? ? ? 实验内容：给定矩阵 ? 4 ?1 ? ? ? ??1 4 ?1 ? ? ? ?1 4 ?1 ? C ?? ?1 4 ?1 ? ? ? ? 1 4 ? 1? ? ? ? ?1 4 ? ? ? 实验要求：

（1） 分别用幂法、反幂法、雅可比方法编写程序求矩阵 A、B、C 的按模最大特征 值、按模最小特征值、全部特征值（要求误差＜10 ?5 ） 。 （2） 用 MATLAB 的内部函数 eig 求矩阵 A、B、C 的全部特征值，并与（1）的结果比 较。

实验 7.2 用 QR 算法求矩阵的特征值 实验目的：通过实验进一步熟悉掌握求矩阵特征值的 QR 方法及原理。 - 33 - 数值分析实验指导书 实验内容：给定矩阵 实验要求： ? 6 2 1? ? ? A ? ? 2 3 1? ? 1 1 1? ? ?, ?2 ? ?4 H ? ?0 ? ?0 ?0 ? 3 4 5 6? ? 4 5 6 7? 3 6 7 8? ? 0 2 8 9? 0 0 1 0? ? （1） 根据 QR 算法原理编写程序求矩阵 A 及矩阵 H 的全部特征值（要求误差＜ 10 ?5 ） 。 （2） 直接用 MATLAB 的内部函数 eig 求矩阵 A 及矩阵 H 的全部特征值，并与（1）的结 果比较。 - 34 - 数值分析实验指导书 第 8 章 常微分方程数值解 一、主要要求 编写 Euler 法、改进的 Euler 法、四阶 Runge-Kutta 法(Matlab 自带)通用子程序； 二、主要结果回顾 解一阶常微分方程初值问题 ? y'

? f ( x, y) a ? x ? b ? y ( x0 ) ? y ? 将区间[a,b]作 n 等份，取步长 h ? b?a n 显式 Euler 公式：

yk ?1 ? yk ? hf ( xk , yk ) 隐式 Euler 公式：

y k ?1 ? y k ? h ( f ( x k , y k ) ? f ( x k ?1 , y k ?1 )) 2 y k ?1 ? y k ? hf ( xk , y k ) ?~ ? 改进的 Euler 公式：

? h y k ?1 ? y k ? ( f ( xk , y k ) ? f ( xk ?1 , ~ y k ?1 )) ? 2 ? ? ? y p ? y i ? hf ( xi , y i ) ? ? ? y c ? y i ? hf ( xi ?1 , y p ) ? ? y i ?1 ? 1 ( y p ? y c ) ? 2 ? 或表示为： 三、数值计算实验（以下实验都需利用 Matlab 软件来完成） 实验目的：比较几种数值解法的精度 实验内容：用 Euler 法和 Runge-Kutta 法求解初值问题： 2x ? ? y'

? y ? y ? ? y ( 0 ) ? 1 ? - 35 - (0 ? x ? 1) 数值分析实验指导书 四、具体操作 例：设有微分方程： ? y , ? x ? y x0 ? x ? T ? ? y ( x0 ) ? y 0 1） 求出上述微分方程的解析表达式（精确解） ； 2） 用 Euler 法解微分方程（步长取 h） 。

3） 将两方法的结果显示在一张图上，比较其精确度。

给出数据：x0 , T, y0 , h 输出(显示器)：数值解。

如 给出：0 输出 x 0.1 0.2 0.3 0.4 0.4 1 0.1 y 1.10000 1.22000 1.36200 1.52820 x 解：1）由微分方程可得精确解为：

y ? 2 * e ? x ? 1 2）先编写 M—文件 euler(a,b,x0,y0,h,n) 再根据题目要求调用函数 Euler(0.1,0.4,0,1,0.1,4) 结果：

ans =1.1000 1.2100 1.3310 1.4641 3）再 Matlab 命令窗口输入 x ? [0.1 0.2 0.3 0.4]； y ? 2 * e ? x ? 1 x 得 y ? [ 1.1103 1.2428 1.3997 1.5836] 1.4909] 而用 Euler 法计算得 y1 =[1.1050 1.2210 1.3492 作图：在 Matlab 命令窗口输入：

plot(x,y) hold on plot(x,y1,'-r') - 36 - 数值分析实验指导书 从图中可看出，用 Euler 法计算出的结果与精确值有较大误差，还需改进。 - 37 - 数值分析实验指导书 第四部分 MATLAB 入门 一、MATLAB 简介 MATLAB 全名叫 Matrix Laboratory，是矩阵实验室的意思。

MATLAB 最初是由 CleveMoler 用 Fortran 语言设计的， 现在的 Matlab 是由 MathWorks 公司用 C 语言开发的。

MATLAB 容易使用、可以由多种操作系统支持、具有丰富的内部函数、具有强大的 图形和符号功能、可以自动选择算法、与其他软件和语言有良好的对接性。MATLAB 是 适用于科学和工程计算的数学软件系统。

1．用户界面介绍 MATLAB 用户界面如图 1 所示，包括主菜单、主工具栏、命令窗口、命令历史窗口、 工作间管理窗口、当前路径窗口、编译窗口、图形窗口、启动按钮、帮助浏览器等。下 面简单介绍命令窗口，其他的内容在使用的过程中会逐渐熟悉。

命令窗口是用于输入数据，运行 MATLAB 函数和脚本并显示结果的主要工具之一， 命令窗口如图 1 右边中下区域所示。 - 38 - 数值分析实验指导书 图 1 MATLAB 用户界面 2．常量 常量是 MATLAB 语言预定义的一些变量，在默认的情况下这些变量的值为常数。下 面介绍常用的一些常量：

pi 表示圆周率； realmax 表示 MATLAB 可以表示的最大浮点数，MATLAB 可以表示的最大浮点数 为 ； Inf 表示无穷大，超过 MATLAB 可以表示的最大浮点数时，系统将视该数为无穷大； realmin 表示 MATLAB 可以表示的最小的正浮点数，MATLAB 可以表示的最小正浮 点数为 ； eps 表示用来判断是否为 0 的误差限，一般情况下，MATLAB 函数的默认误差限为 eps，其值为 ； NaN（Not a number）表示不定值，例类似 0/0，inf/inf 所生成的结果； i 或 j 表示纯虚数单位。

3．变量 变量是 MATLAB 的基本元素之一，与其他常规设计语言不同的是 MATLAB 语言不 要求对所使用的变量进行事先说明，而且他也不需要指定变量的类型，系统会根据该变 量被赋予的值或是对该变量所进行的操作来自动确定变量的类型。

在 MATLAB 语言中，变量的命名有如下规则：

（1） 变量名长度不超过 31 位，超过 31 位的字符系统将忽略不计； （2） 变量名区分大小写； （3） 变量名必须以字母开头，变量名中可以包含字母、数字或下划线。

注：

（1） 使用变量不要求事先声明； - 39 - 数值分析实验指导书 （2） 也存在变量的作用域问题，一般视为局部变量，仅在其调用的 M 文件内有 效。若要定义全局变量，应加上关键字 global，并一般用大写； （3） 如果在对某个变量赋值时，如果该变量已经存在，系统则会自动使用新值 来代替该变量的旧值。例如在命令窗口中输入“a=1;a=2”命令，则会得到 如下结果：

a= 2 4．算术运算 MATLAB 中用“+” 、 “-” 、 “*” “/”和“^”分别表示算术运算中的加、减、乘、除 和乘方。 1 ? 256 ? ( ) 2 ? 5 3 ? 2 3 3 例如：计算 6 1 1 >>256^(1/6)*(1/3)^(-1/2)+5^(1/3)/2^3 ans = 4.5782 5．数字的输入、运算 （1） （2） 对于简单的数字运算，可以直接在命令窗口以惯用的形式输入。

当表达式较复杂或重复次数较多时，可先定义变量，然后由变量表达式计 算。

（3） （4） （5） 若用户没有对表达式设定变量，则 MATLAB 自动将当前结果给 ans 变量。

%以后的内容只起注释的作用。

若不想显示中间的结果，可用‘； ’结束一行；若想再次察看，只需输入变 量名。

（6） （7） 乘幂和开方运算分别由^和函数 sqrt 实现。

数据存储和运算都以双精度进行的。 - 40 - 数值分析实验指导书 6．输出格式 MATLAB 中数值有多种显示形式。

（1） 缺省情况下，若数据为整数，则以整型表示；若为实数，则以保留小数点 后的 4 位浮点数表示。

（2） 输出格式可由 format 控制，该命令只影响在屏幕上的显示结果，不影响在 内部的存储和运算。

可结合 short （4 位小数） 、 long （14 位小数） 、 hex （16 进制） 、bank（ 2 进制） 、short e （4 位小数的指数形式） 、long e（ 14 位 小数） 、rational（分数形式） 。例在命令窗口中输入 format long ，则以 14 位小数显示结果；输入 format rational，则以 14 位小数显示结果。

7．用 MATLAB 编程 将一个完整的命令集合写入 M 文件便是一段 MATLAB 程序，但要注意，编程是在 MATLAB 的编辑窗口而不是命令窗口。

MATLAB 还提供了一般程序语言的基本功能。

(1) for 循环语句 for i=1:n for j=1:m A(i,j)=sin((i+j)/(m+n))+B(i,j);

end end 循环中的步长是可以选择的。如 for i=n:-2:n/2 A(i)=sin(i+n);

end 即循环控制变量从 n 开始，步长是-2 到 n/2 结束。 - 41 - 数值分析实验指导书 (2) while 循环语句 与计算机的其他高级语言一样，while 循环语句要由关系运算和逻辑运算给出的逻辑 控制，该语句的一般形式为 while（逻辑表达式） （一族可执行语句） end MATLAB 中的关系运算有 ＝＝ <＝ >＝ 相等 小于等于， 大于等于， ～＝ 不等， <

小于， >

大于. 而 MATLAB 提供的逻辑运算有 &

与， ∣ 或， ～ 非。 (3) 条件语句 for i=1:n for j=1:m if i==j A(i,j)=1;

else if(i<j)&(i>j/2) A(i,j)=-1;

else A(i,j)=0;

end end end 注意条件语句是以 end 结尾的。 - 42 - 数值分析实验指导书 (4) 内部函数 MATLAB 提供了丰富的运算函数，只需正确调用其形式就可得到满意的结果，常用 的运算函数如表 1 所示。 表 1 MATLAB 常用函数表 函数名 sin cos tan cot sec csc asin acos atan acot asec acsc asinh 函数功能 正弦 余弦 正切 余切 正割 余割 反正弦 反余弦 反正切 反余切 反正割 反余割 反双曲正弦 函数名 acosh atanh acoth asech acsch pow2 exp log log10 log2 sqrt lim fsolve 函数功能 反双曲余弦 反双曲正切 反双曲余切 反双曲正割 反双曲余割 以 2 为底的幂函数 以 e 为底的幂函数 自然对数 以 10 为底的对数 以 2 底的对数 平方根 求极限 求非线性方程的根 函数名 diff int dsolve dot cross fix floor ceil round mod rem sign abs 函数功能 求导 积分 求解微分方程 向量点乘 向量叉乘 向零方向舍入 向负方向舍入 向正方向舍入 四舍五入 有符号求余 无符号求余 符号函数 绝对值 (5) 用户自定义函数 MATLAB 允许用户使用 M 文件定义函数，但必须符合一定的规则。例如下面一段 程序存为 funsim.m function p=funsim(x) % define a simple function p=sqrt(x)-2*x^3+cos(x); - 43 - 数值分析实验指导书 所给出的是函数 x ? 2x ? cos(x) 。函数的输出可以多于一个， 3 例如 function [x1,x2]=quadroot(a,b,c) % solve quadratic equation ax^2+bx+c=0. ds=sqrt(b*b-4*a*c);

x1=(-b+ds)/2/a;

x2=(-b-ds)/2/a;

将它存为 quadroot.m，它给出二次方程的根。用户如此定义的函数可以与 MATLAB 内部函数同样使用. 8. 用 MATLAB 绘制图形 计算的可视化是当今科学与工程计算的主导方向之一，MATLAB 提供了许多可以选 用的图形功能，简介如下：

(1) 二维图形函数 plot 它是最常用和最简单的绘图命令。

例如 plot(x,y,‘—’)； 将向量 x 和 y 对应元素定义的点一次用实线联接 （的维数必须一样） ； 如果 x 和 y 为矩阵， 则按列依次处理。

plot(x1,y1,‘*’ ， x2,y2,‘+’)； 将向量 x 和 y 对应元素定义的点用星号标出，将向量 x2 和 y2 对应元素定义的点用‘+’ 标出。即 MATLAB 可以划线或者点，它提供的点和线类型如表 2 所示。

表 2 点和线对应符号 线 实线 虚线 符号 — -- 44 - 点 实心圆点 加号 符号 . + 数值分析实验指导书 点 虚线间点 : -. 星号 空心圆点 叉号 * 。

× 另外“fplot”用于化已定义函数在指定范围内的图像，它与“plot”类似但差别在于 fplot 可以根据函数的性质自适应的选择取值点。

目前大部分的硬盘拷贝设备尚不能支持颜色，但彩色的显示已经是十分普及的，所 以彩色显示在计算可视化中是十分有益的。颜色在绘图中的使用方法与线形的控制方法 一样，常用的颜色有 蓝色“b” ,黄色“y” ，红色“ r” ，绿色“g” 例如 plot(x,y,,’ —r ’) 将 x 和 y 对应的元素定义的点依次用红色实线联接。

(2) 绘图辅助函数 利用这一族函数可以为画出的图像加上标题等内容。

title(′· · ·′)；在图形的上方显示′′中所指定的内容； xlabel(′· · ·′)； 将′′中所指定的内容标在 x 轴； ylabel(′· · ·′)； 将′′中所指定的内容标在 y 轴； grid； 在图上显示虚线格； text(x,y,′· · ·′)； 将′′中所指定的内容标在由 x,y 所指定的位置上； gtext(′· · ·′)； 运行到该命令，屏幕光标位置显示符号“+”等待，它将′′ 中所指定的内容标在鼠标所指定的位置； axis([x1 xr y1 yr])； 其中的 4 个实数分别定义 x 和 y 方向显示的范围； hold on；后面 plot 的图将叠加在一起； hold off；解除 hold on 命令，plot 将先冲去图形窗口已有图形； 注意上述辅助函数必须放在相应的“plot”之后。

(3) 多窗口绘图函数 subplot 该函数的形式为 subplot(p,q,r) - 45 - 数值分析实验指导书 该命令将图形窗口分成 p 行 q 列共计 p×q 个格子上画图，格子是从上到下按行 一次记数的。

例如考虑 Chebeshev 多项式，它可以用其递推公式定义如下： t 0 ( x) ? 1, t1 ( x) ? x, t k ( x) ? 2 xtk ?1 ( x) ? t k ?2 ( x), k ? 2,3,? 下面的程序将[-1，1]上的前四个 Chebeshev 多项式画在一张图上。

% see how hold works Chebeshev P 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1 y T0 >>

x=-1:0.05:1;

>>

t0=1.0+0*x;t1=x; T3 T1 T2 >>

t2=2*x.*t1-t0;t3=2*x.*t2-t1;

>>

plot(x,t0);gtext('T0');

>>

title('Chebeshev P');

>>

xlabel('x');ylabel('y'); -0.8 -0.6 -0.4 -0.2 0 x 0.2 0.4 0.6 0.8 1 >>

hold on >>

plot(x,t1);gtext('T1'); >>

plot(x,t2);gtext('T2');

>>

plot(x,t3);gtext('T3');

>>

hold off 得到的图形将在图形窗口输出如图 2 所示的图形。如果将上述一段命令在文件编辑 窗口中依次录入，并存为 shd.m，需要运行该程序时，只要在命令窗口输入命令“shd” 即可。

图 2 重叠绘制出 Chebeshev 多项式 而使用 subplot 则程序为 % see how hold works Chebeshev >>

x=-1:0.05:1; - 46 - 数值分析实验指导书 >>

t0=1.0+0*x;t1=x;

>>

t2=2*x.*t1-t0;t3=2*x.*t2-t1;

>>

subplot(2,2,1);plot(x,t0);

>>

title('Chebeshev T0');

>>

xlabel('x');ylabel('y');

>>

subplot(2,2,2);plot(x,t1);

>>

title('Chebeshev T1');

>>

xlabel('x');ylabel('y');

>>

subplot(2,2,3);plot(x,t2);

>>

title('Chebeshev T2');

>>

xlabel('x');ylabel('y');

>>

subplot(2,2,4);plot(x,t3);

>>

title('Chebeshev T3');

>>

xlabel('x');ylabel('y');

得到的图形输出如图 3 所示。 Chebeshev T0 2 1.5 1 0.5 0 -1 1 0.5 0 -0.5 -1 -1 Chebeshev T1 y -0.5 0 0.5 x Chebeshev T2 1 y -0.5 0 0.5 x Chebeshev T3 1 1 0.5 0 -0.5 -1 -1 1 0.5 0 -0.5 -1 -1 y -0.5 0 x 0.5 1 y -0.5 0 x 0.5 1 图 3 重叠绘制出 Chebeshev 多项式 - 47 - 数值分析实验指导书 (4) 三维图形函数 三位的立体图形在 MATLAB 中形成也不难。例如对于函数 f ( x, y ) ? sin( x 2 ? y 2 ) x2 ? y2 下面的程序画出该函数在区域 ( x, y) ? [?18,18] ? [?18,18] 上的三维图形。 1 0.5 z 0 -0.5 80 60 40 20 y 0 0 20 x 60 40 80 图 4 三维图形 二、数值分析中的常用命令 数值分析主要研究数学问题的数值解，涉及的内容包括代数、微积分、微分方程等。

利用 MATLAB 提供的函数，可以实现以下问题的数值求解。 1．求解线性方程组 线性方程组的解法很多，主要介绍几种常用的方法。

(1) 基于矩阵变换的直接解法 ? 使用\或/运算符，可以基于一系列矩阵变换直接求解线性方程组； 例如：求解方程组 Ax ? b - 48 - 数值分析实验指导书 >>

A=[2,1,-5,1;1,-3,0,-6;0,2,-1,2;1,4,-7,6];

>>

b=[8,9,-5,0]';

>>

x=a\b x= 3.0000 -4.0000 -1.0000 1.0000 ? linsolve(A,b)：解线性方程组 Ax=b. 例如：求解方程组 Ax ? b >>

a=[2,1,-5,1;1,-3,0,-6;0,2,-1,2;1,4,-7,6];

>>

b=[8,9,-5,0]';

>>

x=linsolve(a,b) x= 3.0000 -4.0000 -1.0000 1.0000 (2) 迭代法 如 Jocabi 迭代法、Gauss-Seidel 迭代法、SOR（超松弛）迭代法需要根据算法编 写相应的 M 文件来实现。 2.求矩阵的特征值和特征向量 用 eig 或 eigs 函数求矩阵的特征值和特征向量。 - 49 - 数值分析实验指导书 ?1 2 0 ? ? ? A ? ? 2 5 ?1 ? ? 4 10 ?1 ? ? ? 的特征值和特征向量 例如：求矩阵 >>

A=[1 2 0;2 5 -1;4 10 -1];

>>

[b,c]=eig(A) b= -0.2440 -0.3333 -0.9107 c= 3.7321 0 0 0 0.2679 0 0 0 1.0000 -0.9107 0.3333 -0.2440 0.4472 0.0000 0.8944 其中，b 和 c 矩阵分别为特征向量矩阵和特征值矩阵. 3．求解一元非线性方程 (1) 用 fzero 函数求求一元非线性方程的零点。该函数的调用格式为：

x=fzero(fun,x0):如果 x0 为标量，则试图寻找函数 fun 在 x0 附近的零点。

例如：求方程 x ? 4 x ? 5 ? 0 3 >>

f=@(x)x.^3-4*x+5;

>>

z=fzero(f,1) z= -2.4567 (2) 用 roots 函数计算多项式的根，该函数的调用格式为：

r=roots(c):返回一个列矢量，其元素为多项式 c 的解。 - 50 - 数值分析实验指导书 例如：求方程 x ? 4 x ? 5 ? 0 3 >>

p=[1,0,-4,5];

（用行矢量表示多项式） >>

r=roots(p) r= -2.4567 1.2283 + 0.7256i 1.2283 - 0.7256i (3) newton 法需要编写 M 文件来实现。

非线性方程组的求解一般采用迭代法求解，主要有不动点迭代法、Newton 迭代法和 拟 Newton 迭代法，需要编程实现。 4. 插值 (1) 一维线性插值 函数的调用格式为：

yi=interp1(x,y,xi,’method’) 功能：对一组点（x,y)进行插值，计算 xi 的值。xi, yi 可以是矩阵。

method 的取值可以是：

? Nearest（线性最近项插值） ? Linear（线性插值，默认值） ? Spline（三次样条插值） ? Cubic（三次插值） 例如：

>>

x=[0.4:0.1:0.8]; - 51 - 数值分析实验指导书 >>

y=[-0.916291,-0.693147,-0.510826, -0.356675,-0.223144];

>>

interp1(x,y,0.54,'cubic') ans = -0.61610990887357 >>

interp1(x,y,[0.54 0.54],'nearest') ans = -0.69314700000000 -0.69314700000000 (2) 二维线性插值 函数的调用格式为：

zi=interp2(x,y,z,xi,yi’method’) 功能：对一组点（x,y,z)进行插值，计算(xi, yi )的值。

xi, yi 可以是矩阵，method 的取值同上。

(3) 三维线性插值 函数的调用格式为：

vi=interp2(x,y,z,v,xi,yi ,zi’method’) 功能：对一组点（x,y,z,v)进行插值，计算(xi, yi , ,zi)的值。

xi, yi zi 可以是矩阵，method 的取值同上。

(4) 三次样条插值 函数的调用格式为：

yi=spline(x,y,xi) 功能：消除高阶多项式的插值产生的病态。等同于 yi = interp1(x,y,xi,’spline’) 例如：

>>

x=[0.4:0.1:0.8];

>>

y=[-0.916291,-0.693147,-0.510826, -0.356675,-0.223144]; - 52 - 数值分析实验指导书 >>

yi=spline(x,y,0.5) yiz = -0.6931 5．曲线拟合 函数的调用格式为：

p=ployfit(x,y,n) 功能：利用离散的点来生成一条连续的曲线。已知一组数据（xi,yi）,从中找出自变量 x 与因变量 y 之间的函数关系 y=f(x)。并不要求 y=f(x)在每个点上都完全相等，它只要求 在给定点 xi 上使误差 ? [f(xi ) ? yi ] 最小。

例如：

>>

x=[1 3 4 5 6 7 8 9 10];

>>

y=[10 5 4 2 1 1 2 3 4];

>>

[p,s]=polyfit(x,y,4);

>>

y1=polyval(p,x);

>>

plot(x,y,'go',x,y1,'b--') 2 6．数值积分 (1) 定积分 ? 用 trapz 函数进行梯形数值求积. ? 例如：计算 0 πsin xdx >>

X=0:pi/100:pi;

>>

Y=sin(X);

>>

Z=trapz(X,Y) - 53 - 数值分析实验指导书 Z= 1.9998 ? 用 quad 和 quad8 函数进行自适应递归 Simpson 求积 ? 例如：计算 2 0 1 dx x ? 2x ? 5 3 首先编写函数的 M 文件 myfun.m function y=myfun(x) y=1./(x.^3-2*x-5);

然后在命令窗口输入 >>

Q=quad(@myfun,0,2) Q= -0.4605 ? Gauss 求积和 Romberg 求积需要编程序来实现。 (2) 二重积分 用 dblquad 函数对二重积分进行数值计算，该函数语法格式为：

q=dblquad(fun,xmin,xmax,ymin,ymax) q=dblquad(fun,xmin,xmax,ymin,ymax,tol) q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method) 例如：计算 ? ?? 0 ? 2? ( ysin(x)? xcos(y ))dxdy 首先编写函数的 M 文件 integrnd.m function z=integrnd(x,y) z=y*sin(x)+x*cos(y);

然后在命令窗口输入 >>

Q=dblquad(@integrnd,pi,2*pi,0,pi) - 54 - 数值分析实验指导书 Q= -9.8696 (3) 三重积分 用 triplequad 函数对三重积分进行数值计算，该函数语法格式为：

q=triplequad (fun,xmin,xmax,ymin,ymax,zmin,zmax) q=triplequad (fun,xmin,xmax,ymin,ymax, zmin,zmax,tol) q=triplequad (fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method) 例如：计算 ? ? ? 1 1 ? -1 0 0 ( ysin(x)? zcos(y))dxdydz 首先编写函数的 M 文件 integrnd.m function u=integrnd(x,y) u=y*sin(x)+z*cos(y);

然后在命令窗口输入 >>

Q=triplequad(@integrnd,0,pi,0,1,-1,1) Q= 2.0000 7．求解微分方程 用 dsolve 命令可以求出微分方程的解，该函数的调用格式为：

r=dsolve(‘eq1,eq2,…’’cond1,cond2,...’,’v’) dy ? 1 ? y2 dx 例如：求解微分方程 >>

dsolve('Dy=1+y^2') ans = tan(t+C1) - 55 - 数值分析实验指导书 若指定初始条件 y x?0 ? 1 >>

y=dsolve('Dy=1+y^2','y(0)=1') y= tan(t+1/4*pi) - 56 - 

张亚玲 数值分析实验指导书2013 暂无相关推荐文档 如要投诉违规内容,请到 百度文库投诉中心 ;如要提出功能问题或意见建议,请 点击此处 进行反馈. 暂无评价 | 0人阅读 | 0次...

数值分析实验指导书 数值分析实验指导书

《新乡市政协委员张亚玲,石家庄张亚玲,新乡张亚玲,张亚玲,数值分析实验指导书2013》出自:QQ空间素材网
链接地址：http://www.qzoneai.com/sucai/jeVUcRett1Of1fZ9.html