计算方法实验指导书1

第一章 绪论
一、主要要求
通过实验，认真理解和体会数值计算的稳定性、精确性与步长的关系。

二、主要结果回顾：
1、算法：电子计算机实质上只会做加、减、乘、除等算术运算和一些逻辑运算，由这些基本运算及运算顺序规定构成的解题步骤，称为算法．它可以用框图、算法语言、数学语言或自然语言来描述。

用计算机算法语言描述的算法称为计算机程序。

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

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

3、算法的稳定性：一个算法如果输入数据有误差，而在计算过程中舍入误差不增长，则称此算法是数值稳定的，否则称此算法为不稳定的。

换句话说：若误差传播是可控制的,则称此算法是数值稳定的,否则称此算法为不稳定的。

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

三、数值计算实验（以下实验都需利用Matlab 软件来完成） 实验1.1（体会数值计算精度与步长关系的实验）
实验目的：数值计算中误差是不可避免的，要求通过本实验初步认识数值分析中两个重
要概念：截断误差和舍入误差，并认真体会误差对计算结果的影响。

问题提出：设一元函数f :R →R ，则f 在x 0的导数定义为：
h
x f h x f x f h )
()(lim
)('000
0-+=→
实验内容:根据不同的步长可设计两种算法,计算f 在x 0处的导数。

计算一阶导数的算法有两种：
h x f h x f x f )
()()('000-+≈ （1）
h
h x f h x f x f 2)
()()('000--+≈
（2）
请给出几个计算高阶导数的近似算法，并完成如下工作：
1、对同样的h ，比较(1)式和(2)式的计算结果；
2、针对计算高阶导数的算法，比较h 取不同值时(1)式和(2)式的计算结果。

实验要求：选择有代表性的函数f (x )（最好多选择几个），利用Matlab 提供的绘图工具画出该函数在某区间的导数曲线f (s)(x )，再将数值计算的结果用Matlab 画出来，认真思考实验的结果。

实验分析：不论采用怎样的算法，计算结果通常都会有误差。

比如算法(1)，由Taylor 公式，知：
)()(''!
2)(')()(302
000h O x f h x hf x f h x f +++=+ 所以有
)()(''!
2)(')()(20000h O x f h
x f h x f h x f ++=-+
利用上式来计算f ’(x 0)，其截断误差为： )()(''!
2201h O x f h
T += 所以误差是存在的，并且当步长h 越来越小时，(1)的近似程度也越来越好。

类似地可以分析(2)的截断误差为：)()(!
330)
3(22h O x f h T += 上述截断误差的分析表明(2)是比(1)更好的算法，因为对步长h (<<1)，(2)比(1)更接近于f ’(x 0)。

计算方法的截断误差是数值计算中误差的重要来源，但不是唯一的！！如果在实验中确已将h 取到足够小的话，特别在高阶导数的计算中，就会发现当h 小到一定程度后，计算结果的误差不但不再减少，反而会变大！(参见图1)
事实上当步长h 过小时，计算结果的误差变大就是由于舍入误差的缘故。

截断误差用我们原有的数学思维方式就比较容易理解的，而舍入误差则是本课程引入的一个新概念。

要真正理解舍入误差，特别是它在计算中的传播及最终对计算结果的影响，是初步具备科学计算能力的重要标志。

希望大家在完成实验后，认真仔细去体会截断误差和舍入误差的含义及对计算结果的影响。

（图1）
实验1.2（理解误差传播与算法稳定性实验）
实验目的：体会算法的稳定性在选择算法中的地位。

问题提出：考虑一个简单的积分序列 ⎰
-=1
1dx e x I x n n n=1,2,….
显然I n >0,n=1,2,…当n=1时，得：e
dx xe I x 11
1
1==⎰- 当n ≥2时，由分部积分可得：11
11---==⎰
n x n n nI dx e x I n=2,3,……
另外，还有：1
11
1
1
+=
≤=
⎰⎰
-n dx x dx e x I n x n n 实验内容：由递推关系I n =1-nI n-1，可得计算积分序列{I n }的两种算法：
算法一、直接使用递推公式得：e
I 11=
I n =1-nI n-1 n=2,3…
算法二、利用递推公式变形得：
,...
3,210
1=-=
=-n n
I I I n
n N
10-7
10-5
10
-5
10-4 10-3 10-2 10-1 10-3
10-1
实验要求：用上述两种算法分别在计算中采用5位、6位和7位有效数字，请判断哪种算法给出的结果更精确。

实验分析：两种算法的优劣可能与你的第一感觉完全不同。

设算法一中I 1的误差为e 1，由I 1递推计算I n 的误差为e n 算法二中I N 的误差为εN ，由I N 向前递推计算I n (n<N)的误差为εn 如果假定上述两种算法在后面的计算中都不再引入其他误差，则有：
,...,3,2,!1==n e n e n
1,2,3,...,2,1,!
/!--==
N N m m N N
m εε
由此可见，算法一中的e 1可能很小，但在计算中它的影响急剧扩散，结果真实的数据很快被淹没了；而算法二中的εN 可能相对比较大，但在计算中误差影响不扩散，某一步产生误差后，该误差对后面的影响不断衰减。

注：误差扩散的算法是不稳定的，是我们所不期望的；误差衰减的算法是稳定的，是我们努力寻求的，也是贯穿本课程始终的目标。

第二章 一元非线性方程的解法
一、主要要求
编写二分法、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 ]区间二等分的中点2
0b
a x +=
，若f(x 0)=0,则x 0是f (x )=0的实根； 若f (a )f (x 0)<0 成立,则 x * 必在区间(a , x 0)内,取a 1=a ，b 1= x 0；否则 x *必在区间(x 0,b)内，取a 1= x 0，b 1=b , 这样,得到新区间(a 1,b 1),其长度为[a ,b ]的一半, 如此继续下去，进行k 次等分后,得到一系列有根区间：],[],[],[11k k b a b a b a ⊃⊃⊃ ，其中每一个区间长度都是前一个区间长度的一半，从而[a k ,b k ]的长度为k k k a
b a b 2
-=
-如此继续下去,则有这些区间将收敛于一点，该点即为所求的根.当做到第k 步时，有:
ε<-=-≤-+12
2|*|k k k k a
b a b x x （ε为给定的精度）
此时2
k
k k b a x +=
即为所求方程的近似解。

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）； 4）若b -a <ε，则输出满足精度要求的跟c ，否则转向2）
3、二分法程序框图（见图2）
4、不动点迭代法：
不动点迭代法又称简单迭代法，基本思想是构造不动点方程，以求得近似根。

即由方程f (x )=0变换为等价方程 x=ϕ(x ) ，这样原方程的根必满足：
x *= ϕ(x *) ，即ϕ(x ) 作用在x *上，其值不发生变化，因此我们也称x *为ϕ(x ) 的不动点，要求方程f (x )=0的根就转化为求ϕ(x ) 的不动点了。

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

如果{x k }收敛于x *，而ϕ(x )是连续函数时，那么x *即是ϕ(x )的不动点。

也即x *就是方程的根。

5、Newton 迭代法的思想：
给定方程f (x )=0，迭代公式为：)
(')
(1k k k k x f x f x x -=+，应用该公式来解方程的方法
就称为牛顿迭代法。

6、Newton 迭代法的算法：
1）给出初始近似根x 0及精度ε； 2）计算)
(')
(000x f x f x -
→x 1；
3）若|x 1-x 0|<ε或迭代次数达到预定指定的次数N ，则转向 4）；否则x 1→x 0，转向2）； 4）输出满足精度的根x 1，结束。

（程序框图略）
7、收敛性判定定理：
定理1：假定函数ϕ (x ) 满足下列条件： 1、对任意x ∈[a,b ] 有 ϕ(x ) ∈[a,b] 2、存在正数 L<1，使对任意x 1，x 2∈[a,b ] 有
|ϕ(x 1) - ϕ(x 2)|≤L|x 1-x 2| (0≤L<1)
则迭代过程 x k +1= ϕ (x k )对于任意初值x 0∈[a,b ]均收敛于方程x= ϕ(x)的根x* ，且有如下的误差估计式：
01*
1x x L
L x x k
k --≤-
局部收敛性定义： 如果存在不动点x *的某个邻域U(x *,δ)，使得对于任意初值x 0∈ U(x *,δ)，迭代公式x k +1= ϕ(x k ) （k =0,1,2......)产生的数列{xk }均收敛与x *，则称迭代公式x k +1= ϕ(x k )是局部收敛的。

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

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

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

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

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

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

实验2.1（求非线性方程根实验（二）） 实验目的：会求非线性方程的根 实验内容：试有方程：032
=-x e x
1、分别用Newton 迭代法和快速弦截法求上述方程的根；
2、比较两种迭代法的收敛速度。

四、具体操作见下例。

（以下实验都需利用Matlab 软件来完成） 例1：用二分法求方程x 2-x -1=0的正根，要求误差小于0.05 解：1）首先根据程序框图编制二分法的函数文件： erfen.m
2）再编写函数文件f c.m
funtion y=f c(x ) y=x ^2-x -1;
3）最后在Matlab 命令窗口输入调用命令：erfen (‘fc ’,0,10,0.05) 4）输出结果为 n = 56
ans = 1.6180
例2：求解非线性方程5x 2sinx-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)
即可求出 求方程fc 1=0在x0附近的近似解 得出结果
ans =
0.5018 3.1407 6.2832 9.4248
图3
注：用Newton 迭代法可求出多个根。

具体做法：先用fplot 命令作出函数的图形，再根据图形给出不同的初始值x 0，最后通过调用Newton 迭代法程序求出非线性方程的所有根。

五、思考题
思考题1、 在二分法中取区间中点，每次计算一次函数值。

如果每次把区间三等份，计算两个函数值，是否可以找出方程的根？给出它的算法。

和二分法比较它的效率如何？
思考题2、目的：找出一维搜索的最佳途径。

内容： 假设f (x )=0在[a ,b]区间内只有一个根（可以是重根），求解该方程等价于在 [a ,b ] 区间找|f (x )|的极小值点。

设计一种寻找极小值点的方法，使得计算f (x )的次数尽可能的少，并完成数值实验。

你能从理论上证明你的搜索方法最好吗？
思考题3、目的：了解牛顿收敛法的结构和“局部”收敛性。

内容：牛顿法可以直接用来求解复数方程z 3-1=0，在复平面上它的三个根是
z 1*=1 ，.2
1
23*3,2i z ±-
=。

选择中心位于坐标原点边长为2的正方形内的点为初始值，把收敛到三个不同根的初始值分别标上不同的颜色。

只要计算足够多的点，你将得到关于牛顿收敛域的彩色图片。

第三章 线性方程组的解法
一、主要要求
编写方程组求解的通用程序：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)}收敛的充
要条件是.1)(<B ρ。

定理2：若||B||<1,则迭代公式X (k)=BX (k-1)+f 收敛. 3、Jacobi 迭代公式：n i x a a a b x
n i
j j k ii
ij
ii i k i
,,2,11)
()1( =-=∑≠=+
J k J k f X B X Jacobi +=+)()1(:迭代的矩阵格式
其中：B J = -D -1(U+L)，f J =D -1b 。

4、Seidel 迭代公式：n i x a
a x a a a
b x
n
i j k j ii
ij
i j k j ii
ij ii i k i
,,2,11)(11)1()
1( =--=∑∑+=-=++
S k S k f X B X Seidel +=+)()1(:迭代的矩阵格式
其中：B s = -(D+L )-1U ，f s =(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-x 0||<=1.0e -6，则转到5），否则，转到3） 5）输出y 和n 。

Seidel 迭代算法： 1)给出矩阵A ，b ，x0
2)计算B = (D-L )-1U ，f =(D-L )-1
b ，则有：
3) y=BX 0+f
4)若 ||y-x 0||<=1.0e -6，则转到5），否则，转到3） 5）输出y 和n 。

6、程序框图：（见图4）
图4
仿Jacobi 迭代法算法和Seidel 迭代算法可给出超松驰迭代的算法（略）
三、数值计算实验（以下实验都需用Matlab 软件来完成） 实验3.1（求解线性方程组实验）
实验目的：掌握各种迭代法，比较各种迭代法的渐进收敛速度. 实验内容：令 1、计算A 的条件数cond(A)；（可选用任何一种范数）
2、上述方程组是否为病态方程组？若是，如何改变方程组的病态性？
3、分别用Jacobi 迭代、Seidel 迭代与超松驰迭代求方程组AX=b 的近似解；
⎪⎪⎪⎪
⎪⎭⎫
⎝⎛--------=231063
51710101768416104125A ⎪⎪⎪⎪⎪⎭
⎫
⎝⎛=31332332b
4、比较Jacobi 迭代、Seidel 迭代与超松驰迭代的渐进收敛速度。

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

四、具体操作见下例（以下实验都需用Matlab 软件来完成） 例：用Jacobi 方法求解下列方程组，设X(0)=0,精度为10-6
⎪⎩
⎪⎨⎧=+-=-+-=-6
10272109103
13
2
121x x
x x x x x 解：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为迭代次数。

五、思考题
思考题：目的：以Hilbert 矩阵为例，研究处理病态问题可能遇到的困难。

内容： 设Hilbert 矩阵为
⎪⎪
⎪⎪
⎪⎪
⎭
⎫
⎝
⎛-++++==)12/(1)2/(1)1/(1/1)2/(15
/14/13/1)1/(14/13/12
/1/13/12/11
)(n n n n n n n
h H ij n
它是一个对称正定矩阵，而且cond(H n )随着n 的增加迅速增加。

其逆矩阵H n -1=(a ij )，这
里
)!()!(])!1()!1([(1()!
1()!1()1(2
j n i n j i j i j n i n a j i ij -----+-+-+-=+
1） 画出ln(cond(H n )) ~ n 之间的曲线（可以用任何一种范数）。

你能猜出
ln(cond(H n )) ~ n 之间有何种关系吗？提出你的猜测并想法验证。

2） 设D 是H n 的对角元素开方构成的矩阵。

11ˆ--=D H D H n n ，不难看出它仍然
是对称矩阵，而且对角线元素都是1。

把H n 变成n H ˆ的技术称为预处理(preconditioning)。

画出ln(cond(
n H ˆ)/ cond(H n
)) ~ n 之间的曲线（可以用任何一种范数）。

对于预处理你能得
出什么印象？
3） 对于4≤n ≤12，给出不同的右端项b 。

分别用
X 1=H n -1b ，
b D H x n 11ˆˆ--=，
X 2=D -1x
ˆ以及X 3=H n \b(这是Matlab 的一条命令)求解H n X=b ，比较计算结果。

4） 取不同的n 并以H n 的第一列为右端向量b ，同高斯—塞德尔迭代求解H n X=b ，
观察其收敛性。

最后你能对于有关Hilbert 矩阵的计算得出哪些结论。

第四章 插值与拟合
一、主要要求
编写拉格朗日插值法、分段线性插值法、*三次样条插值通用程序(Matlab 自带)。

拉格朗日插值多项式： 二、主要结果回顾
插值法是函数逼近的重要方法之一，有着广泛的应用 。

在生产和实验中，函数f(x)或者其表达式复杂不便于计算或者无表达式而只有函数在给定点的函数值(或其导数值) ，此时我们希望建立一个简单的而便于计算的函数ϕ(x )，使其近似的代替f (x )，这就是所谓的插值法.有很多种插值法,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit 插值,分段插值和样条插值.
1、插值：求近似函数的方法:由实验或测量的方法得到所求函数 y=f (x ) 在互异点x 0 , x 1, ... , x n 处的值 y 0 , y 1 , … , y n ,
构造一个简单函数 ϕ(x) 作为函数 y=f (x ) 的近似表达式 ： y= f (x ) ≈ ϕ(x ) 使 ϕ(x 0)=y 0 , ϕ(x 1)=y 1 , ⋯, ϕ(x n)=y n , (*) 这类问题称为插值问题。

f(x) 称为被插值函数，ϕ(x ) 称为插值函数， x 0 , x 1, ... , x n 称为插值节点。

(*)式称为插值条件。

2、Lagrange 插值公式
∑∏
∑
=≠==+-+---=----------=n
j n j
i i j
i
j
i
j
n
j n j j j j j j j n j j n
y
x
x x y
x x x x x P
x x x x x x x 0
0011101110)()
)...()()...()(())...()()...()(()(
其中x 0 , x 1,... , x n 为插值节点，y j 为插值节点x j 处的函数值（j=1,2,…n ）。

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

4、拉格朗日插值法程序框图：（见图5）
（图5）
三、数值计算实验（以下实验都需利用Matlab 软件来完成） 实验4.1（观察龙格(Runge)现象实验）
实验目的：观察拉格朗日插值的龙格(Runge)现象. 实验内容：对于函数2
25
)(x a x f +=
进行拉格朗日插值，取不同的节点数，在区间[-5,5]
上取等距间隔的节点为插值点，把f (x )和插值多项式的曲线画在同一张图上进行比较。

（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、分析上述曲线图，你可以得出什么结论？
四、具体操作
例1、 给出f (x )=ln x 的数值表，用Lagrange 插值计算ln(0.54)的近似值。

解：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/2小时的温度值。

2）分别用分段线性插值，立方插值，三次样条插值和最邻近插值估计其值。

实验4.3 （插值与拟合实验）
实验目的：求下列数据的多项式拟合函数
x=[0 ,1, 2 ,3, 4 ,5, 6 ,7 ,8, 9 ,10];
y=[1 ,2.3, 2.1 ,2, 4.6 ,4.7, 4.3, 8.1 ,9.2, 9.8, 10.3];
实验内容：
1、做出原始数据的散点图；
2、做出原始数据的折线图（即分段线性插值函数图）；
3、做出原始数据的分段二次拟合曲线图（见图6）；
4、做出原始数据的分段三次拟合曲线图。

（图6）
五、思考题
思考题1、目的：观察最小二乘多项式的数值不稳定性现象。

内容：1、在[-1,1]区间上取n=20个等距节点，计算出以相应节点上e x 的值
做为数据样本，
2、以1,x ,x 2,...,x t 为基函数做l=3,5,7,9次的最小二乘多项式。

3、画出ln(cond(A)) ~n 之间的曲线，其中A 是确定最小二乘多项式系
数的矩阵。

4、计算出不同阶最小二乘多项式给出的最小偏差
.
))(()(21
i n
i i y x y l -=∑=σ
思考题2、目的：观察对于非光滑函数进行多项式插值的可能性。

内容：在[0，1]上取f(x)=|sink πx |，选择不同的k 和n ，用等距节点做拉格朗
日多插值项式，观察误差大小和收敛情况。

第五章 数值积分
一、主要要求
编写定步长复合梯形公式、定步长复合Simpson 公式、变步长梯形法及龙贝格算法*通用子程序；
二、主要结果回顾
1、梯形公式：对于连续函数f (x )，有积分中值定理
)),,()(()()(b a f a b dx x f b
a
∈-=⎰
ξξ
其中f (ξ)是被积函数f (x )在积分区间上的平均值。

因此，如果我们能给出求平均值f (ξ)的一种近似方法，相应地就可以得到计算定积分的一种数值方法。

如果取2
)
()()(b f a f f +≈
ξ，即得下列梯形公式：
)],()([2
)
()(b f a f a b dx x f b a +-≈⎰
二、辛普生(Simpson)公式：先用某个简单函数 近似逼近f (x )，然后用)(x ϕ在[a,b]区间的积分值近似表示f (x )在[a,b]区间上的定积分，即取
⎰⎰
≈b
a
b
a
dx x dx x f )()(ϕ
我们可以取)(x ϕ为前面介绍的f (x )的插值多项式Ln(x )或拟合多项式P(x )进行近似计算。

如取)(x ϕ为插值多项式Ln(x)，则相应得到的数值积分公式称为插值型求积公式。

如：
⎰⎰
≈b
a
n b
a
dx x L dx x f )()(
若考虑过点A,B,C 三点的抛物线段L 代替曲线段y=f (x )（见图7）
（图7）
L 就是以a,b,c 为节点构造二次插值多项式：
)()
)(()
)(()())(())(()())(())(()(b f c b a b c x a x c f b c a c b x a x a f b a c a b x c x x ----+----+----=
ϕ
取2b
a c +=
即得下列抛物线（辛普生）公式： )],()2
(4)([6)()(b f b a f a f a b dx x dx x f b a
b a +++-=≈⎰⎰ϕ
3、误差分析：
梯形公式的截断误差(余项)为:)),(()(''12
)(3
1b a f a b E ∈--
=ξξ 辛普生(Simpson)公式的截断误差(余项)：)),(()('''2280
)(5
2b a f a b E ∈--
=ξξ 注： 1)抛物线求积公式的计算精度高于梯形求积公式的计算精度。

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

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

• 从余项的讨论看到，积分区间越小，也可使求积公式的截断误差变小。

因此，我
们经常把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式，这就是复合求积公式的基本思想。

• 常用的复合求积公式是复合梯形公式和复合辛普生(Simpson)公式
4、复合梯形公式：
把区间[a,b] n 等分，取节点x i =a+ih ，i=0,1,...n, h=(b-a)/n ，对每个小区间[x i ,x i+1]用梯形求积公式，再累加起来得：
∑-==++≈1
1
)](2)()([2n i n i T x f b f a f h
I
5、复合辛普生(Simpson)公式：
把区间[a,b] n 等分，取节点x i =a+ih ，，i=0,1,...n, h=(b-a)/n ，对每个小区间[x i ,x i+1]用抛物线求积公式，再累加起来得：
)]()2
(4)([6110+-=+++=∑k n k k k n x f h
x f x f h I
（其中：2
h
x i +
为节点x i 与x i+1的中点，i=0,1,...n ） 6、变步长复合积分法
基本思想：是在求积过程中，通过对计算结果精度的不断估计，逐步改变步长（逐次分半），直到满足精度要求为止。

即按照给定的精度实现步长的自动选取。

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

利用梯形公式可得积分近似值：)],()([2
)
()(b f a f a b dx x f b
a
+-≈
⎰
记为T 0，即：)],()([2
0b f a f a
b T +-≈
此时步长h 0=b-a ，再将[a,b]二等分，即取步长2
)
(1a b h -=
使用n=2时的梯形公式可得积分近似值：
)]()2
(2)([4)(1b f b a f a f a b T +++-≈
),2
(2210a b a f a b T -+-+=
再将[a,b]四等分，即取步长2
22
)
(a b h -= 使用n=4时的梯形公式可得积分近似值：
)}()2)((2)({2)(2
3
`
132b f a b i a f a f a b T i +-++-≈∑= ,)](212[2)(2121
2
21∑=--+-+=i a b i a f a b T
再将[a,b]八等分，即取步长3
32)
(a b h -=
使用n=8时的梯形公式可得积分近似值：
)}()2)((2)({2)(3
7
`
143b f a b i a f a f a b T i +-++-≈∑= ,)](212[2)(212
213
32∑=--+-+=i a b i a f a b T 依次类推可得变步长梯形公式：
,...2,1,)](212[2)(211
21
1=--+-+≈∑-=-k a b i a f a b T T k i k
k k k 即：
)],()([2)
(0b f a f a b T +-≈
),2(2)(2101a b a f a b T T -+-+≈
,)](212[2)(21212212∑=--+-+≈i a b i a f a b T T
,)](212[2)(212
21
3
323∑=--+-+≈i a b i a f a b T T (*),...2,1,)](212[2)(211
21
1=--+-+≈∑-=-k a b i a f a b T T k i k
k k k 若预定精度为ε，以公式(*)计算积分的近似值，直到|T k+1-T k |<ε时终止计算,并以当前值
T k+1为所求近似值.
对Simpson 公式也可类似构造出相应的变步长积分公式。

（略）
7、算法：
a) 复合梯形公式：∑-==++≈1
1
)](2)()([2n i n i T x f b f a f h
I
b) 复合辛普生(Simpson)公式：)]()2
(4)([6110+-=+++=∑k n k k k n x f h
x f x f h I
算法：
1）令0),2
(,21=+=-=
I h
a f I n a
b h 2）对k=1,2,...,n-1，计算)(),2
(2211kh a f I I h
kh a f I I ++=+++= 3）))(24)((6
21b f I I a f h
I +++=。

c ）变步长梯形公式公式：)(2
21
2
12∑+
+=k n n n x
f h T T ，这里n
a
b h n -=。

算法：
1）输入精度ε，n=1，)]()([2
1b f a f a
b T +-=； 2）s=0；
3）对k=0,1,2,...,n-1计算：2
*),(h h k a x x f s s ++=+=； 4）)*(2
1
12s h T T +=
，n=2n ； 5）如果|T2-T1|<ε，则结束；否则执行6）； 6）2
h
h =
，T1=T2，转2）。

三、数值计算实验（以下实验都需利用Matlab 软件来完成） 实验5.1 数值积分实验
实验内容：用数值积分法求解下面的积分方程：
t
t d y t y 220
2
sin 24252
cos 2.38.6)
()(-=+-+⎰π
ττ
τ
1、 用梯形积分公式求解上述积分方程；
2、 用抛物线积分公式求解上述积分方程；
3、比较不同数值积分公式对方程解的影响。

实验5.2 数值微分实验
实验目的：熟悉MATLAB 中实现两点公式的函数diff ，称为微分或差分函数。

其主要调用格式有如下三种：
i ）diff （S ） 对由findsym 返回的自变量，求符号表达式S 的微分。

ii ）diff （S,‘v ’）或ddif （S,sym （‘v ’）） 对自变量v ，求表达式S 的微分。

iii ）diff （S,n ） 对整数n ，对表达式S 微分n 次。

diff （S,‘v ’,n ）或diff （S,‘v ’,n ） 这两种格式都可以被识别。

如：求sinx 2的一阶导数和t 6
的6阶导数。

输入
syms x,diff(sin(x^2))
syms t,diff(t^6,6)
ans =
2*cos(x^2)*x
ans =
720
实验内容：已知函数表Array试求f (x)在各节点的一、二阶导数的近似值。

四、思考题
思考题、目的：比较不同数值积分公式对方程解的影响。

内容：用n=4，利用复合梯形公式和复合Simpson公式，重复解实验中
的积分方程，比较其结果。

第六章 微分方程数值解
一、主要要求
编写Euler 法、改进的Euler 法、四阶Runge-Kutta 法(Matlab 自带)通用子程序；
二、主要结果回顾
解一阶常微分方程初值问题
⎩
⎨⎧=≤≤=y x y b
x a y x f y )(),('0
将区间[a,b]作n 等份，取步长n
a
b h -=
显式Euler 公式：),(1k k k k y x hf y y +=+ 隐式Euler 公式：)),(),((2
111+++++
=k k k k k k y x f y x f h
y y
改进的Euler 公式：⎪⎩⎪
⎨⎧++=+=++++))~,(),((2),(~1111k k k k k k k k k k y x f y x f h
y y y x hf y y
或表示为： ⎪⎪⎪
⎩
⎪
⎪⎪⎨⎧
+=+=+=++)(21),(),(11c p i p i i c i i i p y y y y x hf y y y x hf y y
三、数值计算实验（以下实验都需利用Matlab 软件来完成） 实验目的：比较几种数值解法的精度
实验内容：用Euler 法和Runge-Kutta 法求解初值问题：
)10(1
)0(2'≤≤⎪⎩⎪⎨⎧
=-
=x y y
x y y
四、具体操作
例：设有微分方程：
⎩⎨
⎧=≤≤+=0
00,)( y x y T x x y x y
1） 求出上述微分方程的解析表达式（精确解）； 2） 用Euler 法解微分方程（步长取h ）。

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

给出数据：x 0 , T, y 0 , h 输出(显示器)：数值解。

如 给出：0 0.4 1 0.1 输出
x y 0.1 1.10000 0.2 1.22000 0.3 1.36200 0.4 1.52820 解：1）由微分方程可得精确解为：1*2--=x e y x 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]； 1*2--=x e y x 得 =y [ 1.1103 1.2428 1.3997 1.5836]
而用Euler 法计算得 1y =[1.1050 1.2210 1.3492 1.4909] 作图：在Matlab 命令窗口输入： plot(x,y) hold on
plot(x,y1,'-r')
从图中可看出，用Euler 法计算出的结果与精确值有较大误差，还需改进。

五、思考题
思考题1、目的：用数值微分的数值解解决实际问题：
内容：1. 小型火箭初始质量为900千克，其中包括600千克燃料。

火箭竖直
向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力。

当燃料用尽时引擎关闭。

设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）。

重力加速度取9.8米/秒2.
A. 建立火箭升空过程的数学模型（微分方程）；
B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间
和高度。

思考题2、目的：用数值微分的数值解解决实际问题：
内容： 种群的数量（为方便起见以下指雌性）因繁殖而增加，因自然死亡
和人工捕获而减少。

记x k (t )为第t 年初k 岁（指满k-1岁，未满k 岁，下同）的种群数量，b k 为k 岁种群的繁殖率（1年内每个个体繁殖的数量），d k 为k 岁种群的死亡率（1年内死亡数量占总量的比例），h k 为k 岁种群的捕获量（1年内的捕获量）。

今设某种群最高年龄为5岁（不妨认为在年初将5岁个体全部捕获），b 1=b 2=b 5=0，b 3=2，b 4=4，d 1=d 2=0.3，d 3=d 4=0.2，h 1=400，h 2=200，h 3=150，h 4=100。

A. 建立x k (t+1)与x k (t )的关系（k=1,2,⋯5, t=0,1,⋯），如
11112)()()1(h t x d t x t x --=+。

为简单起见，繁殖量都按年初的种群数量
x k (t )计
算，不考虑死亡率。

B. 用向量T
t x t x t x ))(),(()(51 =表示t 年初的种群数量，用b k 和d k 定义
适当的矩阵L ，用h k 定义适当的向量h ，将上述关系表成h t Lx t x -=+)()1(的形
式。

C. 设t=0种群各年龄的数量均为1000，求t=1种群各年龄的数量。

又问设
定的捕获量能持续几年。

D. 种群各年龄的数量等于多少，种群数量x (t )才能不随时间t 改变。

E. 记D 的结果为向量x *, 给x * 以小的扰动作为x (0)，观察随着t 的增加x (t )
是否趋于x *, 分析这个现象的原因。

案例分析 ——估计水塔的水流量
该题是1991年美国大学生数学建模竞赛A 题（数据略有改动）
1、问题的提出
某地区用水管理机构需要对居民的用水速度（单位时间的用水量）和日总用水量进行估计。

现有一居民区，某自来水是由一个圆柱行水塔提供，水塔高12.2m ，塔的 直径为17.4m 。

水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。

按照设计，当水塔中的水位降至最低水位，约8.2m ，水泵自动启动加水；当水位升高到最高水位，约10.8m ，水泵停止工作。

表3.4给出的是某一天的测量记录数据，测量了28个时刻的数据，当由于水泵正向水塔供水，有3个时刻无法侧到水位（表中为—）。

表3.4 水塔中水位原始数据
时刻(t)/h 0 0.921 1.843 2.949 3.871 4.978 5.900 水位/m
9.677
9.479 9.308 9.125 8.982 8.814 8.686 时刻(t)/h 7.006 7.928 8.967 9.981 10.925 10.954 12.032 水位/m
8.525
8.388 8.220 — — 10.820 10.500 时刻(t)/h 12.954 13.875 14.982 15.903 16.926 17.931 19.037 水位/m
10.210
9.936 9.653 9.409 9.180 8.921 8.662 时刻(t)/h 19.959 20.839 22.015 22.958 23.880 24.986 25.908 水位/m
8.433
8.220
—
10.820
10.591
10.354
10.180
试建立数学模型，来估计居民的用水速度和日总用水量。

2、 问题的求解
1)水塔中水的体积的计算
由问题的要求，解决问题的关键是求出居民的用水速度，即单位时间内所用水的体积。

假设水塔是标准圆柱体，因此 h D V 24
π
=
（*）
式中，D 为水塔的直径（17.4m ）。

h 为水塔的水位高度。

计算水的流量，首先需要计算出水塔中水的体积，然后再由式（*）可计算出不同时刻水塔中水的体积V 。

输入
t=[ 0 0.921 1.843 2.949 3.871 4.978 5.900….。