数值分析2016上机实验报告

序言数值分析是计算数学的范畴，有时也称它为计算数学、计算方法、数值方法等，其研究对象是各种数学问题的数值方法的设计、分析及其有关的数学理论和具体实现的一门学科，它是一个数学分支。

是科学与工程计算（科学计算）的理论支持。

许多科学与工程实际问题（核武器的研制、导弹的发射、气象预报）的解决都离不开科学计算。

目前，试验、理论、计算已成为人类进行科学活动的三大方法。

数值分析是计算数学的一个主要部分,计算数学是数学科学的一个分支,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。

现在面向数值分析问题的计算机软件有：C,C++,MATLAB,Python,Fortran等。

MATLAB是matrix laboratory的英文缩写，它是由美国Mathwork公司于1967年推出的适合用于不同规格计算机和各种操纵系统的数学软件包，现已发展成为一种功能强大的计算机语言，特别适合用于科学和工程计算。

目前，MATLAB应用非常广泛，主要用于算法开发、数据可视化、数值计算和数据分析等，除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。

本实验报告使用了MATLAB软件。

对不动点迭代，函数逼近（lagrange插值，三次样条插值，最小二乘拟合），追赶法求解矩阵的解，4RungeKutta方法求解，欧拉法及改进欧拉法等算法做了简单的计算模拟实践。

并比较了各种算法的优劣性，得到了对数值分析这们学科良好的理解，对以后的科研数值分析能力有了极大的提高。

目录序言 (1)问题一非线性方程数值解法 (3)1.1 计算题目 (3)1.2 迭代法分析 (3)1.3计算结果分析及结论 (4)问题二追赶法解三对角矩阵 (5)2.1 问题 (5)2.2 问题分析（追赶法） (6)2.3 计算结果 (7)问题三函数拟合 (7)3.1 计算题目 (7)3.2 题目分析 (7)3.3 结果比较 (12)问题四欧拉法解微分方程 (14)4.1 计算题目 (14)4.2.1 方程的准确解 (14)4.2.2 Euler方法求解 (14)4.2.3改进欧拉方法 (16)问题五四阶龙格-库塔计算常微分方程初值问题 (17)5.1 计算题目 (17)5.2 四阶龙格-库塔方法分析 (18)5.3 程序流程图 (18)5.4 标准四阶Runge-Kutta法Matlab实现 (19)5.5 计算结果及比较 (20)问题六舍入误差观察 (22)6.1 计算题目 (22)6.2 计算结果 (22)6.3 结论 (23)7 总结 (24)附录问题一 非线性方程数值解法1.1 计算题目编写不动点迭代法求根程序：把方程010423=-+x x 写成至少四种x=g （x ）的形式，取初值5.1x 0=，进行不动点迭代求根，并比较收敛性及收敛速度。

1.2 迭代法分析将方程f （x ）=0改写成其等价形式)(x x φ=取方程根的某一近似值x0作为初始点，由函数)(01x x φ=可计算出x 1，如此下去，设当前点为x k ，有）（x φ计算出x k+1，即)(1k k x x φ=+ k=0,1，.....称为迭代公式。

（收敛条件）设)(x φ在[a ，b]上有连续的一阶导数，并满足：].,[,1)()2(;)(],,[)1('b a x L x b x a b a x ∈∀≤≤≤≤∈∀φφ有则有函数φ(x)在区间[a,b]上存在唯一的不动点（方程的根）x*；对任何x0属于[a,b]，可由迭代公式得到序列{x k }均收敛到方程的根x*。

设上述条件成立时，算法的中止条件为：k k k x k kx k x x L L x x x x L L x x --≤---≤-+--1*01*11 现将方程010423=-+x x 改写成如下四种x=φ（x ）形式并计算收敛性。

(1)21031k k x x -=+ 计算得φ’(x0) 不收敛 (2)321410x x k -=+计算得φ’(x0) 不收敛 (3)k k k x x x 41021+=+ 计算得φ’(x0) 不收敛 (4)4101+=+k k x x 计算得φ’(x0) 收敛 1.3计算结果分析及结论在command 窗口输入func=inline(‘φ(x)’)；[y k ]=StablePoint(1.5,func) 函数相对误差为3-10*1。

计算结果如下表1-1：从计算结果看到，迭代法（2）（3）均不收敛，因为他们不满足局部收敛条件，而迭代法（4）比迭代法（1）收敛快，只需要四步就可以计算得到近似值。

在做不动点迭代时，为使误差尽可能小且数据稳定。

由局部收敛性定理，在将函数f(x)化作x=φ（x ）时，应尽可能构造函数使φ（x ）收敛。

表1-1 迭代法计算结果问题二追赶法解三对角矩阵2.1 问题编写有效程序解线性方程组Ax=b，其中2.2 问题分析（追赶法）我们利用矩阵的直接三角分解法来推导三对角的计算公式，由系数矩阵A的特点，可以将方程分解为两个三角矩阵的乘积，即A=LU。

其中L 为下三角矩阵，U为单位上三角矩阵。

用追赶法求解严格占优矩阵Ax=f 等价于解两个三角方程组：①Ly=f，求y。

②Ux=y，求x。

从而得到解三对角方程组的追赶法公式：1)计算{βi}的递推公式β1=c1/b1，Βi=ci/（bi-aiβi-1），i=2,3,4，n-1；2)解Ly=fy1=f1/b1,Yi=(fi-aiyi-1)/(bi-aiβi-1)，i=2,3，...，n；3)解Ux=yXn=yn ，Xi=yi-βixi+1，i=n-1，n-2，...，2,12.3 计算结果当b=（1,2,3,4,5，...，50）T 时，即初始值为matrix=[1:1:50]，计算得x=（0.3333,0.6667，....，13.1186）。

只需要不断更换matrix 控制输入b 向量，就可以得到解x 。

MA TLAB 程序为附录程序2。

问题三 函数拟合3.1 计算题目 对函数22511)(xx f +=在区间[-1,1]上取xi=-1+0.2i （i=0.1....，10） （a ）对函数做多项式插值和三次样条插值，并画出插值函数及f （x ）的函数；（b ）对函数求其三次拟合曲线并画出拟合曲线的图像，与(a)中结果进行比较。

3.2 题目分析3.2.1 Lagrange 插值对于插值函数)(x ϕ，我们通常可以选择多种不同的函数类型，但由于代数多项式具有简单和一些良好的特性，我们常选用代数多项式作为插值函数.首先我们来看这样一个问题：给定两个插值点),(),,(1,100y x y x 其中,10x x ≠怎样做通过这两点的一次插值函数？过两点作一条直线，这条直线就是通过这两点的一次多项式插值函数，简称线性插值.下面先用待定系数法构造插值直线.设直线方程为,)(101x a a x L +=将)(),,(1,100y x y x 分别代入直线方程)(1x L ，得⎩⎨⎧=+=+11100010y x a a y x a a , 当10x x ≠时,因01110≠x x 所以方程组有解，且解唯一.这也表明，平面上两个点有且仅有一条直线通过，用待定系数法构造插值多项式的方法简单直观，容易看到解的存在性和唯一性，但要解一个方程组才能得到插值函数的系数，因工作量大且不便向高阶推广，故这种构造方法不宜采用.当10x x ≠时，若用两点式表示这条直线，则有： 101001011)(y x x x x y x x x x x L --+--= 这种形式称为Lagrange 插值多项式. 记)(),(,)(,)(1001011010x l x l x x x x x l x x x x x l --=--=称为插值基函数，计算),(),(10x l x l 的值,可知 .,0,1)(⎩⎨⎧≠===ji j i x l ij j i δ 在Lagrange 插值多项式中，可将)(1x L 看作两条直线0101y x x x x --与1010y x x x x --的叠加，并可看到两个插值点的作用和地位是平等的。

如果我们给定三个插值点2,1,0)),(,(=i x f x i i ，其中i x 互不相等，那么该怎样构造函数)(x f 的二次（抛物线）插值多项式呢？仿照线性插值的Lagrange 插值，我们可设)(),()()()()()()(2211002x l x f x l x f x l x f x l x L i ++=为二次函数。

对)(0x l 来说，要求21,x x 是它的零点，因此可设),)(()(210x x x x A x l --=同理)(),(21x l x l 也有相应形式。

),())(()())(()())(()(2101200212x f x x x x C x f x x x x B x f x x x x A x L --+--+--=∴ 将210,,x x x x =分别代入，可得 ))((1,))((1,))((1120221012010x x x x C x x x x B x x x x A --=--=--= 有)())(())(()())(())(()())(())(()(2120210121012002010212x f x x x x x x x x x f x x x x x x x x x f x x x x x x x x x L ----+----+----=一般地，当给定n+1个互不相同的插值节点时，就可得出函数的n 次插值多项式: )())(()()())(()()()()(11011000n i i i i i i n i i n i i n i i i n x x x x x x x x x x x x x x x x x f x f l x L --------==+-+-==∑∑ 下面我们以定理的形式来给出Lagrange 插值多项式的误差估计。

设)(x f 在区间[]b a ,上有直到n+1阶导数，n x x x ,,,10 是[]b a ,上n+1个互异节点，)(x P n 满足)()(i i n x f x p =的n 次插值多项式，则对[]b a x ,∈∀，有)()!1()()(1)1(x n f x R n n n +++=ωξ，其中()b a x x x n i i n ,,)()(01∈-=∏=+ξω，且依赖于.x 。

3.2.2 三次样条插值所谓三次样条插值多项式)(x S n 是一种分段函数，它在节点i x 011()n n a x x x x b -=<<⋅⋅⋅<<=分成的每个小区间1[,]i i x x -上是3次多项式，其在此区间上的表达式如下：22331111111()[()()]()()666[,]1,2,,.i i i i i i i i i i i i i i ii i h x x h x x S x x x M x x M y M y M h h h x x x i n --------=-+-+-+-∈=⋅⋅⋅，， 因此，只要确定了i M 的值，就确定了整个表达式，i M 的计算方法如下： 令：11111111116()6(,,)i i i i i i i i i i i i i i i i i i i i i h h h h h h y y y y d f x x x h h h h μλμ++++--+++⎧===-⎪++⎪⎨--⎪=-=⎪+⎩， 则i M 满足如下n-1个方程：1121,2,,1i i i i i i M M M d i n μλ-+++==⋅⋅⋅-，对于第一种边界条件下有⎪⎪⎩⎪⎪⎨⎧-=+-=+---000110111)'],([62]),['(62h f x x M M h x x f f M M n n n n n n 如果令,]),['(6,1,)'],[(6,111000100---==-==n n n n n n h x x f f d h f x x f d μλ那么解就可以为 ⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d M M M M 110110111102222 μλμλμλ3.2.3 曲线拟合由已知的离散数据点选择与实验点误差最小的曲线)(...)()()(1100x a x a x a x S n n ϕϕϕ+++=称为曲线拟合的最小二乘法。