数值分析课程

数值分析课程
实验指导书
太原科技大学
应用科学学院数学系
目录
前言 (1)
第一部分数值实验报告格式 (1)
第二部分数值实验报告范例 (2)
第三部分数值实验 (6)
数值实验一 (6)
数值实验二 (8)
数值实验三 (10)
数值实验四 (12)
数值实验五 (13)
数值实验六 (16)
数值实验七 (17)
第四部分MATLAB入门 (19)
前言
该实验指导书是《数值分析》课程的配套数值实验教材。

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

数值实验是数值分析课程中不可缺少的部分，利用计算机进行数值实验，以消化巩固所学的内容，增加对算法的可靠性、收敛性、稳定性及效率的感性认识，体会和重视算法在计算机上实验时可能出现的问题。

学生通过选择算法、编写程序、分析数值结果、写数值实验报告等环节的综合训练，逐步掌握数值实验的方法和技巧，获得各方面的数值计算经验，培养学生运用所学算法解决实际问题和进行理论分析的能力。

该实验指导书由王希云、刘素梅、王欣洁、李晓峰等老师编写。

第一部分数值实验报告格式
一个完整的实验，应包括数据准备、理论基础、实验内容及方法，最终对实验结果进行分析，以达到对理论知识的感性认识，进一步加深对相关算法的理解，数值实验以实验报告形式完成，实验报告格式如下：
一、实验名称
实验者可根据报告形式需要适当写出。

二、实验目的及要求
首先要求做实验者明确，为什么要做某个实验，实验目的是什么，做完该实验应达到什么结果，在实验过程中的注意事项，实验方法对结果的影响也可以以实验目的的形式列出。

三、算法描述（实验原理与基础理论）
数值实验本身就是为了加深对基础理论及方法的理解而设置的，所以要求将实验涉及到的理论基础，算法原理详尽列出。

四、实验内容
实验内容主要包括实验的实施方案、步骤、实验数据准备、实验的算法以及可能用到的仪器设备。

五、程序流程图
画出程序实现过程的流程图，以便更好的对程序执行的过程有清楚的认识，在程序调试过程中更容易发现问题。

六、实验结果
实验结果应包括实验的原始数据、中间结果及实验的最终结果，复杂的结果可以用表格形式列出，较为简单的结果可以与实验结果分析合并出现。

七、实验结果分析
（对算法的理解与分析，包括改进与建议）
实验结果分析是数值实验的重要环节，只有对实验结果认真分析，才能对实验目的、实验方法进一步理解，对实验的重要性充分认识，明确数值分析的实用范围及其优缺点。

第二部分数值实验报告范例
为了更好地做好数值实验并写出规范的数值实验报告，下面给出一简单范例供读者参考。

数值实验报告
一、实验名称
误差传播与算法稳定性
二、实验目的
1．理解数值计算稳定性的概念。

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

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

三、实验内容
计算
dx
x x I n
n ⎰
+=1
10，1,2,,10n =
四、算法描述
由 dx x x I n
n ⎰
+=1
10，知 dx
x x I n n ⎰+=--101110
则
n dx x dx x x x I I n n n n n 11010101011
11=
=++=+⎰⎰
---
可得递推关系
1．=n I 1
101
--n I n ，10,,2,1 =n
2．
)1
(1011n n I n I -=
- ，1,,9,10 =n
下面分别以1，2递推关系求解
方案1 =n I 1
101
--n I n ，10,,2,1 =n
当0=n 时，
=
+=⎰
dx x I 1
0101
㏑=1011㏑1.1 , 递推公式为
⎪⎩⎪⎨
⎧
==-=-1
.1ln 10
,,2,1,10101I n I n I n n （1）
方案2
)1
(1011n n I n I -=
- ，1,,9,10 =n
当10<<x 时， n
n n x x x x 101
1011
1≤+≤ 则
dx
x dx x x dx x n
n n 1011011110101
⎰⎰⎰
≤+≤
即
)1(101)1(111+≤
≤+n I n n
取递推初值
)110(22021])110(101)110(111[2110+=
+++≈I 递推公式为⎪⎪⎩⎪⎪⎨⎧
+==-=-)110(220211,,9,10),1(101101I n I n I n n （2）
取递推公式（1）中的初值
095310.01.1ln 0≈=I ，得
⎪⎩⎪⎨
⎧
≈=-=-095310
.010
,,2,1 ,10101I n I n I n n
取递推公式（2）中的初值
008678.010≈I ，得
⎪⎩⎪⎨⎧
≈=-=-008678
.01,,9,10),1(101101
I n I n I n n
五、程序流程图
由于实验方案明显、简单，实现步骤及流程图省略。

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

由此看出，当n 较大时，用递推公式（1）中的n I ~
近似n I 是不正确的。

主要原因是初值095310.0~0=I 不是精确值，设有误差)~
(0I e ，由递推公式（1）知
)~
(10)~(1--=n n I e I e 则有
)~()10()~(100)~(10)~(021I e I e I e I e n n n n -=-=-=--
误差
)~(n I e 随n 的增大而迅速增加，增加到)~(0I e 的n )10(-倍。

由此可见，递推公式计算的误差不仅取决于初值的误差，公式的精确性，还依赖于误差的传递即递推计算的稳定性。

由递推公式（2）知
008678.010≈I ，n I 为估计值，并不精确，有
12101
)(10≤
I e ，
而由
)(101)(*
*
1n n I e I e -
=- 得
)()101()(*
*
0n n I e I e -
=
误差)(*
0I e 随递推公式逐步缩小。

综上所述，在递推计算中，数值计算方法是非常重要的，
误差估计、误差传播及递推计算的稳定性都会直接影响递推结果。

第三部分 数值实验
本部分列出了七个数值实验类型，每个数值实验都应在计算机上实现或演示，由实验者独立编程实现。

要求：
（1） 用MATLAB 语言或其他算法语言编程，使之尽量具有通用性。

（2） 上机前充分准备，复习有关算法，画出程序流程图。

（3） 完成实验后写出实验报告。

（4） 编程语言的种类、运行环境及程序清单以附录形式给出。

以下分别列出七个数值实验的内容。

数值实验一
1． 误差传播与算法稳定性
实验目的：体会稳定性在选择算法中的地位。

误差扩张的算法是不稳定的，是我们所不期望的；误差衰竭的算法是稳定的，是我们努力寻求的，这是贯穿本课程的目标。

实验内容：计算
1,2,n ,11
==-⎰dx e x E x n n
算法一：
e E 1
1=
， 2,3,n ,11=-=-n n nE E
算法二：
0=N E ，
,3,22,-N 1,-N n ,11 =-=
-n E E n
n
实验要求：
（1） 分别用算法一、算法二采用6位有效数字计算n E ，请判断哪种算法能给出更
精确的结果。

（2） 请从理论上证明你实验得出的结果，解释实验的结果。

设算法一中1E 的计算误差为e 1，由1E 递推计算到n E 的误差为e n ；算法二中N E 的计算误差为N ε，由N E 向
前递推计算到
n E （N n <）的误差为n ε。

如果在上述两种算法中都假定后面的计算不再
引入其他误差，试给出n e 与1e 的关系和
n ε与N ε关系。

（3） 算法一中通常1e 会很小，当n 增大时，n e 的变化趋势如何？算法二中N ε通
常相对比较大，当n 减小时，误差
n ε又是如何传播的？即比较两个算法，当某一步产生
误差后，该误差对后面的影响是衰减还是扩张的。

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

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

实验内容：三种求㏑2的算法比较。

方案一：利用级数`
∑
∞
=--=+-+-=11
)1(.41312112ln k k k ， 设
∑
=--=n
k k n k S 11
)1(， 则n S ≈2ln
方案二：对上述
∑
=--=n
k k n k S 11
)1(， 按)
3,4,(n 2)s (2
12
1n =+---=---∧
n n n n n n s s s s S s
生成新数列∧n s ，则∧
≈n s 2ln
方案三：利用级数 ∑∞
==+⋅+⋅+⋅+⋅143221241231221211k k k
设
∑
==n
k k
n k S 121，则n S ≈2ln
实验要求：
分别用三种方案求出2ln 的近似值，要求
51021
-⨯=
ε，观察比较三种计算方案的收
敛速度。

在MATLAB 命令窗口输入log(2)求解，并与三种方案计算的结果进行比较。

数值实验二
1. 解线性方程组的列主元素高斯消去法和LU 分解法
实验目的：通过数值实验，从中体会解线性方程组选主元的必要性和LU 分解法的优点，以及方程组系数矩阵和右端向量的微小变化对解向量的影响。

实验内容：解下列两个线性方程组
（1） ⎪
⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝
⎛--11134.981.4987.023.116.427.199.103.601.3321x x x （2） ⎪⎪⎪⎪
⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫
⎝
⎛----15900001.58201
2151526099999.2310
710
4321x x x x 实验要求：
（1） 用你熟悉的算法语言编写程序用列主元高斯消去法和LU 分解求解上述两个方程组，输出Ax=b 中矩阵A 及向量b, A=LU 分解的L 及U ，detA 及解向量x.
（2） 将方程组（1）中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解变换后的方程组，输出列主元行交换次序，解向量x 及detA ，并与（1）中结果比较。

（3） 将方程组（2）中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的方程组，输出解向量x 及detA ，并与（1）中的结果比较。

（4）用MATLAB 的内部函数inv 求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个方程组的解，并与列主元高斯消去法和LU 分解法求出的解进行比较，体会选主元的方法具有良好的数值稳定性。

用MATLAB 的内部函数det 求出系数行列式的值，并与（1）、（2）、（3）中输出的系数行列式的值进行比较。

2. 研究解线性方程组b Ax =迭代法收敛速度
实验目的：认识迭代法收敛的含义以及迭代法初值和方程组系数矩阵性质对收敛速度的影响。

实验内容：用迭代法求解b Ax =，其中20
20⨯∈R A 为五对角矩阵
20201132
411132*********
24111134
22
411134
221134
2
A ⨯⎛⎫-- ⎪ ⎪ ⎪---
⎪ ⎪ ⎪----
⎪ ⎪⋅
⋅⋅⋅ ⎪=⋅⋅⋅⋅
⎪ ⎪⋅⋅⋅⋅ ⎪ ⎪-
--- ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪-- ⎪
⎝
⎭
实验要求：
（1）选取不同的初始向量X
)
0(及右端向量b ，给定迭代误差要求，用Jacobi 迭代法
和Gauss-Seidel 迭代法求解，观察得到的序列是否收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论。

（2）用SOR 迭代法求上述方程组的解，松驰系数ω取21<<ω的不同值，在
5
)
1()(10-∞
+≤-k k X X 时停止迭代，记录迭代次数，分析计算结果与松驰系数ω的关
系并得出你的结论。

（3）用MATLAB 的内部函数inv 求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个方程组的解，并与上述方法求出的解进行比较。

3．病态线性方程组的求解
实验目的：通过实验体会病态线性方程组的性态，了解求解病态线性方程组的方法。

实验内容：考虑方程组b Hx =的求解，其中系数矩阵H 为Hilbert 阵，
n
j i j i h h H j i n n j i ,,2,1,,11
,)(,, =-+=
=⨯
这是一个著名的病态问题。

通过首先给定解（例如取各个分量均为1）再计算出右端的办法给出确定的问题。

实验要求：
（1）选择问题的维数为6，分别用列主元Gauss消去法、Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？
（2）逐步增大问题的维数，仍然用上述方法来解它们，计算的结果如何？计算的结果说明了什么？分析产生结果的原因。

（3）讨论病态方程组的解法，并求解维数为6、10、100的上述方程组。

数值实验三
1．矩阵特征值的求法
实验目的：通过实验进一步熟悉掌握求矩阵特征值的各种方法的使用与比较。

实验内容：给定矩阵
⎪⎪
⎪
⎪
⎪
⎭
⎫
⎝
⎛
-
-
=
1
1
1
2
1
1
1
1
1
1
2
2
2
1
2
1
A
⎪
⎪
⎪
⎭
⎫
⎝
⎛
-
-
=
3
1
1
2
6
1
2
B
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭
⎫
⎝
⎛
-
-
-
-
-
-
-
-
-
-
=
4
1
1
4
1
1
4
1
1
4
1
1
4
1
1
4
C
实验要求：
（1）分别用幂法、反幂法、雅可比方法编写程序求矩阵A、B、C的按模最大特征值、按模最小特征值、全部特征值（要求误差＜105-）。

（2）用MATLAB的内部函数eig求矩阵A、B、C的全部特征值，并与（1）的结果比较。

2．用QR算法求矩阵的特征值
实验目的：通过实验进一步熟悉掌握求矩阵特征值的QR方法及原理。

实验内容：给定矩阵
⎪
⎪
⎪
⎭
⎫
⎝
⎛
=
1
1
1
1
3
2
1
2
6
A
,
⎪
⎪
⎪
⎪
⎪
⎪
⎭
⎫
⎝
⎛
=
1
9
8
2
8
7
6
3
7
6
5
4
4
6
5
4
3
2
H
实验要求：
（1）根据QR算法原理编写程序求矩阵A及矩阵H的全部特征值（要求误差＜105-）。

（2）直接用MATLAB的内部函数eig求矩阵A及矩阵H的全部特征值，并与（1）的结果比较。

数值实验四
1．迭代函数对收敛性的影响
实验目的：初步了解非线性方程的简单迭代法及其收敛性，体会迭代函数对收敛性的影响，知道当迭代函数满足什么条件时，迭代法收敛。

实验内容：用简单迭代法求方程
1
2
)
(3=
-
-
=x
x
x
f的根。

方案一：化
1
2
)
(3=
-
-
=x
x
x
f为等价方程
)
(
2
1
3x
x
xφ
==
∆
+
=
方案二：化
1
2
)
(3=
-
-
=x
x
x
f为等价方程)
(
1
23x
x
xφ
==
∆
-
=
实验要求：
（1）分别对方案一、方案二取初值
=
x
，迭代10次，观察其计算值，并加以分
析。

（2）用MATLAB内部函数solve直接求出方程的所有根，并与（1）的结果进行比较。

2. 初值的选取对迭代法的影响
实验目的：通过具体的数值实验，体会选取不同的初值对同一迭代法的影响。

实验内容：用牛顿迭代法求方程 013
=--x x 在x =1.5附近的根。

实验要求：
（1）对牛顿迭代公式：
1312
3
1
----=+k k k k k x x x x x ，编写程序进行实验，分别取00
=x ，5.10=x 迭代10次，观察比较其计算值，并分析原因。

（2）用MATLAB 内部函数solve 直接求出方程的所有根，并与（1）的结果进行比较。

3． 收敛性与收敛速度的比较
实验目的：通过用不同迭代法解同一非线性方程，比较各种方法的收敛性与收敛速度。

实验内容：求解非线性方程 0232=-+-x
e x x 的根，准确到106
-。

实验要求：
（1） 用你自己设计的一种线性收敛的迭代法求方程的根，然后用斯蒂芬森加速迭代计算。

输出迭代初值、各次迭代值及迭代次数。

（2）用牛顿迭代法求方程的根，输出迭代初值、各次迭代值及迭代次数，并与（1）的结果比较。

（3）用MATLAB 内部函数solve 直接求出方程的所有根，并与（1）、（2）的结果进行比较。

数值实验五
1． 观察直接利用拉格朗日插值多项式的病态性
实验目的：如果直接利用插值多项式的定义，即
i n
i n i i i y x a x a x a a x y =++++= 2210)(
确定系数，就必须解方程
⎪⎪⎪⎪⎪⎭⎫
⎝
⎛n n n n
n
x x
x x x x 111110
0⎪⎪⎪⎪
⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n n y y y a a a 1010
该方程的系数矩阵为范德蒙（Vandermonde ）矩阵，它是严重病态的，因此这种方法不能使用。

本实验的目的是观察这种病态性。

实验内容：取不同的n ，并在[0，1]区间上取n 个等分点，计算上述系数矩阵A 的条件数，画出㏑(cond(A)）随着n 变化的图形。

（可利用MATLAB 的内部函数cond （A ）求出矩阵的条件数，当n 取不同的值时，可分别用内部函数plot 画出所要求的图形）。

2． 龙格现象的发生、防止，插值效果的比较
实验目的：观察拉格朗日插值的龙格（Runge ）现象，探索避免此现象发生的方法，比较插值效果。

实验内容：对区间[-5，5]作10等分，并给出以下三个函数：
（1）
215x y +=
（2）x y arctan = （3）41x x
y +=
分别对上述函数计算点k x 上的值，作出插值函数的图形并与)(x f y =的图形比较。

实验要求：
（1） 编写程序对上述函数作拉格朗日插值，并在MATLAB 中用内部函数plot 利用插值点绘制函数的图形。

（2） 编写程序对上述函数作分段线性插值，并在MATLAB 中用内部函数plot 利用插值点绘制函数的图形。

（3） 编写程序对上述函数作三次样条插值，并在MATLAB 中用内部函数plot 利用插值点绘制函数的图形。

（4） 在MATLAB 中用内部函数ezplot 可直接绘制出以上函数的图形，并与作出的插值函数的图形进行比较。

3． 牛顿插值法
实验目的：学会牛顿插值法，并应用算法于实际问题。

实验内容：给定函数 x x f =
)(，已知：
414214.1)0.2(=f 449138.1)1.2(=f 483240.1)2.2(=f
516575.1)3.2(=f 549193.1)4.2(=f
实验要求：
（1）用牛顿插值法求4次Newton 插值多项式在2.15处的值，以此作为函数的近似值)15.2(15.2N ≈。

在MATLAB 中用内部函数ezplot 绘制出4次Newton 插值多项式的函数图形。

（2） 在MATLAB 中用内部函数ezplot 可直接绘制出以上函数的图形，并与作出的4次Newton 插值多项式的图形进行比较。

4． 拟合方式实验
实验目的：学会用最小二乘法求拟合数据的多项式，并应用算法于实际问题。

实验内容：给定数据点
),(i i y x 如下：
i x 0 0．5
0．6
0．7
0．8
0．9
1．0
i y
1
1．75 1．96 2．19 2．44 2．71 3．00
实验要求：
（1） 编写程序用最小二乘法求拟合数据的多项式，并求平方误差，作出离散函数
),(i i y x 和拟合函数的图形。

（2） 用MATLAB 的内部函数polyfit 求解上面最小二乘法曲线拟合多项式的系数及平方误差，并用MATLAB 的内部函数plot 作出其图形，并与（1）的结果进行比较。

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

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

实验内容：
（1）在[-1，1]区间上取n=20个等距节点，计算出以相应点上x e的值作为数据样本，
编写程序以1，
)
(
,
),
(
),
(
2
1
x
P
x
P
x
P
l
为基函数作出9,7,5,3
=
l次的最小二乘多项式，其
中，
)
(x
P
l是勒让德多项式。

画出
))
(
ln(A
cond随着n变化的曲线，其中A是确定最小二
乘多项式系数的矩阵。

计算出不同阶最小二乘多项式给出的最小偏差
∑=-
=
n
i
i
i
y x
y
l
1
2
)
)
(
(
)(
σ
（2）用MATLAB的内部函数polyfit求解上面最小二乘法曲线拟合多项式的系数，并用MATLAB的内部函数plot作出其图形。

数值实验六
1. Newton-cotes型求积公式
实验目的：学会Newton-cotes型求积公式,并应用该算法于实际问题.
实验内容：求定积分⎰π0cos xdx
e x
实验要求：选择等分份数n,用复化Simpson求积公式求上述定积分的误差不超过8
10-的近似值,用MATLAB中的内部函数int求此定积分的准确值,与利用复化Simpson求积公式计算的近似值进行比较。

2. Romberg算法
实验目的：学会数值求积的Romberg算法,并应用该算法于实际问题.
实验内容：求定积分
⎰15.0dx x
实验要求：
（1）要求程序不断加密对积分区间的等分,自动地控制Romberg算法中的加速收敛过程,直到定积分近似值的误差不超过6
10-为止,输出求得的定积分近似值。

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

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

（2）可用MATLAB 中的内部函数int 求得此定积分的准确值与Gauss 型求积公式求得的值进行比较。

数值实验七
1. 常微分方程性态和R-K 法的稳定性
实验目的：学习掌握4阶经典R-K 方法，体会参数和步长对问题的影响。

实验内容：求解常微分方程初值问题
1,01(0)1dy
y x x dx y αα⎧=-+<<⎪⎨⎪=⎩
其中5050α-≤≤。

实验要求：
（1）对参数α取不同的值，取步长0.01h =，编写程序用4阶经典R-K 法计算，将计算结果画图比较，并分析相应的初值问题的性态。

（2）取参数α为一个绝对值不大的负值和两个不同的步长，一个步长使参数,h α在经典R-K 法的稳定域内，另一步长在经典R-K 法的稳定域外，分别用R-K 法编写程序计算并比较计算结果，取全域等距的10个点上的计算值。

（3）用MATLAB 中的内部函数dsolve 求此常微分方程初值问题的解并与（1）与（2）中的结果进行比较。

2. 解初值问题各种方法比较
实验目的：掌握了解各种解初值问题的方法，体会步长对问题解的影响。

实验内容：给定初值问题
22,12(0)1x
dy y x e x dx x y ⎧=+<≤⎪⎨⎪=⎩
其精确解为
2()x
y x e e =- 实验要求： 分别按
（1）欧拉法，步长0.025,0.1h h ==； （2）改进的欧拉法，步长0.05,0.1h h ==； （3）四阶标准龙格-库塔法，步长0.1h =； 编写程序求在节点
10.1(1,2,10)k x k k =+= 处的数值解及误差，并比较各方法的
优缺点。

用MATLAB 中的内部函数dsolve 求此常微分方程初值问题的解并与上述结果进行比较。

3. 认识刚性微分方程
实验目的：了解认识刚性微分方程。

实验内容：⎪
⎪
⎪
⎩
⎪
⎪
⎪
⎨
⎧
-
=
=
-
=
+
+
-
=
2
)0(
,0
)0(
25
.
1000
75
.
999
2000
v
u
v
u
dt
dv
v
u
dt
du
实验要求：
任选一显式方法，编写程序取不同步长求解上述方程。

第四部分MATLAB入门
一、MATLAB简介
MATLAB全名叫Matrix Laboratory，是矩阵实验室的意思。

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

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

MATLAB是适用于科学和工程计算的数学软件系统。

1．用户界面介绍
MATLAB用户界面如图1所示，包括主菜单、主工具栏、命令窗口、命令历史窗口、工作间管理窗口、当前路径窗口、编译窗口、图形窗口、启动按钮、帮助浏览器等。

下
面简单介绍命令窗口，其他的内容在使用的过程中会逐渐熟悉。

命令窗口是用于输入数据，运行MA TLAB函数和脚本并显示结果的主要工具之一，命令窗口如图1右边中下区域所示。

图1 MA TLAB用户界面
2．常量
常量是MA TLAB语言预定义的一些变量，在默认的情况下这些变量的值为常数。

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

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

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

注：
（1）使用变量不要求事先声明；
（2）也存在变量的作用域问题，一般视为局部变量，仅在其调用的M文件内有效。

若要定义全局变量，应加上关键字global，并一般用大写；
（3）如果在对某个变量赋值时，如果该变量已经存在，系统则会自动使用新值来代替该变量的旧值。

例如在命令窗口中输入“a=1;a=2”命令，则会得到
如下结果：
a =
2
4．算术运算
MATLAB中用“+”、“-”、“*”“/”和“^”分别表示算术运算中的加、减、乘、除和乘方。

例如：计算
3
3
1
2
1
62
5
)
3
1
(
256÷
+
⨯-
>>256^(1/6)*(1/3)^(-1/2)+5^(1/3)/2^3
ans =
4.5782
5．数字的输入、运算
（1）对于简单的数字运算，可以直接在命令窗口以惯用的形式输入。

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

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

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

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

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

（7）数据存储和运算都以双精度进行的。

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结束。

(2) while循环语句
与计算机的其他高级语言一样，while循环语句要由关系运算和逻辑运算给出的逻辑控制，该语句的一般形式为
while（逻辑表达式）
（一族可执行语句）
end
MATLAB中的关系运算有
＝＝相等～＝不等，
<＝小于等于，< 小于，
>＝大于等于，> 大于.
而MA TLAB提供的逻辑运算有
& 与，∣或，～非。

(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结尾的。

(4) 内部函数
MATLAB提供了丰富的运算函数，只需正确调用其形式就可得到满意的结果，常用的运算函数如表1所示。

表1 MA TLAB常用函数表
函数名函数功能函数名函数功能函数名函数功能
sin 正弦acosh 反双曲余弦diff 求导
cos 余弦atanh 反双曲正切int 积分
tan 正切acoth 反双曲余切dsolve 求解微分方程
cot 余切asech 反双曲正割dot 向量点乘
sec 正割acsch 反双曲余割cross 向量叉乘
csc 余割pow2 以2为底的幂函数fix 向零方向舍入
asin 反正弦exp 以e为底的幂函数floor 向负方向舍入
acos 反余弦log 自然对数ceil 向正方向舍入
atan 反正切 log10 以10为底的对数 round 四舍五入 acot 反余切 log2 以2底的对数 mod 有符号求余 asec 反正割 sqrt 平方根 rem 无符号求余 acsc 反余割 lim 求极限 sign 符号函数 asinh
反双曲正弦
fsolve
求非线性方程的根
abs
绝对值
(5) 用户自定义函数
MATLAB 允许用户使用M 文件定义函数，但必须符合一定的规则。

例如下面一段程序存为funsim.m
function p=funsim(x) % define a simple function p=sqrt(x)-2*x^3+cos(x);
所给出的是函数)cos(23
x x x +-。

函数的输出可以多于一个， 例如
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 点和线对应符号
线符号点符号
实线—实心圆点.
虚线- - 加号+
点: 星号*
虚线间点-. 空心圆点。

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

目前大部分的硬盘拷贝设备尚不能支持颜色，但彩色的显示已经是十分普及的，所以彩色显示在计算可视化中是十分有益的。

颜色在绘图中的使用方法与线形的控制方法一样，常用的颜色有
蓝色“b”,黄色“y”，红色“r”，绿色“g”
例如plot(x,y,,’—r’) 将x和y对应的元素定义的点依次用红色实线联接。

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

title(′···′)；在图形的上方显示′′中所指定的内容；
xlabel(′···′)；将′′中所指定的内容标在x轴；
ylabel(′···′)；将′′中所指定的内容标在y轴；。