1、1第一章 绪论一、主要要求通过实验,认真理解和体会数值计算的稳定性、精确性与步长的关系。 二、主要结果回顾:1、算法:电子计算机实质上只会做加、减、乘、除等算术运算和一些逻辑运算,由这些基本运算及运算顺序规定构成的解题步骤,称为算法它可以用框图、算法语言、数学语言或自然语言来描述。用计算机算法语言描述的算法称为计算机程序。 (如 c语言程序,c+语言程序,Matlab 语言程序等) 。2、最有效的算法:应该运算量少,应用范围广,需用存储单元少,逻辑结构简单,便于编写计算机程序,而且计算结果可靠。3、算法的稳定性:一个算法如果输入数据有误差,而在计算过程中舍入误差不增长,则称此算法是数值稳定的,
2、否则称此算法为不稳定的。换句话说:若误差传播是可控制的,则称此算法是数值稳定的,否则称此算法为不稳定的。4、控制误差传播的几个原则:1)防止相近的两数相减;2)防止大数吃小数;3)防止接近零的数做除数;4)要控制舍入误差的累积和传播;5)简化计算步骤,减小运算次数,避免误差积累。三、数值计算实验(以下实验都需利用 Matlab 软件来完成)实验 1.1(体会数值计算精度与步长关系的实验)实验目的:数值计算中误差是不可避免的,要求通过本实验初步认识数值分析中两个重要概念:截断误差和舍入误差,并认真体会误差对计算结果的影响。问题提出:设一元函数 f :RR,则 f 在 x0 的导数定义为: hxf
3、h)(lim)( 00实验内容:根据不同的步长可设计两种算法 ,计算 f 在 x0 处的导数。计算一阶导数的算法有两种:(1)hxff )()00(2)fff2( 000请给出几个计算高阶导数的近似算法,并完成如下工作:1、对同样的 h,比较(1)式和(2)式的计算结果;2、针对计算高阶导数的算法,比较 h 取不同值时(1)式和 (2)式的计算结果。2实验要求:选择有代表性的函数 f(x)(最好多选择几个) ,利用 Matlab 提供的绘图工具画出该函数在某区间的导数曲线 f(s)(x),再将数值计算的结果用 Matlab 画出来,认真思考实验的结果。实验分析:不论采用怎样的算法,计算结果通常
4、都会有误差。比如算法(1),由 Taylor 公式,知: )(!2)()( 30000 hOxfhxffhxf 所以有 200ff利用上式来计算 f (x0),其截断误差为: )(!221hxhT所以误差是存在的,并且当步长 h 越来越小时,(1)的近似程度也越来越好。类似地可以分析(2)的截断误差为: )(!330)(2Of上述截断误差的分析表明(2)是比(1)更好的算法,因为对步长 h(0,n=1,2,当 n=1 时,得: x1当 n2 时,由分部积分可得: n=2,3,10nnIeI另外,还有: 110dxexIn实验内容:由递推关系 In=1-nIn-1,可得计算积分序列I n的两种算
5、法:算法一、直接使用递推公式得: eI1In=1-nIn-1 n=2,3算法二、利用递推公式变形得: ,.3210nIInN实验要求:用上述两种算法分别在计算中采用 5 位、6 位和 7 位有效数字,请判断哪种4算法给出的结果更精确。实验分析:两种算法的优劣可能与你的第一感觉完全不同。设算法一中 I1 的误差为 e1,由 I1 递推计算 In 的误差为 en算法二中 IN 的误差为 N,由 IN 向前递推计算 In (n0输出失败信息n=1y1*y newton(x0) 即可求出 求方程 fc1=0 在 x0 附近的近似解 得出结果ans =0.5018 3.1407 6.2832 9.424
6、89图 3注:用 Newton 迭代法可求出多个根。具体做法:先用 fplot 命令作出函数的图形,再根据图形给出不同的初始值 x0,最后通过调用 Newton 迭代法程序求出非线性方程的所有根。五、思考题思考题 1、 在二分法中取区间中点,每次计算一次函数值。如果每次把区间三等份,计算两个函数值,是否可以找出方程的根?给出它的算法。和二分法比较它的效率如何?思考题 2、目的:找出一维搜索的最佳途径。内容: 假设 f(x)=0 在a,b区间内只有一个根(可以是重根) ,求解该方程等价于在 a,b 区间找|f(x)| 的极小值点。设计一种寻找极小值点的方法,使得计算f(x)的次数尽可能的少,并完
7、成数值实验。你能从理论上证明你的搜索方法最好吗?思考题 3、目的:了解牛顿收敛法的结构和“局部”收敛性。内容:牛顿法可以直接用来求解复数方程 z3-1=0,在复平面上它的三个根是z1*=1 , 。选择中心位于坐标原点边长为 2 的正方形内的点为.213*3,2iz初始值,把收敛到三个不同根的初始值分别标上不同的颜色。只要计算足够多的点,你将得到关于牛顿收敛域的彩色图片。10第三章 线性方程组的解法一、主要要求编写方程组求解的通用程序:Jacobi 迭代、Seidel 迭代以及 SOR 迭代程序 二、主要结果回顾1、迭代法: 它的基本思想是将线性方程组 AX=b 化为 : X=BX+f 再由此构
8、造向量序列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.0e-6x0=y y=B*x0+fn=n+1输出x结束图 4仿 Jacobi 迭代法算法和 Seidel 迭代算法可给出超松驰迭代的算法(略)三、数值计算实验(以下实验都需用 Matlab 软件来完成
9、)实验 3.1(求解线性方程组实验)实验目的:掌握各种迭代法,比较各种迭代法的渐进收敛速度.实验内容:令1、计算 A 的条件数 cond(A);(可选用任何一种范数)2、上述方程组是否为病态方程组?若是,如何改变方程组的病态性? 3、分别用 Jacobi 迭代、Seidel 迭代与超松驰迭代求方程组 AX=b 的近似解;23106578462 312b124、比较 Jacobi 迭代、Seidel 迭代与超松驰迭代的渐进收敛速度。注:上述实验分两次完成。四、具体操作见下例(以下实验都需用 Matlab 软件来完成)例:用 Jacobi 方法求解下列方程组,设 X(0)=0,精度为 10-661
10、027291031xx解: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.99380.93810.6938n = 9注: n=9 为迭代次数。五、思考题思考题:目的:以 Hilbert 矩阵为例,研究处理病态问题可能遇到的困难。内容: 设 Hilbert 矩阵为 )12/()/(1)/(1/5/4/3/ )(2/13/1/1)( nnnnhHijn 它是一个对称正定矩阵,而且 cond(Hn)随着 n 的
11、增加迅速增加。其逆矩阵 Hn-1=(aij),这里 )!(!1()1(2jijiji jajiij 1) 画出 ln(cond(Hn) n 之间的曲线(可以用任何一种范数) 。你能猜出ln(cond(Hn) n 之间有何种关系吗?提出你的猜测并想法验证。2) 设 D 是 Hn 的对角元素开方构成的矩阵。1DHnn,不难看出它仍13然是对称矩阵,而且对角线元素都是 1。把 Hn 变成 n的技术称为预处理(preconditioning)。画出 ln(cond( n)/ cond(Hn) n 之间的曲线(可以用任何一种范数) 。对于预处理你能得出什么印象?3) 对于 4n12,给出不同的右端项 b
12、。分别用 X1=Hn-1b, bDHx1,X 2=D-1 x以及 X3=Hnb(这是 Matlab 的一条命令)求解 HnX=b,比较计算结果。4) 取不同的 n 并以 Hn 的第一列为右端向量 b,同高斯塞德尔迭代求解HnX=b,观察其收敛性。最后你能对于有关 Hilbert 矩阵的计算得出哪些结论。14第四章 插值与拟合 一、主要要求编写拉格朗日插值法、分段线性插值法、*三次样条插值通用程序(Matlab 自带) 。拉格朗日插值多项式:二、主要结果回顾插值法是函数逼近的重要方法之一,有着广泛的应用 。在生产和实验中,函数f(x)或者其表达式复杂不便于计算或者无表达式而只有函数在给定点的函数
13、值(或其导数值) ,此时我们希望建立一个简单的而便于计算的函数 (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 , (*
14、)这类问题称为插值问题。 f(x) 称为被插值函数,(x) 称为插值函数, x0 , x1, . , xn 称为插值节点。(*)式称为插值条件。2、Lagrange 插值公式 njji ji jnj njjjjjjjyxyPxx00 1110)( ).()(.()其中 x0 , x1,. , xn 为插值节点, yj 为插值节点 xj 处的函数值(j=1,2,n ) 。3、Lagrange 插值的截断误差(插值余项)定理: 设 Ln(x)是过点 x0 ,x 1 ,x2 ,xn 的 n 次插值多项式, f (n)(x)在a ,b上连续,f (n+1)(x)在a ,b 上存在,其中a,b 是包含点
15、 x0 ,x1 ,x2 ,,xn 的任一区间,则对任意给定的 xa,b,总存在一点 (a,b) (依赖于 x)使)()!()() 1)(nxLfRnnnf其中: 。.()101nx4、拉格朗日插值法程序框图:(见图 5)15开始输入x i ,yii =0,1,2,.,nL=0,k=0T=1T=T(x-xi)/(xk-xi) i=0,1,2,.,k-1, k+1,.,nL=L+ykTk=n k=k+1输出L结束NY(图 5)三、数值计算实验(以下实验都需利用 Matlab 软件来完成)实验 4.1(观察龙格(Runge)现象实验)实验目的:观察拉格朗日插值的龙格(Runge)现象.实验内容:对于
16、函数 进行拉格朗日插值,取不同的节点数,在区间-25)(xaf5,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)=lnx 的数值表,用 Lagrange 插值计算 ln(0.54)的近似值。16x 0.4 0.5 0
17、.6 0.7 0.8Lnx -0.916291 -0.693147 -0.510826 -0.357765 -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 分段插值实验实验目的:分段插值计算,会使用一维插值函数 ,寻找最佳的插值方法。
18、实验内容:在 112 的 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、做出
19、原始数据的分段二次拟合曲线图(见图 6) ;4、做出原始数据的分段三次拟合曲线图。(图 6)17五、思考题思考题 1、目的:观察最小二乘多项式的数值不稳定性现象。内容:1、在-1,1区间上取 n=20 个等距节点,计算出以相应节点上 ex 的值做为数据样本,2、以 1,x,x2,.,xt 为基函数做 l=3,5,7,9 次的最小二乘多项式。3、画出 ln(cond(A) n 之间的曲线,其中 A 是确定最小二乘多项式系数的矩阵。4、计算出不同阶最小二乘多项式给出的最小偏差.)()(21iniiyxl思考题 2、目的:观察对于非光滑函数进行多项式插值的可能性。内容:在0,1上取 f(x)=|si
20、nkx| ,选择不同的 k 和 n,用等距节点做拉格朗日多插值项式,观察误差大小和收敛情况。18第五章 数值积分一、主要要求编写定步长复合梯形公式、定步长复合 Simpson 公式、变步长梯形法及龙贝格算法*通用子程序;二、主要结果回顾1、梯形公式:对于连续函数 f (x),有积分中值定理 ),()bafbdba 其中 f ( )是被积函数 f (x)在积分区间上的平均值。因此,如果我们能给出求平均值 f ()的一种近似方法,相应地就可以得到计算定积分的一种数值方法。如果取 ,即得下列梯形公式:2f),(2)()bfabdxfba 二、辛普生(Simpson)公式:先用某个简单函数 近似逼近
21、f (x),然后用 在a,b区)(x间的积分值近似表示 f (x)在a,b区间上的定积分,即取 babad我们可以取 为前面介绍的 f (x)的插值多项式 Ln(x)或拟合多项式 P(x)进行近似计算。如取 为插值多项式 Ln(x),则相应得到的数值积分公式称为插值型求积公式。如:)(xbanbadxLxf)()(若考虑过点 A,B,C 三点的抛物线段 L 代替曲线段 y=f (x)(见图 7)19ABC a b2o xy(图 7)L 就是以 a,b,c 为节点构造二次插值多项式: )()()()()() bfcabxcfacxfbacxx 取 即得下列抛物线(辛普生)公式:2bac ),(2
22、(4)6)()( fffdxxfbaa3、误差分析:梯形公式的截断误差(余项)为: ),()(12(31 bafabE辛普生(Simpson)公式的截断误差 (余项): ),()(8052f注: 1)抛物线求积公式的计算精度高于梯形求积公式的计算精度。2)当积分区间较大时,利用上述公式计算积分产生的误差也较大。所以在建立积分公式时,应选择小区间和抛物线求积公式。 从余项的讨论看到,积分区间越小,也可使求积公式的截断误差变小。因此,我们经常把积分区间分成若干小区间,在每个小区间上采用次数不高的插值公式,如梯形公式或抛物线公式,构造出相应的求积公式,然后再把它们加起来得到整个区间上的求积公式,这就
23、是复合求积公式的基本思想。20 常用的复合求积公式是复合梯形公式和复合辛普生(Simpson)公式4、复合梯形公式:把区间a,b n 等分,取节点 xi=a+ih,i=0,1,.n, h=(b-a)/n,对每个小区间x i,xi+1用梯形求积公式,再累加起来得:1)(2)(2ni niTfbfahI5、复合辛普生(Simpson)公式:把区间a,b n 等分,取节点 xi=a+ih, ,i=0,1,.n, h=(b-a)/n,对每个小区间x i,xi+1用抛物线求积公式,再累加起来得: )()2(4)(6110 knkkkfhffhI(其中: 为节点 xi 与 xi+1 的中点,i=0,1,.
24、n )2hxi6、变步长复合积分法基本思想:是在求积过程中,通过对计算结果精度的不断估计,逐步改变步长(逐次分半) ,直到满足精度要求为止。即按照给定的精度实现步长的自动选取。下面以梯形公式为例来说明这种方法。利用梯形公式可得积分近似值: ),(2)()bfabdxfba 记为 T0,即: ,(20f此时步长 h0=b-a,再将a,b二等分,即取步长 )(1h使用 n=2 时的梯形公式可得积分近似值: )()2()(4(1 bfaffabT,20b再将a,b四等分,即取步长 )(ah使用 n=4 时的梯形公式可得积分近似值: )(2)(2)(2(313 bfabiffbTi ,)(1122ia
25、f21再将a,b八等分,即取步长 32)(abh使用 n=8 时的梯形公式可得积分近似值: )(2)()(2( 37143 bfabiffbTi ,)(12133iaf依次类推可得变步长梯形公式: ,.21,)(212)(11 kabifbTkikkk即: ),(0fa2(21b,)(1)121i abiafT ,)(2)(213323iifb (*),.21,)(2)(11 kabiafTkikk若预定精度为 ,以公式(*)计算积分的近似值,直到|T k+1-Tk|espif feval(fun,a)*feval(fun,c)0 b=c;c=(a+b)/2; elseiffeval(fun,c)*feval(fun,b)0 a=c;c=(a+b)/2;