数值计算方法I试验手册

数值计算方法I试验手册
计算方法实验手册
江西财经大学信息管理学院数学与决策科学系
2003.1
目录
预备实验MATLAB使用练习 (2)
实验一数值计算中误差的传播规律 (6)
实验二数值计算中的算法稳定性 (8)
实验三Gauss消去法主元素的选取与算法的稳定性 (10)
实验四观察Runge现象和对非光滑函数进行插值的可能性 (13) 实验五观察及改善最小二乘拟合的数值不稳定现象 (15)
实验六曲线逼近方法的比较 (17)
实验七利用数值微分方法求解偏微分方程 (18)
实验八比较不同数值积分公式对积分方程解的影响 (20)
实验九高维积分数值计算的Monte-Carlo方法 (23)
实验十观察前进Euler公式的收敛性和数值不稳定性 (26)
实验十一利用实际例子筛选有效算法 (28)
实验十二估计线性方程组的性态与条件数 (29)
实验十三迭代法的算法设计与比较 (31)
实验十四迭代法收敛阶的估计 (33)
实验十五迭代法过程中的分支与混沌 (34)
综合实验一 (35)
综合实验二 (36)
预备实验MATLAB使用练习
一、实验目的
MATLAB是本实验课的主要软件，尽快熟悉其基本命令和语法，以适应后续实验对MATLAB的使用要求．
１．掌握MATLAB中常用的基本命令和函数；
２．掌握MATLAB的简单作图方法；
３．掌握MATLAB的简单编程技巧．
二、实验内容
１.矩阵和向量的输入和运算
在MATLAB的Commond窗口中输入以下命令，并仔细观察执行后的结果．
（1）矩阵的输入：
A=[1 2 3;4 5 6]?（分号表示矩阵换行，?表示回车）
A=[1,2,3;4,5,6]?（空格或逗号表示同一行元素的分隔）
A=[1 2;3 4];?（在一行的末尾加分号，运行结果将不显示）
Z=zeros(2,3)?（3
2?阶零矩阵）
E=ones(2)?（2
2?阶全1矩阵）
I=eye(2,4)?（4
2?阶对角线元素为1的矩阵）
R=rand(2,3)?（3
2?阶(0,1)均匀分布随机矩阵）
H=hilb(4)?（4阶hilbert矩阵）
另外还有标准正态分布矩阵randn(m,n),n阶幻方矩阵magic(n)等.
（2）矩阵元素的调用和赋值：
a=A(1,3)?（矩阵A的第１行第３列元素赋值给a）
A(1,3)?（查看矩阵A的第１行第３列元素）
A(1,3)=0?（将A的第１行第３列元素赋值为０）
A(3,4)=1?（注意矩阵A原本是3
2?阶矩阵）
B(2,3)=2?（注意前面没有给出矩阵B）
（3）矩阵的裁剪和拼接：
A?（显示矩阵A）
A(2,:)?（A的第２行）
A(:,3)?（A的第3列）
C1=A(1:2,:)?（将A的第１,２行赋值给矩阵C1）
C２=A(:,２:４)?（将A的第２～４列赋值给矩阵C2）
C３=A(２:３,２:４)?（将A的第２,３行,第２～４列赋值给C3）C4= A(1:２:３,4:-1:2)?（将A的第1、３行,第４、３、２列赋值给C3）
A?（显示矩阵A为4
3?）
B1=[A,C2]?（左右拼接矩阵A和C2为B1）
B2=[A；C1]?（上下拼接矩阵A和C1为B2）
（４）矩阵的运算
MATLAB为矩阵的运算提供了下列矩阵运算符：
＋（加法）；－（减法）；＇（转置）；＊（乘法）；＾（乘幂）；＼（左除）；/（右除）.＊（点乘法）；.＾（点乘幂）；.＼（点左除）；./（点右除）
矩阵的“＋”、“－”、“＇”、“＊”、“＾”必须遵循矩阵的运算规则，点运算“.＊”、“.＾”、“.＼”、“./”是指矩阵对应元素之间进行相应的运算．
在MATLAB的Commond窗口中继续输入以下命令，并仔细观察结果：
C=zeros(3)?
D=ones(3)?
比较C+1和C+D．
A=[1 2 ２;4 5 6;7 8 9]?
B=[1 2 3]?（B为一列矩阵即列向量）
X=A\B＇?（左除表示左边乘以A的逆矩阵,X为方程组AX=B’的解）X=B/A?（右除表示左边乘以A的逆矩阵,X为矩阵方程XA=B的解）C=A?
A*C?
A .*C?（A,C的对应分量相乘）
A^2?
A .^2?（A的每个分量平方）
A/A?（相当与1-
AA）
A./A?（结果为全１矩阵）
在MATLAB中，数组、向量与矩阵的形式完全一致，并且运算与矩阵也基本类似，但要注意行列向量之分．
除了与矩阵一样的输入方式外，数组（向量）还有以下的输入方式：
a=1:5?
b=1:2:8?
c=10:-2:0?
b=[0:2:8,zeros(1,3)]? (相当与矩阵的拼接)
c=linspace(1,10,10)? (从１到10的共10个数值的等差数组) x=0:pi/4:pi ? (pi 是圆周率) ２.变量和函数
MATLAB 中的变量区分大小写字母，第一个必须是字母，不需要说明变量类型或维数，系统自动分配内存空间,字符串是用单引号括起来的字符集合，可以像向量或矩阵一样进行拼接与剪裁．
MATLAB 有几个特殊的变量：pi(圆周率)；eps （最小浮点数）；inf(正无穷大，特指
1)；
NaN(不定值，特指0
0)；i,j(虚数单位)．观察下列结果：
a=[0 1 0],b=[1 0 0],c=a./b ?
MATLAB 提供了大量的函数，按照分类分为标量函数、向量函数和矩阵函数三种类型．请用help 命令查看下列各种函数的调用格式：标量函数：sin cos tan cot sec csc asin acos atan acot asec acsc
sqrt exp log log10 abs round floor ceil fix sign real imag 向量函数：max min sum length mean median prod sort
矩阵函数：zeros ones eye rand randn diag det rank inv eig trace expm poly norm cond lu qr svd ３.基本图象功能
二维图象最常用的命令有：
基本作图命令：plot fplot polar
图象注释命令：title xlabel ylabel text gtext legend 坐标管理命令：axis
窗口管理命令：hold on hold off grid subplot figure 二维图象最常用的线型和颜色有：
线型（线方式）：－实线：点线－．虚点线――波折线线型（点方式）：．圆点＋加号＊星号 xX 形 o 小圆
颜色： y 黄 r 红 g 绿 b 蓝 w 白 k 黑 m 紫 c 青４.程序设计
关系运算符：<小于>大于<=小于或等于>=大于或等于==等于 ~=不等于
关系运算比较两个数值，当指出的关系成立时结果为１，否则为０，关系运算可以作用于两个同型的矩阵或向量，比较时逐个元素或分量比较，结果是一由０和１组成的矩阵或向量．
逻辑运算符：&与运算（表示且）|或运算（表示或）~非运算（表示不）
流控制语句：if（条件语句）；while（循环语句）；for（循环语句）；switch（条件语句），它们都用end结束，具体形式用help命令查看．
文本M－文件：由语句组成的MATLAB命令集，运行时不必输入参数．
函数M－文件：可以方便自己调用的函数文件，第一行有特殊要求，其形式必须为function <变量名>＝＜函数名＞（＜自变量＞）该文件的存储名必须与函数名相同．
三、实验要求
本实验主要是为了熟悉MA TLAB而设计的．MATLAB是本实验课的主要软件，它将计算、可视化和编程功能集成在非常便于使用的环境中，是一个交互式的、以矩阵计算为基础的科学和工程软件，具
有编程效率高、计算功能强、使用简便和易于扩充等特点．要求通过本实验，尽快熟悉其基本命令和语法，以适应后续实验对MA TLAB的使用要求．
实验一数值计算中误差的传播规律
一、相关原理
1．设*x 是x 的一个近似值，则有如下概念绝对误差x x x E -=**)(，简称误差，简记为E 误差限*)(*)(x E x ≥ε，简记为ε 相对误差 x x x x x E x E r -==
**)(*)(，简记为r E
或 *
**
*)(*)(*x x x x x E x E r -=
=
，简记为*r E
相对误差限*)(*)(x E x r r ≥ε，简记为r ε
或*)(**)(*x E x r r ≥ε，简记为*r ε
2．设),(21x x f y =，*1x 是1x 的一个近似值，*
2x 是2x 的一个近似值，),(**
2*
1x x f y =，
则*y 为相应的),(21x x f y =的近似值．设
21,E E 分别为*
2*
1,x x 的绝对误差；
21,εε分别为*
2*1,x x 的绝对误差限；
*
2*
1,r r E E 分别为*
2*
1,x x 的相对误差；
*2*1
,r r ε
ε
分别为*
2*1,x x 的相对误差限；则 *y 的绝对误差为 2* 2
1*
1*)(E x f E x f
y E
+
≈ *y 的绝对误差限为 2* 2
1*1εεε
+
=x f x f
*y 的相对误差为 *
2*
2
*
2*
1*
1
*
1***)(r r r E x f
y x E x f
y x y E
≈
*y 的相对误差限为 *
2*
2*
2*
1
*
1*
1*
**r r r
x f
y x x f
y x εεε
+
=
3. 设k
m a a a x 10.0*21?±= （01≠a ）是x 的一个近似值，如果*x 有n （m n ≤）
位有效数字，则
n
r a x E -?≤
11
10
21*)(
如果*x 的相对误差满足
n
r a x E -?+≤
)
1(21*)(
则有*x 至少有n （m n ≤）位有效数字．
二、实验目的
１．观察并初步分析数值计算中误差的传播；２．观察有效数字与误差传播的关系．
三、实验内容
１．使用MATLAB 的help 命令学习MATLAB 命令digtis 和vpa 的用途和使用格式；２．在4位浮点数下解二次方程01622 =++x x ；３．计算下列５个函数在点2=
x 处的近似值
（１）6
0)1(-=x y ，
（２）6
1)
1(1+=
x y ，
（３）3
2)23(x y -=，（４）3
3)
23(1x y +=
，
（５）x y 70994-=．
四、实验要求
本次实验包含三个相对独立的内容．
１．在内容１中，请解释两个命令的格式和作用；
２．求解方程01622
=++x x 时，分别使用求根公式和韦达定理两种方法，并比较其有效数字和相对误差；
３．实验内容３中的５个函数在2=
x 处的精确值都是相等的，若取
4.12≈进行计
算，计算各函数的结果，作图观察并比较它们的绝对误差（作图区间可取]42.1,4.1[甚至更小）,并从算法设计原则上说明原因．实验二数值计算中的算法稳定性
一、实验目的
误差扩张的算法是不稳定的，是我们所不期望的；误差衰竭的算法是稳定的，是我们努力寻求的，这也是贯穿本课程的目标．本实验的目的是通过实验，体会稳定性在选择算法中的地位．
二、实验内容
考虑一个简单的由积分定义的序列 ,2,1,0,
1
1
==
-n dx e
x I x n n （E １）
显然n I n ?>,0．当0=n 时，显然有 1
1
1
01---==
e
dx e
I x
而对1≥n ，由分部积分法不难得到 ?
-----
==1
1
1
10
1
1
1
dx e
nx
e
x dx e
x I x n x n x n n
即
,3,2,1,
11
=-=-n nI
I n n （E ２）
同样由（E １）式，我们还可以得到
1
11
1
1
+=
≤
=
-n dx x dx e
x I n
x n n （E ３）
由递推关系（E ２）式，可以得到计算积分序列（E １）的两种算
法：１．第一种是（E ２）的直接应用，即
=-=-=--
,3,2,1,111
1
0n nE
E e E n n （E ４）
２．另一种是利用（E ２）和（E ３）的变形得到
--=-==-1
,2,3,,2,1,101
N N n n I I I n
n N （E ５）
三、实验要求
１．分别利用算法（E ４）和（E ５）式计算，并在计算中分别采用５位、６位和７位有效数字，请判断哪种算法能给出更精确的结果；
２．两种算法的优劣，与你的第一感觉是否吻合？请从理论上证明你实验得出的结果，解释实验的结果．设算法（E ４）中0I 的计算误差为0e ，由0I 递推计算n I 的误差为n e ；设算法（E ５）中N I 的计算误差为N ε，由N I 向前递推计算n I （N n <）的误差为n ε；如果在上面两种算法中都假定后面的计算不会再引入其他误差，试给出n e 与0e 、n ε与N ε的关系；
３．算法（E ４）中0e 通常会很小，当n 增大时，误差n e 的变化趋势如何？算法（E ５）中N ε通常相对比较大，当n 减小时，误差n ε又是如何传播的？即比较上述两个算法，当某一步产生误差后，该误差对后面的影响是衰减还是扩张？
４．通过理论分析和数值实验，针对算法（E ４）和（E ５）的稳定性，给出你的结论．
实验三 Gauss 消去法主元素的
选取与算法的稳定性
一、相关原理
解线性方程组
n
n
n R b R A b Ax ∈∈=?,
,
其中
=
nn n n n n a a a a a a a a a A
2
1
22221
11211，
=n b b b b 21，记A A =)1(，b b
=)
1(．Gauss 列主元消去法通常包含下面四个步骤：
１．选取主元素：在矩阵)
(k A
的第k 列元素)
(k kk a 及其以下的元素中按某种规则（如按模
最大或最小），寻找到列主元，记录其行号，记为k i ，并且当
0)(,≠k k i k
a 时进行换行，否则退出消去法．
２．换行：交换矩阵),()
()
(k k b
A
的第k 行与k i 行，即取
行—行—k k
i
i k I k
=11
1
1
1
1
1
1
,
)(k i k ≥ 作乘积),()
()
(,k k k
i
b
A
I k
，即为交换矩阵),()
()
(k k b
A
的第k 行与k i 行，当k i k =时，表示不需
换行．换行后，仍将),() ()
(,k k k
i
b
A
I k
的元素记为)
(k ij
a ．
３．消元过程：定义行乘数) ()(k kk
k ik
ik a a m =
（n k i ,,1 +=），令
---=++11
1
11,,2,1k
n k
k k k k
m m m L
作乘积),()
()
(,k k k
i
k b
A
I L k
),()
1()
1(++=k k b
A ．
重复以上过程，1,,2,1-=n k ，当1-=n k 时，Gauss 选主元的消元过程结束．４．回代：通过上述过程，增广矩阵),()
1()
1(b
A →),()
()
(n n b
A
，并且)
(n A
是上三角矩
阵
=)()
2(2)
2(22
)1(1)
1(12
)
1(11)
(n nn n n n a a a a a a A
，
=)()2(2)
1(1)
(n n n b b b b 当0)
(≠n nn a 时，利用下面的回代公式即可得到方程组的解：
--=-==∑+=)1,2,,2,1(,)
(1
)
()()
()
( n n i a a b x a b x i ii n
i j i ij
i i i n nn
n n n Gauss 列主元消去法结束．
二、实验目的
１．简单了解条件数的概念（条件数的理论将在教材的第六章介绍）；２．通过实验了解算法的稳定性与主元素及条件数的关系．
三、实验内容
Gauss 消去法是我们在线性代数中已经熟悉的．但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元素的选择．主元素的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题．
考虑线性方程组
n
n
n R b R
A b Ax ∈∈=?,
编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss 消去过程．
三、实验要求
１．取矩阵
=68
168
168
16
A ，???????? ??=1415157 b ，则方程组有解精确解
=1111 x ，取10=n ，用MATLAB 的cond 命令计算矩阵A 的条件数．让程序自动选取主元，结果如何？
２．现选择程序中手动选取主元素的功能．每步消去过程总选取按模最小或按模尽可能小的元素作为主元素，观察并记录结果．若每步消去过程总选取按模最大的元素作为主元素，结果又如何？分析实验的结果．
３．取矩阵阶数20=n 或者更大，重复上述实验过程，观察并记录分析不同的问题及消去过程中选择不同的主元素时计算结果的差异，说明主元素的选取在消去过程中的作用．
４．选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数．重复上述实验，观察记录并分析实验的结果．
５．选择著名的Hilbert 矩阵H 作方程组b Hx =的系数矩阵，其中
n
n ij h H ?=)(，1
1-+=
j i h ij ，n j i ,,2,1, =
n b b b b ),,,(21 =，∑
==
n
j ij i h b 1
，n i ,,2,1 =
该方程组的精确解为T
x )1,,1,1( =．分别选择方程的维数为3，6，10，计算矩阵H 的条件数，并用Gauss 列主元消去法求解，记录结果，并给出你的结论．实验四观察Runge 现象和对非光滑函数进行插值的可能性
一、相关原理
)(x f 在节点b x x x a n ≤<<<≤ 10处的函数值为n y y y ,,,10 ，构造其Lagrange 插值多项式)(x L n 的插值基函数为
∏
≠=--=
n
j
i i i
j
i j x x
x x x l 0
)(，n j ,,2,1,0 =
Lagrange 插值多项式为
∑
==
n
j j j n x l y x L 0
)()(
其截断误差为
)()!
()(1)
1(x n f
x R n n n +++=
ωξ
其中),(b a ∈ξ，∏=+-=n
i i n x x x 0
1)()(ω
二、实验目的
１．观察高次Lagrange 插值多项式)(x L n 的Runge 现象；２．观察非光滑函数进行多项式插值的可能性．
三、实验内容
１．考虑在一个固定的区间上用Lagrange 插值逼近一个函数．显然Lagrange 插值中使
用的节点越多，插值多项式似的次数就越高．我们自然关心插值多项式的次数增加时，
)(x L n 是否也更加靠近被逼近的函数．设区间]1,1[-上函数
2
2511)(x
x f +=
考虑区间]1,1[-上的一个等距分割，节点为
n i n
i x i ,,2,1,0,21 =+
-=
作)(x f 在]1,1[-上的Lagrange 插值多项式
∑
=+=
n
i i i
n x l x x L 0
2
)(2511)(
其中)(x l i 为Lagrange 插值基函数．
２．连续非光滑函数的几何特性非常差，在几何图象上一般会出现大量的尖点．在构造非光滑函数的多项式插值时，由于多项式具有高阶光滑度，两者之间会产生怎样的现象？选择区间]1,0[上的连续非光滑函数
x k x g πsin )(=
作)(x g 区间]1,0[上的Lagrange 插值多项式．
四、实验要求
１．选择不断增大的节点数目 ,3,2=n ，画出原函数)(x f 及插值多项式)(x L n 在区间的]1,1[-上的图象，比较并分析实验结果．２．选择其他的函数，例如定义在区间]5,5[-上的函数
)arctan()(,1)(4
x x v x
x x u =+=
重复上述的实验过程，观察其结果又将如何．
３．如果不取等距节点，而改为取如下节点
++-+-=
π)
1(21
2cos 2
2
n i a b a b x i ，n i ,,2,1,0 = 以n x x x ,,,10 为插值节点构造上述函数)(),(),(x v x u x f 的Lagrange 插值多项式，比较其结果．４．选择不同的k 和n ，用等距节点作)(x g 的n 次Lagrange 插值多项式，观察其误差大小及收敛情况．
实验五观察及改善最小二乘拟合的
数值不稳定现象
一、相关原理
１．最小二乘拟合：设)(x f 是定义在区间],[b a 上的函数，m x x x ,,,10 是区间],[b a 上的一组节点，)(i i x f y =（m i ,,2,1,0 =）是函数)(x f 对应于节点的函数值．)(,),(),(10x x x n 是定义在区间],[b a 上的线性无关的连续函数，Φ是由拟合基函数)(,),(),(10x x x n 生成的函数类，即
)}(,),(),({10x x x span n =Φ
在Φ中必存在一个函数∑
==n
j j j x a x S 0
*
)()(*?是函数)(x f 在],[b a 上的最小二乘拟合
函数，其中拟合系数*j a 满足法方程组
=
),(),(),(),()
,()
,(),()
,(),(),(),(),(101
1021110101000f f f a a a n n n n n n n ，记为b Aa =，其中
=
,()
,(),()
,(),(),(),(),(1021110101000n n n n n A ，??????? ??=n a a a a 10，
=),(),(),(1
0f f f b n ∑==
m
k k j k i
j i x x 0
)()(),(??
，∑==
m
k i k i
i y x f 0
)(),(?
最小二乘解的平方误差为
∑
=-=
m
i i i y x S 0
2
22
)
)(*(δ
)*,(),(f S f f -=
)*,(),(b a f f -=
∑
==-
=
n
j j j m
i i f a y 0
*
2
),(?
２．Legendre （勒让德）正交多项式在区间]1,1[-上正交多项式
1)(0=x P
,4,3,2,1,])1[(!21
)(2
=-=
n x
dx
d
n x P n
n
n n
n
称为Legendre 正交多项式．Legendre 正交多项式具有递推公式：x x P x P ==)(,1)(10
)()()12()()1(11x nP x xP n x P n n n n -+-+=+， ,2,1=n
二、实验目的
1．观察最小二乘拟合多项式的数值不稳定现象；2．探索改善最小二乘拟合多项式的不稳定现象的可能性．
三、实验内容
１．在区间]1,1[-上取20=m 个等距节点，计算出相应节点上函数x
e 的值做为样本数据，以n
x x x ,,,,12
为基函数作出9,7,5,3=n 次的最小二乘多项式；
２．在区间]1,1[-上取20=m 个等距节点，计算出相应节点上函数x
e 的值做为样本数据，以Legendre 多项式)(,),(),(10x P x P x P n 为基函数作出9,7,5,3=n 次的最小二乘多项式．
四、实验要求
１．在区间]1,1[-上取20=m 个等距节点，计算出相应节点上函数x
e 的值做为样本数据，以n
x x x ,,,,12
为基函数作出9,7,5,3=n 次的最小二乘多项式．用MATLAB 的＂cond ＂命令，求出相应的法方程组的系数矩阵A 的条件数，画出n A cond ~))(ln(之间的曲线．计算出不同次数的最小二乘拟合多项式的平方误差22
)
(n δ．
２．仍在区间]1,1[-上取20=m 个等距节点，计算出相应节点上函数x
e 的值做为样本数据，以Legendre 多项式)(,),(),(10x P x P x P n 为基函数作出9,7,5,3=n 次的最小二乘多项式．用MATLAB 的＂cond ＂命令求出相应的法方程组的系数矩阵A 的条件数，画出n A cond ~))(ln(之间的曲线．计算出不同次数的最小二乘拟合多项式的平方误差22
)
(n δ．
３．比较使用两类不同基函数得到的拟合结果，给出你的结论．实验六曲线逼近方法的比较
一、实验目的
掌握各种函数逼近方法的适用范围及实际应用中方法选择时应注意的问题．
二、实验内容
曲线的拟合和插值，是逼近函数的基本方法，每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的．现在仍然考虑实验四中的著名问题，即用不同的插值和拟合方法作区间]1,1[-上函数 2
2511)(x
x f +=
的逼近函数，并加以比较．下面的MATLAB 程序给出了该函数的二次和三次拟合多项式．
x=-1:0.2:1; y=1./(1+25*x.*x); xx=-1:0.02:1; p2=polyfit(x,y,2); yy2=polyval(p2,xx); plot(x,y,’0’,xx,yy2) xlabel(‘x’); ylabel(‘y’); hold on
p3=polyfit(x,y,3); yy3=polyval(p3,xx); plot(xx,yy3) hold off
适当修改上述MATLAB 程序，也可以拟合其他你感兴趣的函数．
三、实验要求
１．将拟合的结果与Lagrange 插值及样条插值的结果比较；
２．归纳总结数值实验的结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题．。