数值分析课程报告
插值法和多项式拟合的研究
摘要
在科研和生产实践中，常常需要通过一组测量数据来寻找变量x与y的函数关系近似表达式。解决这类问题的方法有两种：一种是插值法,另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。拟合法能够是从给定的一组实验数据出发，寻找函数的一个近似表达式,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合。本文主要介绍拉格朗日插值法、埃尔米特插值法、三次样条插值法以及基于最小二乘法的多项式拟合。
关键词：拉格朗日插值，埃尔米特插值，样条插值，多项式拟合
1方法的意义
在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达，通常只是由观察与测试得到一些离散数值。有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析。解决这类问题的方法有两种：一种是插值法,另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。它要求给出函数的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定一个简单的函数()x ?作为()f x 的近似，概括地说，就是用简单函数为离散数组建立连续模型。插值法在实际应用中非常广泛，但是它也有明显的缺陷，一是测量数据常常带有测试误差，而插值多项式又通过所有给出的点，这样就是插值多项式保留了这些误差；二是如果实际得到的数据过多，则必然得到次数较高的插值多项式，这样近似的效果并不理想。拟合法能够很好的解决这些问题，它从给定的一组实验数据出发，寻找函数的一个近似表达式y=()x ?,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合的问题，函数的近似表达式y=()x ?称为拟合曲线。常用最小而二乘法来确定拟合曲线。
2插值法的介绍
2.1 插值法定义
设 f (x )为[a ，b ]上的函数，在互异点n x x x ,...,,10处的函数值分别为
)(),...,(),(10n x f x f x f ,构造一个简单函数 ?(x ) 作为函数 f (x ) 的近似表达式y = f (x ) ≈ ?(x ),使
)()(i i x f x =? , i ＝0, 1, 2, …,n (1.0) 则称?(x ) 为关于节点n x x x ,...,,10的插值函数；称n x x x ,...,,10 为插值节点；称
))((i i
x f x , i =1,2,… , n 为插值点；f (x ) 称为被插值函数。式（1.0）称为插值条
件。这类问题称为插值问题。插值的任务就是由已知的观测点,为物理量(未知量)建立一个简单的、连续的解析模型，以便能根据该模型推测该物理量在非观测点
处的特性。常用的插值函数类{()x ?}是代数多项式，相应插值问题是代数插值，本文主要介绍三种代数差值：拉格朗日插值，埃尔米特插值和样条插值。 2.2拉格朗日插值 2.2.1两点插值问题
已知)(i i x f y = )1,0(=i ，求满足插值条件i i y x P =)(1 )1,0(=i 的插值多项式
)(1x P 。
根据解析几何知识可知，所求的)(1x P 为过点),(),,(1100y x y x 的直线，即：
)()(00
10
101x x x x y y y x P ---+
= 上式经整理可改写为：
式中101
0)(x x x x x l --=
，0
101)(x x x x x l --=。显然{)1(0)0(1)(0===i i x l i ,{)1(1)0(0)(1===i i x l i ,
且)(0x l ，)(1x l 为由插值节点唯一确定的线性函数。
)(0x l ，)(1x l 为节点10,x x ，上的一次插值基函数。可以看出，节点10,x x 上的插值基函数的次数为插值节点个数减一，基函数组中所含的函数个数与插值节点数相同。而满足i i y x P =)(1 )1,0(=i 的插值多项式)(1x P 就是节点10,x x 上插值基函数的线性组合，其组合系数分别为10,y y 。这种表示为插值基函数线性组合的一次插值多项式也就是一次拉格朗日插值多项式。
当给定n+1个插值节点后，可类似定义n 次插值基函数，并以此构造n 次拉格朗日多项式。
2.2.2 n 次拉格朗日插值多项式
n 次插值基函数：
)
)...()()...(()
)...()()...(()(110110n k k k k k k n k k k x x x x x x x x x x x x x x x x x l --------=
+-+- )...,1,0(n k =
显然)(x l k 具有以下性质： 性质1，{
)
(0)
(1)(k i k i x l i k === )...,1,0(n k =
性质2，)(x l k )...,1,0(n k =为由插值节点n x x x ,...,,10唯一确定的n 次函数。
性质3， 基函数组所含的基函数个数与插值节点个数相同。
以n 次插值基函数为基础，可得拉格朗日插值多项式为：
k n
k n
k k k k k k n k k n
k k k n y x x x x x x x x x x x x x x x x y x l x L ∑
∑=+-+-=--------==01101100
))...()()...(()
)...()()...(()()(
()n
n i
k k i k i
i k
x x y x x ==≠-=-∑∏
若记∏=+-=n
i i n x x x 0
1)()(ω，则有∏≠=+-='
n
k
i i i k k n x x x 0
1)()(ω。于是)(x L n 也可写成：
k n
k k n k n n y x x x x x L ∑
=++'
-=0
11)
()()
()(ωω
2.2.3拉格朗日插值多项式的余项
)()()(x L x f x R n n -=称为拉格朗日插值多项式的余项，也叫截断误差。用简单的插值多项式)(x L n 代替复杂的函数)(x f ，这种做法是否有效，取决于截断误差是否满足精度要求。拉格朗日型余项：
)()!
1()
()()()(11x n f x L x f x R n n n n +++=-=ωξ
其中∏=+-=n
i i n x x x 01)()(ω，),(b a ∈ξ且依赖于0x 。一般情况下，),(b a ∈ξ的具体
数值无法知道，但是若能够求出1)1()(max ++≤≤=n n b
x a M x f ，则可以得出插值多项式的
截断误差限为：
)(max )!1()(11
x n M x R n b
x a n n +≤≤++≤
ω 由此看出，)(x R n 的大小除了与1+n M 有关外，还与插值节点有密切关系。当给定M 个点出的函数值，但仅选用其中)1(1m n n <++个作为插值条件而求某点
x 处函数值时，n+1个节点n x x x ,...,,10的选取应该尽可能的接近x
，使得计算的函数值的误差限尽可能的小。
2.3埃尔米特插值
许多实际问题不但要求插值多项式与被插值函数在节点处的函数值相当，而且还需要求其导数值相等。满足这种要求的插值多项式就是埃尔米特插值多项式。
一般情形的埃尔米特插值问题
一般情形的埃尔米特插值问题是指所满足的插值条件中函数值的个数与导数值的个数相等。即当函数)(x f 在区间],[b a 上n+1个节点i x ),...,1,0(n i =处的函数值)(i i x f y =及导数值)(i i x f m '=给定时，要求一个次数不超过2n+1的多项式
)(12x H n +，使之满足：
),...,1,0()()(1212n i m x H y x H i i n i i n =?
??
='
=++
这里给出个2n+2个插值条件，可唯一确定一个形式为
12121012...)(++++++=n n n x a x a a x H 的多项式，但是要确定2n+2个系数非常复杂，因此此时可以借用构造拉格朗日插值多项式的基函数。
设)(),(x x j j βα ),...,1,0(n j =为次数不超过2n+1的多项式，且满足：
),...,1,0,()(,0)(0)(,)(n j i x x x x ij i j i j i j ij i j =??
?
??='
=='=δββαδα
则满足上述条件的埃尔米特插值多项式可以写成用插值基函数表示的形式：
∑=++=n
j j j i j n m x y x H 012])()([βα
上式满足插值条件。
可以得到基函数)(),(x x j j βα的解析式为：
)(1)(21)(20x l x x x x x j n
j k k k
j j j ????
??????---=∑≠=α ),...,1,0(n j = )()()(2x l x x x j j j -=β ),...,1,0(n
j = 式中)(x l j 为拉格朗日插值基函数。 因此埃尔米特插值多项式即为：
j j n
j j i n
j j n j k k k
j j n m x l x x y x l x x x x x H )()()(1)(21)(2002012∑∑∑==≠=+-+????
??????---= 一般情形下的埃尔米特插值多项式的余项为：
)()!
22()()()()(12
)22(1212x n f x H x f x R n n n n +++++=-=ωξ
式中),(b a ∈ξ，且与x 有关。
埃尔米特插值的几何意义：曲线)(12x H y n +=与曲线)(x f y =在插值节点处有相同的公共切线。
在带导数的插值问题中，有时插值条件中的函数值个数与导数值个数不等。这时可以以一般情况的埃尔米特插值多项式为基础，运用待定系数法求出满足插值条件的多项式。
2.4样条插值
在上述方法中，我们根据区间],[b a 上给出的节点得以得到函数)(x f 的插值多项式，但是并非插值多项式的次数越高，逼近函数)(x f 的精度越好，主要原因是因为对于任意的插值节点，当∞→n 时，插值多项式)(x P n 不一定收敛到
)(x f 。这种高次插值不准确的现象称为龙格现象。为了避免高次插值的缺点，
人们常常采用分段插值的方法，即将插值区间分为若干个小区间，在每个小区间上运用前面介绍的插值方法构造低次插值多项式。采用分段现象插值与分段二次插值，可以构造一个整个连续的函数，而采用分段三次埃尔米特插值则可以构造一个整体上具有一节连续导数的插值函数。实际问题中，很少给出插值点上的导数值。三次样条插值就是在只给出插值点上函数值的情况下，构造一个整体上具有二阶连续导数的插值函数。
2.4.1三次样条插值函数的定义
设区间],[b a 上取n+1个节点b x x x x x a n n =<<<<<=+-11210...，若函数)(x S 满足条件：
（1）在整个区间],[b a 上具有二阶连续导数
（2）在每个小区间],[1i i x x - ),...,1,0(n i =上是x 的三次多项式
（3）i i y x S =)( ),...,1,0(n
i = 则称)(x S 为)(x f 的三次样条插值。)(x S 的边界条件为：
（1）给定两端点处的一阶导数值，记为：
00)(m x S ='，n n m x S =')(
（2）给定两端点处的二阶导数值，记为：
00)(M x S =''，n n M x S ='')(
此外，对0)()(0=''=''n x S x S 的边界条件称为自然边界条件。
2.4.2三次样条插值函数的构造
构造三次样条插值函数，就是要写出它在子区间],[1i i x x - ),...,1,0(n i =上的表达式，记为)(x S i ),...,1,0(n i =。
一、用节点处一阶导数表示的三次样条插值函数
记节点处的一阶导数值为i i m x S =')( ),...,1,0(n i =，若已知i m 后，则)(x S 在],[1i i x x - ),...,1,0(n i =上就是满足条件
11)(--=i i y x S ，i i y x S =)(，),...,1,0(n i = 11)(--='i i m x S ，i i m x S =')( ),...,1,0(n i =
的三次埃尔米特插值多项式。步骤如下：
（1）根据)(x S ，)(x S '在内阶段的连续性及插值条件，运用],[1i i x x -上的二点三次埃尔米特插值多项式，写出)(x S 用i m ),...,1,0(n i =表示的形式。
（2）利用)(x S ''在内节点i x )1,...,1,0(-=n i 的连续性及边界条件，导出含i m ),...,1,0(n i =的n+1阶线性方程组。
（3）求解含i m ),...,1,0(n i =的线性方程组，将得到的i m 代入],[1i i x x -上的二点三次埃尔米特插值多项式。 具体公式为：
i i i i i i i i i i i i i i i y x x x x x x x x y x x x x x x x x x S 2
11112
1112121)(????
??--???? ?
?--++???? ??--???? ??--+=------- +i i i i i i i i i i m x x x x x x m x x x x x x 2
1112
1)()(?
??
?
??---+???? ??------- i m 的确定由)(x S 得二阶导数在内节点连续的条件得到，即：
i i i i i i f m M m m =+++-112λ )1,...,1,0(-=n
i 该式即为n+1阶线性方程组。其中：
???
?
?
????
-+-=-=+=-=+-++++-)(31111111
1i i i i i i i i i i
i i i i i
i i i h y y h y y M f M h h h x x h λλλ )1,...,1,0(-=n i 求解该线性方程组还需两个边界条件。
对于第一种边界条件，将其代入，得到只含n-1个未知数的线性方程组：
????????
??--=???????? ?????
??
?????
??????--------n n n n n n n n n m M f f f m f m m m m M M M 112201112211222212222 λλλλ 对于第二种边界条件，将其代入，得到n+1阶线性方程组：
?????????
?
??=?????????? ???????
????
?
?
?----n n n n n n f f f f f m m m m m M M M 1210121011
2211212221
2
λλλ 其中：
???
?
???+-=--=-n n
n n n n M h h y y f M h h y y f 23231011010
二、用节点处二阶导数表示的三次样条插值函数
记节点处的二阶导数值为i i m x S ='')( ),...,1,0(n i =，由于)(x S ''在],[1i i x x -
),...,1,0(n i =上是x 的线性函数，因此构造以节点处二阶导数表示的三次埃尔
米特插值多项式步骤如下：
（1）根据)(x S ''在内节点的连续性及为线性函数的特点，将)(x S ''表示为线性函数。在根据)(x S 在内节点的连续性及插值条件，写出)(x S 用i M
),...,1,0(n i =表示的形式。
（2）利用)(x S '在内节点i x )1,...,1,0(-=n i 的连续性及边界条件，导出
含 i M ),...,1,0(n i =的n+1阶线性方程组。
（3）求解含i M ),...,1,0(n i =的线性方程组，将得到的i M 代入],[1i i x x -上的 )(x S 的表达式，即得到二点三次埃尔米特插值多项式。
具体公式为：
i
i i i i i i i i i i i i i i i i h x x h M y h x
x h M y M h x x M h x x x S 122113113)
6()6(6)(6)()(-------+--+-+-= i M 的确定由)(x S 得二阶导数在内节点连续的条件得到，即：
i i i i i i f M M M M =+++-112λ )1,...,1,0(-=n
i 该式即为n+1阶线性方程组。其中：
???
?
?
????---+=-=+=-=+-+++++-)(611111111
1i i i i i i i i i
i
i i i i i
i i i h y y h y y h h f M h h h M x x h λ )1,...,1,0(-=n i 解该方程还需要两个边界条件。
对于第一种边界条件，将其代入，得到只含n-1个未知数的线性方程组：
???????? ??=???????? ?????
??
?????
??????-----n n n n n n f f f f M M M M M M 11012211111212212
λλ 式中：
???
?
???
--=--=-)(6
)(61010110n n n n n n h y y m h f m h y y h f 对于第二种边界条件，将其代入，得到n+1阶线性方程组：
???????
?
??--=???????? ?????
??
?????
??????--------n n n n n n n n n M f f f M M f M M M M M M M 112201112211222212222λλλλ 对于三次样条插值函数来说，当插值节点逐渐加密时，可以证明：不但样条
插值函数收敛于函数本身，而且其导数也收敛于函数的导数。
3.曲线拟合
3.1 用最小二乘法求解矛盾方程组
工程实际中许多问题都可归结为矛盾方程组，实际中需要寻求矛盾方程组的未知数n x x x ,...,,21的一组取值，它使得偏差的绝对值之和∑=N
i i 1δ尽可能的小，为了
便于分析计算和应用，常采用使偏差的平方和2
1
1
1
2
)(∑∑∑===-==N i n
j i j ij N i i b x a Q δ达到
最小值，这一条件称为最小二乘原则。按最小二乘原则来选择n x x x ,...,,21的一组取值的方法称为求解矛盾方程组的最小二乘法，符合条件的n x x x ,...,,21的一组取值称为矛盾方程组的最小二乘解。
把Q 看出n 个自变量n x x x ,...,,21的二次函数，记为),...,,(21n x x x f Q =，因此，求矛盾方程组的最小二乘解就是二乘函数的最小值点。若矛盾方程组i
j n
j ij b x a =∑=1的系数矩阵A 的秩为n ，则二次函数),...,,(21n x x x f Q =一定存在最小值。即矛盾方程组的最小二乘解存在，且正则方程组b A Ax A T T =有唯一解，此解就是矛盾方程组的最小二乘解。
3.2多项式拟合
设通过测量得到函数)(x f y =的一组数据为：)(i i x f y = ),...,1,0(n i =，求一个次数低于N-1的多项式m m x a x a x a a x y ++++==...)(2210? )1(- ，使其最好的拟合这组数据，最好的标准时：使得)(x ?在i x 的偏差i i i y x -=)(?δ的平方和2 1 1 2 )]([∑∑==-==N i i i N i i y x Q ?δ达到最小。m a a a ,...,,10的解 即正则方程组b A Ax A T T =的解。 ????????? ??=++++=++++=++++∑∑∑∑∑∑∑∑∑∑∑∑∑∑===+=+===+=======N i i m i N i m i m N i m i N i m i N i m i N i i m i N i m i m N i i N i i N i i N i i N i m i m N i i N i i y x x a x a x a x a y x x a x a x a x a y x a x a x a N a 1121221111 01 11132121101 1122110......... 由该正则方程组求得的唯一解代入拟合多项式)(x y ?=，即为所求。 5.数值试验 表1是1971年到1990年我国总人口的统计数字，试根据1971年到1985年这15年人口统计数字用下面几种方法预测未来20年的人口数字，并用图示的方法比较1986年到1990年间预测人口数字与实际统计数字的差异，在你所使用的几种预测方法中找出一组较为合理的预测方法。 （1）指数形式bx ae y =； （2）拉格朗日插值、埃尔米特插值、样条插值； （3）三次多项式拟合； （4）四次多项式拟合 表1 人口统计数字 年份 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 人口 8.5229 8.7177 8.9211 9.0859 9.2420 9.3717 9.4974 9.6259 9.7542 9.8705 年份 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 人口 10.0072 10.1654 10.3008 10.4357 10.5851 10.7507 10.9300 11.1026 11.2704 11.4333 数值试验结果：1.指数形式拟合图形为： 2.三次多项式拟合的效果： 3.四次多项式拟合的效果 结论：基于最小二乘法的多项式拟合能取得良好的函数逼近效果。阶数越高，拟合效果越好。 参考文献 [1] 王正林，龚纯，等．精通MATLAB科学计算．北京：电子工业出版社，2009:135 —157 [2]任玉杰．数值分析及其MATLAB实现［M］．北京：高等教育出版社，2007： 584—642 [3]李信真，车刚明，等．计算方法[M]．西安：西北工业大学出版社，2010： 150—167 [4] 金一庆，陈越．数值方法[M].北京：机械工业出版社，2000：165—176 附录 1.按指数形式拟合的程序： x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542 ,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1 1.1026,11.2704,11.4333]; y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19 83,1984,1985,1986,1987,1988,1989,1990]; p=polyfit(x,log(y),1); a=exp(p(2)); b=p(1); y_n=a*exp(b*x); plot(x,y,x,y,'*',x,y_n),grid,axis([8.5,11.8,1970,1991]),legend('实际统计曲线','实际统计值','指数曲线') 2.按三次多项式拟合的程序： x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542 ,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1 1.1026,11.2704,11.4333]; y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19 83,1984,1985,1986,1987,1988,1989,1990]; p=polyfit(x,y,3); y1=polyval(p,x); plot(x,y,x,y,'*',x,y1),grid,axis([8.5,11.8,1970,1991]),legend('实际统计曲线','实际统计值','三次多项式拟合') 3.按四次多项式拟合的程序： x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542 ,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1 1.1026,11.2704,11.4333]; y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19 83,1984,1985,1986,1987,1988,1989,1990]; p=polyfit(x,y,4); y1=polyval(p,x); plot(x,y,x,y,'*',x,y1),grid,axis([8.5,11.8,1970,1991]),legend('实际统计曲线','实际统计值','四次多项式拟合') 数值分析课程设计题目与要求 （10级应数及创新班） [设计题一] 编写顺序Gauss消去法和列主元Gauss消去法的函数，再分别调用这两个函数求解下面的84阶方程组： = ， 然后考虑将方程组的阶数取为10至100之间多个值进行求解。将你的计算结果与方程组的精确解进行比较。从“快”、“准”、“省”三个方面分析以上两个算法，试提出改进的算法并加以实现和验证。 [设计题二] 编写平方根法和改进的平方根法（参见教材《计算方法》P54的例题2.5）的函数，然后分别调用这两个函数求解对称正定方程组Ax=b，其中A和b分别为： （1）系数矩阵A为矩阵（阶数取为10至100之间多个值）: ， 向量b随机地选取； （2）系数矩阵A为Hilbert矩阵（阶数取为5至40之间多个值），即A的第i行第j列元素，向量b的第i个分量取为。将你的计算结果与方程组的精确解进 行比较。 若出现问题，分析其原因，提出改进的设想并尝试实现之。 对于迭代法 ,......)2,1,0(99.021=-=+k x x x k k k ， 它显然有不动点0* =x 。试设计2个数值实验 得到收敛阶数的大概数值（不利用判定收敛阶的判据定理）： （1） 直接用收敛阶的定义； （2） 用最小二乘拟合的方法。 [设计题四] 湖水在夏天会出现分层现象，接近湖面温度较高，越往下温度变低。这种上热下冷的现象影响了水的对流和混合过程，使得下层水域缺氧，导致水生鱼类的死亡。如果把水温T 看成深度x 的函数T(x),有某个湖的观测数据如下： 环境工程师希望： 1） 用三次样条插值求出T(x)。 2） 求在什么深度处dx dT 的绝对值达到最大（ 即02 2=dx T d ）。 [设计题五] 某飞机头部的光滑外形曲线的型值点坐标由下表给出: ．．．值y 及一阶、二阶导数值y ’，y ”。绘出模拟曲线的图形。 淮海工学院计算机工程学院课程设计报告书 课程名：《数值分析》 题目：数值分析课程设计 班级： 学号： 姓名： 数值分析课程设计 课程设计要求 1、研究第一导丝盘速度y与电流周波x的关系。 2、数据拟合问题运用样条差值方法求出温度变化的拟合曲线。 课程设计目的 1、通过编程加深对三次样条插值及曲线拟合的最小二乘法的理解； 2、学习用计算机解决工程问题，主要包括数据处理与分析。 课程设计环境 visual C++ 6.0 课程设计内容 课程设计题目1： 合成纤维抽丝工段中第一导丝盘的速度对丝的质量有很大的影响，第一丝盘的速度和电流周波有重要关系。下面是一组实例数据： 其中x代表电流周波，y代表第一导丝盘的速度 课程设计题目3： 在天气预报网站上获得你家乡所在城市当天24小时温度变化的数据，认真观察分析其变化趋势，在此基础上运用样条差值方法求出温度变化的拟合曲线。然后将该函数曲线打印出来并与原来的温度变化数据形成的曲线进行比较，给出结论。写出你研究的心得体会。 课程设计步骤 1、利用最小二乘法写出题1的公式和算法； 2、利用excel表格画出数据拟合后题1的图像； 3、在Visual C++ 6.0中编写出相应的代码； 4、搜索11月12日南通当地一天的温度变化数据； 5、在Visual C++ 6.0中编写出相应的代码； 6、利用excel表格画出数据拟合后题3的图像 课程设计结果 课程设计题目1 数值拟合 解：根据所给数据，在excel窗口运行： x=[49.2 50.0 49.3 49.0 49.0 49.5 49.8 49.9 50.2 50.2] y=[16.7 17.0 16.8 16.6 16.7 16.8 16.9 17.0 17.0 17.1] 课程设计题目3 数据为：X=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]; Y=[12,12,11,12,12,12,12,12,13,15,16,17,17,18,17,17,17,16,15,15,15,15,14,14]; 源代码为： 第一题： #include 数值计算课程设计任务书 学院信息与计算科学/应用数学专业班级学生： 题目：典型数值算法的C++语言程序设计 课程设计从2017 年 6 月12 日起到2017 年7月 1 日 1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）： 每人需作10个算法的程序、必做6题、自选4题。 对每个算法要求用C++语言进行编程。 必选题： 1、高斯列主元法解线性方程组 2、牛顿法解非线性方程组 3、经典四阶龙格库塔法解一阶微分方程组 4、三次样条插值算法（压紧样条）用C++语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。 5、龙贝格求积分算法 6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。 自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次. 2、对课程设计成果的要求〔包括图表、实物等硬件要求〕： 2.1 提交课程设计报告 按照算法要求，应用C++语言设计和开发算法程序，提交由： 1）每个算法的原理与公式说明； 2）每个算法相应的程序设计说明（程序中的主要变量语义说明，变量的数据类型说明，数据在内存中组织和存储结构说明，各函数的输入形参和输出形参说明，函数功能说明，函数中算法主要流程图，函数的调用方法说明）； 3）每个程序使用的实例(引用的实例可以自拟，也可以借用相关数值计算参考书中的例题作为作为验证程序是否正确的实例，无论是自拟实例还是引用实例，实例都应详细写入报告的正文中)； 4）每个算法的调试记录(包括程序调试（静态调试和动态调试）和程序修改记录、程序测试（可以手工计算进行测试、也可以利用Matlab的函数或 研究生《数值分析》教学大纲 课程名称：数值分析 课程编号：S061005 课程学时：64 学时 课程学分： 4 适用专业：工科硕士生 课程性质：学位课 先修课程：高等数学，线性代数，计算方法，Matlab语言及程序设计 一、课程目的与要求 “数值分析”课是理工科各专业硕士研究生的学位课程。主要介绍用计算机解决数学问题的数值计算方法及其理论。内容新颖，起点较高，并加强了数值试验和程序设计环节。通过本课程的学习，使学生熟练掌握各种常用的数值算法的构造原理和过程分析，提高算法设计和理论分析能力，并且能够根据数学模型，提出相应的数值计算方法编制程序在计算机上算出结果。力求使学生掌握应用数值计算方法解决实际问题的常用技巧。 二、教学内容、重点和难点及学时安排： 第一章? 数值计算与误差分析( 4学时) 介绍数值分析的研究对象与特点，算法分析与误差分析的主要内容。 第一节数值问题与数值方法 第二节数值计算的误差分析 第三节数学软件工具----MATLAB 语言简介 重点：误差分析 第二章? 矩阵分析基础( 10学时) 建立线性空间、赋范线性空间、内积空间的概念，为学习以后各章打好基础。矩阵分解是解决数值代数问题的常用方法，掌握矩阵的三角分解、正交分解、奇异值分解，并能够编写算法程序。 第一节? 矩阵代数基础 第二节? 线性空间 第三节? 赋范线性空间 第四节? 内积空间和内积空间中的正交系 第五节矩阵的三角分解 第六节矩阵的正交分解 第七节矩阵的奇异值分解 难点：内积空间中的正交系。矩阵的正交分解。 重点：范数，施密特(Schmidt) 正交化过程，正交多项式，矩阵的三角分解, 矩阵的正交分解。 第三章? 线性代数方程组的数值方法( 12学时) 了解研究求解线性代数方程组的数值方法分类及直接法的应用范围。高斯消元法是解线性代数方程组的最常用的直接法，也是其它类型直接法的基础。在此方法基础上加以改进，可得选主元的高斯消元法、按比例增减的高斯消元法，其数值稳定性更高。掌握用列主元高斯消元法解线性方程组及计算矩阵的行列式及逆，并且能编写算法程序。掌握矩阵的直接三角分解法：列主元LU 分解，Cholesky分解。了解三对角方程组的追赶法的分解形式及数值稳定性的充分条件。掌握矩阵条件数的定义，并能利用条件数判别方程组是否病态以及对方程组的直接方法的误差进行估计。 迭代解法是求解大型稀疏方程组的常用解法。熟练掌握雅可比迭代法、高斯- 塞德尔迭代法及SOR 方法的计算分量形式、矩阵形式，并能在计算机上编出三种方法的程序用于解决实际问题。了解极小化方法：最速下降法、共轭斜量法。迭代法的收敛性分析是研究解线性代数方程组的迭代法时必须考虑的问题。对于上述常用的迭代法，须掌握其收敛的条件。而对一般的迭代法，掌握其收敛性分析的基本方法和主要结果有助于进一步探究新的迭代法。 第一节求解线性代数方程组的基本定理 第二节高斯消元法及其计算机实现 第三节矩阵分解法求解线性代数方程组 第三节? 误差分析和解的精度改进 第四节? 大型稀疏方程组的迭代法 第五节? 极小化方法 难点：列主元高斯消元法，直接矩阵三角分解。迭代法的收敛性，雅可比迭代法，高斯-塞德尔迭代法，SOR 迭代法。 数学与计算科学学院实验报告 实验项目名称常微分方程数值解 所属课程名称数值方法B 实验类型验证 实验日期 2013.11.11 班级 学号 姓名 成绩 【实验过程】（实验步骤、记录、数据、分析） 注：以下图形是通过Excel 表格处理数据得出，并未通过MATLAB 编程序所得！ 1、1(0)1dy y x dx y ?=-++???=? 由题可知精确解为：x y x e -=+，当x=0时，y(x)=0。 h=0.1 表1 h=0.1时三个方法与精确值的真值表 图1 h=0.1时三个方法走势图 步长 Euler 法 预估校正法 经典四阶库 精确值 0.1 1.010000 1.005000 1.004838 1.249080 0.2 1.029000 1.019025 1.018731 1.055455 0.3 1.056100 1.041218 1.040818 1.091217 0.4 1.090490 1.070802 1.070320 1.131803 0.5 1.131441 1.107076 1.106531 1.176851 0.6 1.178297 1.149404 1.148812 1.226025 0.7 1.230467 1.197211 1.196586 1.279016 0.8 1.287421 1.249975 1.249329 1.335536 0.9 1.348678 1.307228 1.306570 1.395322 1.0 1.413811 1.368541 1.367880 1.458127 h=0.05（此时将源程序中i的围进行扩大，即for(i=0;i<20;i++)） 表2 h=0.05时三个方法与精确值的真值表步长Euler法预估校正法经典四阶库精确值 0.05 1.002500 1.001250 1.001229 1.011721 0.10 1.007375 1.004877 1.004837 1.024908 0.15 1.014506 1.010764 1.010708 1.039504 0.20 1.023781 1.018802 1.018731 1.055455 0.25 1.035092 1.028885 1.028801 1.072710 0.30 1.048337 1.040915 1.040818 1.091217 0.35 1.063421 1.054795 1.054688 1.110931 0.40 1.080250 1.070436 1.070320 1.131801 0.45 1.098737 1.087752 1.087628 1.153791 0.50 1.118800 1.106662 1.106531 1.176851 0.55 1.140360 1.127087 1.126950 1.200942 0.60 1.163342 1.148954 1.148812 1.226025 0.65 1.187675 1.172193 1.172046 1.252062 0.70 1.213291 1.196736 1.196585 1.279016 0.75 1.240127 1.222520 1.222367 1.306852 0.80 1.268121 1.249485 1.249329 1.335536 0.85 1.297215 1.277572 1.277415 1.365037 0.90 1.327354 1.306728 1.306570 1.395322 0.95 1.358486 1.336900 1.336741 1.426362 1.00 1.390562 1.368039 1.367880 1.458127 图2 h=0.05时三个方法走势图 数值分析学习感想 一个学期的数值分析，在老师的带领下，让我对这门课程有了深刻的理解和感悟。这门 课程是一个十分重视算法和原理的学科，同时它能够将人的思维引入数学思考的模式，在处 理问题的时候，可以合理适当的提出方案和假设。他的内容贴近实际，像数值分析，数值微 分，求解线性方程组的解等，使数学理论更加有实际意义。 数值分析在给我们的知识上，有很大一部分都对我有很大的帮助，让我的生活和学习有 了更加方便以及科学的方法。像第一章就讲的误差，在现实生活中，也许没有太过于注意误 差，所以对误差的看法有些轻视，但在学习了这一章之后，在老师的讲解下，了解到这些误 差看似小，实则影响很大，更如后面所讲的余项，那些差别总是让人很容易就出错，也许在 别的地方没有什么，但是在数学领域，一个小的误差，就很容易有不好的后果，而学习了数 值分析的内容，很容易就可以将误差锁定在一个很小的范围内，在这一范围内再逼近，得出 的近似值要准确的多，而在最开始的计算中，误差越小，对后面的影响越小，这无疑是好的。 数值分析不只在知识上传授了我很多，在思想上也对我有很大的影响，他给了我很多数 学思想，很多思考的角度，在看待问题的方面上，多方位的去思考，并从别的例子上举一反三。像其中所讲的插值法，在先学习了拉格朗日插值法后，对其理解透彻，了解了其中 的原理和思想，再学习之后的牛顿插值以及三次样条插值等等，都很容易的融会贯通，很容 易的就理解了其中所想，他们的中心思想并没有多大的变化，但是使用的方式却是不同的， 这不仅可以学习到其中心内容，还可以去学习他们的思考方式，每个不同的思考方式带来的 都是不同的算法。而在看待问题上，不同的思考方式总是可以快速的全方位的去看透彻问题， 从而知道如何去解决。 在不断的学习中，知识在不断的获取，能力在不断的提升，同时在老师的不懈讲解下， 我逐渐的发现数值分析所涵盖的知识面特别的广泛，而我所需要学习的地方也更加的多，自 己的不足也在不断的体现，我知道这只是我刚刚接触到了数学的那一角，在以后我还会接触 到更多，而这求知的欲望也在不停的驱赶我，学习的越多，对今后的生活才会有更大的帮助。 计算132 2013014923 张霖篇二：数值分析学习报告 数值分析学习心得报告 班级：11级软工一班 姓名： * * * 学号: 20117610*** 指导老师：* * * 学习数值分析的心得体会 无意中的一次选择，让我接触了数值分析。 作为这学期的选修课，我从内心深处来讲，数值分析真的有点难。感觉它是在高等数学 和线性代数的基础上，又加深了探讨。虽然这节课很难，我学的不是很好，但我依然对它比 较感兴趣。下面就具体说说我的学习体会，让那些感兴趣的同学有个参考。 学习数值分析，我们首先得知道一个软件——matlab。matrix laboratory，即矩阵实验 室，是math work公司推出的一套高效率的数值计算和可视化软件。它是当今科学界最具影 响力、也是最具活力的软件，它起源于矩阵运算，并高速发展成计算机语言。它的优点是强 大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面、便捷的与其他程序和语 言接口。 根据上网搜集到的资料，你就会发现matlab有许多优点： 首先，编程简单使用方便。到目前为止，我已经学过c语言，机器语言，java语言，这 《数值分析》课程设计实验报告 龙格—库塔法分析Lorenz 方程 200820302033 胡涛 一、问题叙述 考虑著名的Lorenz 方程 () dx s y x dt dy rx y xz dt dz xy bz dt ?=-???=--???=-?? 其中s ，r ，b 为变化区域内有一定限制的实参数，该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。 二、问题分析 Lorenz 方程实际上是一个四元一阶常微分方程，用解析法精确求解是不可能的，只能用数值计算，最主要的有欧拉法、亚当法和龙格- 库塔法等。为了得到较高精度的，我们采用经典四阶龙格—库塔方法求解该问题。 三、实验程序及注释 （1）算法程序 function [T]=Runge_Kutta(f,x0,y0,h,n) ％定义算法，其中f 为待解方程组， x0是初始自变量，y0是初始函数 值，h 是步长，n 为步数 if nargin<5 n=100; ％如果输入参数个数小于5，则步数 n=100 end r=size(y0);r=r(1); ％返回初始输出矩阵的行列数，并将 值赋给r(1) s=size(x0);s=s(1); ％返回初始输入矩阵的行列数，并 将值赋给s(1) r=r+s; T=zeros(r,n+1); T(:,1)=[y0;x0]; for t=2:n+1 ％以下是具体的求解过程 k1=feval(f,T(1:r-1,t-1)); k2=feval(f,[k1*(h/2)+T(1:r-1,t-1);x0+h/2]); k3=feval(f,[k2*(h/2)+T(1:r-1,t-1);x0+h/2]); k4=feval(f,[k3*h+T(1:r-1,t-1);x0+h]); x0=x0+h; T(:,t)=[T(1:r-1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0]; end 学生实验报告实验课程名称 开课实验室 学院年级专业班 学生姓名学号 开课时间至学年学期 if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;%设置为默认格式显示,显示5位 （2）建立MATLAB界面 利用MA TLAB的GUI建立如下界面求解线性方程组： 详见程序。 五、计算实例、数据、结果、分析 下面我们对以上的结果进行测试，求解： ? ? ? ? ? ? ? ? ? ? ? ? - = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - - - - 7 2 5 10 13 9 14 4 4 3 2 1 13 12 4 3 3 10 2 4 3 2 1 x x x x 输入数据后点击和，得到如下结果： 更改以上数据进行测试，求解如下方程组： 1 2 3 4 43211 34321 23431 12341 x x x x ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果： 课 程 设 计 我再也回不到大二了， 大学是那么短暂 设计题目 数值分析 学生姓名 李飞吾 学 号 x x x x x x x x 专业班级 信息计x x x x x 班 指导教师 设 计 题 目 共15题如下 成绩 数值分析课程设计 1．1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？（15621） 试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题 解：算法分析：解该问题主要使用递推算法，关于椰子数目的变化规律可以设起初的椰子数为0p ，第一至五次猴子在夜里藏椰子后，椰子的数目分别为01234,,,,p p p p p 再设最后每个人分得x 个椰子,由题: 14 (1)5 k k p p +=- （k=0,1,2,3,4）51(1)5 x p =- 所以551p x =+,11k k p p +=+利用逆向递推方法求解 15 1,4 k k p p +=+ （k=0,1,2,3,4） MATLAB 代码： n=input('n= '); n= 15621 for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end if p==fix(p), break end end disp([x,p]) 1．2 设，1 5n n x I dx x =+? （1）从0I 尽可能精确的近似值出发，利用递推公式： 11 5(1,2,20)n n I I n n -=-+= 计算机从1I 到20I 的近似值； （2）从30I 较粗糙的估计值出发，用递推公式： 数值分析课程设计题目与要求 （12级应数及创新班） [设计题一] 编写顺序Gauss消去法和列主元Gauss消去法的函数，再分别调用这两个函数求解下面的84阶方程组： = ， 然后考虑将方程组的阶数取为10至100之间多个值进行求解。将你的计算结果与方程组的精确解进行比较。从“快”、“准”、“省”三个方面分析以上两个算法，试提出改进的算法并加以实现和验证。 [设计题二] 编写平方根法和改进的平方根法（参见教材《计算方法》P54的例题2.5）的函数，然后分别调用这两个函数求解对称正定方程组Ax=b，其中A和b分别为： （1）系数矩阵A为矩阵（阶数取为10至100之间多个值）: ， 向量b随机地选取； （2）系数矩阵A为Hilbert矩阵（阶数取为5至40之间多个值），即A的第i行第j列元素，向量b的第i个分量取为。将你的计算结果与方程组的精确解进 行比较。 若出现问题，分析其原因，提出改进的设想并尝试实现之。 对于迭代法 ,......)2,1,0(99.02 1=-=+k x x x k k k ， 它显然有不动点0*=x 。试设计2个数值实验 得到收敛阶数的大概数值（不利用判定收敛阶的判据定理）： （1） 直接用收敛阶的定义； （2） 用最小二乘拟合的方法。 [设计题四] 湖水在夏天会出现分层现象，接近湖面温度较高，越往下温度变低。这种上热下冷的现象影响了水的对流和混合过程，使得下层水域缺氧，导致水生鱼类的死亡。如果把水温T 看成深度x 的函数T(x),有某个湖的观测数据如下： 环境工程师希望： 1） 用三次样条插值求出T(x)。 2） 求在什么深度处dx dT 的绝对值达到最大（ 即02 2=dx T d ）。 [设计题五] 某飞机头部的光滑外形曲线的型值点坐标由下表给出: ．．．值y 及一阶、二阶导数值y ’，y ”。绘出模拟曲线的图形。 机电工程学院 机械工程 陈星星 6720150109 《数值分析》课程设计实验报告 实验一 函数插值方法 一、问题提出 对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n ==。试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。 数据如下： （1 求五次Lagrange 多项式5L ()x ，计算(0.596)f ,(0.99)f 的值。（提示：结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈） 实验步骤： 第一步：先在matlab 中定义lagran 的M 文件为拉格朗日函数 代码为： function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end 第二步：然后在matlab命令窗口输入： >>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >>p = lagran(x,y) 回车得到： P = 121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845 由此得出所求拉格朗日多项式为 p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步：在编辑窗口输入如下命令： >> x=[0.4 0.55 0.65 0.80,0.95 1.05]; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.0845; >> plot(x,y) 命令执行后得到如下图所示图形，然后 >> x=0.596; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.084 y =0.6257 得到f（0.596）=0.6257 同理得到f（0.99）=1.0542 插值法和多项式拟合的研究 摘要 在科研和生产实践中，常常需要通过一组测量数据来寻找变量x与y的函数关系近似表达式。解决这类问题的方法有两种：一种是插值法,另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。拟合法能够是从给定的一组实验数据出发，寻找函数的一个近似表达式,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合。本文主要介绍拉格朗日插值法、埃尔米特插值法、三次样条插值法以及基于最小二乘法的多项式拟合。 关键词：拉格朗日插值，埃尔米特插值，样条插值，多项式拟合 1方法的意义 在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达，通常只是由观察与测试得到一些离散数值。有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析。解决这类问题的方法有两种：一种是插值法,另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。它要求给出函数的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定一个简单的函数()x ?作为()f x 的近似，概括地说，就是用简单函数为离散数组建立连续模型。插值法在实际应用中非常广泛，但是它也有明显的缺陷，一是测量数据常常带有测试误差，而插值多项式又通过所有给出的点，这样就是插值多项式保留了这些误差；二是如果实际得到的数据过多，则必然得到次数较高的插值多项式，这样近似的效果并不理想。拟合法能够很好的解决这些问题，它从给定的一组实验数据出发，寻找函数的一个近似表达式y=()x ?,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合的问题，函数的近似表达式y=()x ?称为拟合曲线。常用最小而二乘法来确定拟合曲线。 2插值法的介绍 2.1 插值法定义 设 f (x )为[a ，b ]上的函数，在互异点n x x x ,...,,10处的函数值分别为 )(),...,(),(10n x f x f x f ,构造一个简单函数 ?(x ) 作为函数 f (x ) 的近似表达式y = f (x ) ≈ ?(x ),使 )()(i i x f x =? , i ＝0, 1, 2, …,n (1.0) 则称?(x ) 为关于节点n x x x ,...,,10的插值函数；称n x x x ,...,,10 为插值节点；称 ))((i i x f x , i =1,2,… , n 为插值点；f (x ) 称为被插值函数。式（1.0）称为插值条 件。这类问题称为插值问题。插值的任务就是由已知的观测点,为物理量(未知量)建立一个简单的、连续的解析模型，以便能根据该模型推测该物理量在非观测点 第一章 1霍纳（horner）方法： 输入=c + bn*c bn?1*c b3*c b2*c b1*c an an?1 an?2 ……a2 a1 a0 bn bn?1 bn?2 b2 b1 b0 answer p（x）=b0 该方法用于解决多项式求值问题=anxn+an?1xn?1+an?2xn?2+……+a2x2+a1x+a0 ? 2 注：p为近似值 p（x） 绝对误差： ?|ep?|p?p ?||p?p rp? |p| 相对误差： ?|101?d|p?p rp?? |p|2 有效数字: （d为有效数字，为满足条件的最大整数） 3 big oh(精度的计算)： o(h?)+o(h?)=o(h?); o(hm)+o(hn)=o(hr) [r=min{p,q}]; o(hp)o(hq)=o(hs) [s=q+p]; 第二章 2.1 求解x=g(x)的迭代法用迭代规则 ，可得到序 列值｛｝。设函数g 满足 y 定义在得 。如果对于所有 x ，则函数g 在 ，映射y=g(x)的范围 内有一个不动点; 此外，设 ，存在正常数k<1，使 内，且对于所有x，则函数g 在 内有唯一的不动点p。 ，（ii）k是一个正常数， 。如果对于所有 定理2.3 设有（i）g，g ’（iii ） 如果对于所有x在 这种情况下，p成为排斥不动点，而且迭代显示出局部发散 性。波理 尔 查 . 诺 二 分 法 ( 二 分 法 定) <收敛速度较慢> 试值（位）法：<条件与二分法一样但改为寻求过点(a,f(a))和(b,f(b))的割线l与 x轴的交点(c,0)> 应注意 越来越 小，但可能不趋近于0，所以二分法的终止判别条件不适合于试值法 . f(pk?1) 其中k=1,2,……证明：用 f(pk?1) 牛顿—拉夫森迭代函数：pk?g(pk?1)?pk?1? 泰勒多项式证明 第三章线性方程组的解法对于给定的解线性方程组ax=b a11x1 ? a12x2 ? ? ? a1nxn ? b1 a21x1 ? a22x2 ? ? ? a2nxn ? b2 ? an1x1 ? an2x2 ? ? ? annxn ? bn 一gauss elimination （高斯消元法第一步forward elimination 第二步 substitution 二lu factorization 第一步 a = lu 原方程变为lux=y ； 第二步令ux=y,则ly = b由下三角解出y；第三步 ux=y,又上三角解出x ； 三iterative methods（迭代法） a11x1 ? a12x2 ? ? ? a1nxn ? b1 a21x1 ? a22x2 ? ? ? a2nxn ? b2? ） back 初始值 0,x0,?,x0x1n2 四 jacobi method 1.选择初始值 2.迭代方程为 0,x0,?,x0x1n2 k?1? x1k?1 ? x2 本文主要通过Matlab 软件，对数值分析中的LU 分解法、最小二乘法、复化Simpon 积分、Runge-Kutta 方法进行编程，并利用这些方法在MATLAB 中对一些问题进行求解，并得出结论。 实验一线性方程组数值解法中，本文选取LU 分解法，并选取数据于《数值分析》教材第5章第153页例5进行实验。所谓LU 分解法就是将高斯消去法改写为紧凑形式，可以直接从矩阵A 的元素得到计算L 、U 元素的递推公式，而不需要任何步骤。用此方法得到L 、U 矩阵，从而计算Y 、X 。 实验二插值法和数据拟合中，本文选取最小二乘拟合方法进行实验，数据来源于我们课堂学习该章节时的课件中的多项式拟合例子进行实验。最小二乘拟合是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。利用excel 的自带函数可以较为方便的拟合线性的数据分析。 实验三数值积分中，本文选取复化Simpon 积分方法进行实验，通过将复化Simpson 公式编译成MATLAB 语言求积分∫e ;x dx 1 0完成实验过程的同时，也对复化Simpon 积分章节的知识进行了巩固。 实验四常微分方程数值解，本文选取Runge-Kutta 方法进行实验，通过实验了解Runge-Kutta 法的收敛性与稳定性同时学会了学会用Matlab 编程实现Runge-Kutta 法解常微分方程，并在实验的过程中意识到尽管我们熟知的四种方法，事实上，在求解微分方程初值问题，四阶法是单步长中最优秀的方法，通常都是用该方法求解的实际问题，计算效果比较理想的。 实验五数值方法实际应用，本文采用最小二乘法拟合我国2001年到2015年的人口增长模型，并预测2020年我国人口数量。 关键词：Matlab ；LU 分解法；最小二乘法；复化Simpon 积分；Runge-Kutta 课程设计报告 课程名称数值分析 课题名称数值积分 专业信息与计算科学 班级 学号 姓名 指导教师 2015 年12 月20 日 湖南工程学院 课程设计任务书 课程名称数值分析 课题数值积分 专业班级信息与计算科学0901班 学生姓名 学号 指导老师辉 审批 任务书下达日期2015 年12 月7 日任务完成日期2015 年12 月20日 设计内容与设计要求 1. 设计内容： 非奇异矩阵矩阵A ∈R n*n ，已知A -1的一个近似矩阵D （0）∈R n*n ，则由矩阵公式： ?????+=-=--)()1()1(K K K K K F I D D AD I F , K=0，1，2，3........... （1）.已知矩阵A 及其逆矩阵的一个近似D (k)为： A=?? ??? ?? ?? ???--------7.49.43.49.19.47.11.88.78.26 .21.27.07.37.08.38.1 D= ???? ? ???? ???---------185.0061.0388.0293.0199.0009.0046.0230.0089.0016.0169.0035.0270.0163.0460.0211.0 用以上方法计算序列{D (k)}迭代次数超过100次时结束。 （2）分析最后得到的D (k)是否A 的一个较好的近似逆矩阵 2.设计要求： ● 课程设计报告正文内容 a. 问题的描述及算法设计； b. 算法的流程图（要求画出模块图）； c. 算法的理论依据及其推导； d. 相关的数值结果（通过程序调试），； e. 数值计算结果的分析； f. 附件（所有程序的原代码，要求对程序写出必要的注释）。 ● 书写格式 计算方法教案新疆医科大学 数学教研室 张利萍 一、课程基本信息 1、课程英文名称：Numerical Analysis 2、课程类别：专业基础课程 3、课程学时：总学时54 4、学分：4 5、先修课程：《高等数学》、《线性代数》、《Matlab 语言》 二、课程的目的与任务： 计算方法是信息管理与信息系统专业的重要理论基础课程，是现代数学的一个重要分支。其主要任务是介绍进行科学计算的理论方法，即在计算机上对来自科学研究和工程实际中的数学问题进行数值计算和分析的理论和方法。通过本课程的学习，不仅使学生初步掌握数值分析的基本理论知识，而且使学生具备一定的科学计算的能力、分析问题和解决问题的能力，为学习后继课程以及将来从事科学计算、计算机应用和科学研究等工作奠定必要的数学基础。 三、课程的基本要求： 1．掌握计算方法的常用的基本的数值计算方法 2．掌握计算方法的基本理论、分析方法和原理 3．能利用计算机解决科学和工程中的某些数值计算应用问题，增强学生综合运用知识的能力 4．了解科学计算的发展方向和应用前景 四、教学内容、要求及学时分配： (一) 理论教学： 引论（2学时） 第一讲（1-2节） 1．教学内容： 计算方法(数值分析)这门课程的形成背景及主要研究内容、研究方法、主要特点；算法的有关概念及要求；误差的来源、意义、及其有关概念。数值计算中应注意的一些问题。 2．重点难点： 算法设计及其表达法；误差的基本概念。数值计算中应注意的一些问题。3．教学目标： 了解数值分析的基本概念；掌握误差的基本概念：误差、相对误差、误差限、相对误差限、有效数字；理解有效数字与误差的关系。学会选用相对较好的数值计算方法。 数值分析实验指导书 实验一 1．1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？ 试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题（15621）。 1．2 设，1 05n n x I dx x =+? （1）从0I 尽可能精确的近似值出发，利用递推公式： 11 5(1,2,20)n n I I n n -=-+= 计算机从1I 到20I 的近似值； （2）从30I 较粗糙的估计值出发，用递推公式： 111 (30,29,,3,2)55n n I I n n -=-+= 计算从1I 到20I 的近似值； （3）分析所得结果的可靠性以及出现这种现象的原因。 1．3 绘制Koch 分形曲线 问题描述：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形（图1-4）；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形（图1-5），这时，图形中共有17个结点。这种迭代继续进行下去可以形成Koch 分形曲线。在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。Koch 分形曲线的绘制与算法设计和计算机实现相关。 数值分析实验 ——MATLAB实现 姓名sumnat 学号2013326600000 班级13级应用数学2班 指导老师 2016年1月 一、插值：拉格朗日插值 (1) 1、代码： (1) 2、示例： (1) 二、函数逼近：最佳平方逼近 (2) 1、代码： (2) 2、示例： (2) 三、数值积分：非反常积分的Romberg算法 (3) 1、代码： (3) 2、示例： (4) 四、数值微分：5点法 (5) 1、代码： (5) 2、示例： (6) 五、常微分方程：四阶龙格库塔及Adams加速法 (6) 1、代码：四阶龙格库塔 (6) 2、示例： (7) 3、代码：Adams加速法 (7) 4、示例： (8) 六、方程求根：Aitken 迭代 (8) 1、代码： (8) 2、示例： (9) 七、线性方程组直接法：三角分解 (9) 1、代码： (9) 2、示例： (10) 八、线性方程组迭代法：Jacobi法及G-S法 (11) 1、代码：Jacobi法 (11) 2、示例： (12) 3、代码：G-S法 (12) 4、示例： (12) 九、矩阵的特征值及特征向量：幂法 (13) 1、代码： (13) 2、示例： (13) 一、插值：拉格朗日插值 1、代码： function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例： >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5 实验报告 实验项目名称数值积分与数值微分实验室数学实验室 所属课程名称数值逼近 实验类型算法设计 实验日期 班级 学号 姓名 成绩 实验概述： 【实验目的及要求】 本次实验的目的是熟练《数值分析》第四章“数值积分与数值微分”的相关内容，掌握复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式。 本次试验要求编写复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的程序编码，并在MATLAB软件中去实现。 【实验原理】 《数值分析》第四章“数值积分与数值微分”的相关内容，包括：复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的相应算法和相关性质。 【实验环境】（使用的软硬件） 软件： MATLAB 2012a 硬件： 电脑型号：联想 Lenovo 昭阳E46A笔记本电脑 操作系统：Windows 8 专业版 处理器：Intel（R）Core（TM）i3 CPU M 350 @2.27GHz 2.27GHz 实验内容： 【实验方案设计】 第一步，将书上关于复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的内容转化成程序语言，用MATLAB实现；第二步，分别用以上求积公式的程序编码求解不同的问题。 【实验过程】（实验步骤、记录、数据、分析） 实验的主要步骤是：首先分析问题，根据分析设计MATLAB程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。 实验：用不同数值方法计算积分 (1) 取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？ (2) 用龙贝格求积计算完成问题（1）。 (3)用勒让德多项式确定零点，再代入计算高斯公式，使其精度达到10-4 （1）在MATLAB的Editor中建立一个M-文件，输入程序代码，实现复合梯形求积公式的程序代码如下：数值分析课程设计题目与要求
数值分析课程设计
数值计算课程设计任务书
研究生《数值分析》教学大纲
数值分析实验报告之常微分方程数值解
数值分析学习心得体会.doc
《数值分析》课程设计报告
数值分析实验报告
数值分析课程课程设计汇总
12级数值分析课程设计
数值分析实验报告资料
数值分析课程报告
数值分析学习方法
数值分析课程设计(最终版)
数值分析-课程设计doc
数值分析每节课的教学重点、难点
数值分析课程设计实验指导书
数值分析实验— MATLAB实现
数值分析实验报告3