数值分析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

x k k

x k x x L

x x x x L

L x x --≤---≤-+--1*0

11 现将方程010423=-+x x 改写成如下四种x=φ（x ）形式并计算收敛性。 (1)21031k k x x -=+ 计算得φ’(x0) 不收敛 (2)321410x x k -=+ 计算得φ’(x0) 不收敛

(3)k

k k x x x 410

1+=+ 计算得φ’(x0) 不收敛 (4)4

1+=

+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 迭代法计算结果

k Φ1(x) φ2(x) φ3(x) φ4(x)

0 1 2 3 4 5 6 7 8 9 10 ... 1.500000 1.500000 1.500000 1.500000 1.286953 0.333333 1.212121 1.348399 1.402540 3.185185 1.582848 1.367376 1.345458 -10.1938 1.131630 1.364957 1.375170 -135.220 1.722026 1.365264 1.360094 -24375.9 1.014869

1.367846 ... 1.964853

1.363887 ... 0.853237

1.365916 ...

2.414895

1.364878 ... 0.645523

1.365410 ... 3.334673

问题二追赶法解三对角矩阵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=f

y1=f1/b1,

Yi=(fi-aiyi-1)/(bi-aiβi-1)，i=2,3，...，n；

3)解Ux=y

Xn=yn ，

Xi=yi-βixi+1，i=n-1，n-2，...，2,1

2.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 计算题目

对函数2

2511

)(x

x 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 ，得

???=+=+1

1100

010y x a a y x a a , 当10x x ≠时,因01110≠x x

所以方程组有解，且解唯一.这也表明，平面上两个点有且仅有一条直线通过，用待定系数法构造插值多项式的方法简单直观，容易看到解的存在性和唯一性，但要解一个方程组才能得到插值函数的系数，因工作量大且不便向高阶推广，故这种构造方法不宜采用.

当10x x ≠时，若用两点式表示这条直线，则有： 10

01011)(y x x x x y x x x x x L --+--=

这种形式称为Lagrange 插值多项式. 记)(),(,)(,)(100

101101

0x l x l x x x x x l x x x x x l --=--=称为插值基函数，计算),(),(10x l x l 的值,可知

.,0,1)(???≠===j

i j i x l ij j i δ

在Lagrange 插值多项式中，可将)(1x L 看作两条直线0101

y x x x x --与10

10y 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 次插值多项式:

)

())(()()

()()(1101100

n 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 ,,)()(0

1∈-=∏=+ξω，且依赖于.x 。

3.2.2 三次样条插值

所谓三次样条插值多项式)(x S n 是一种分段函数，它在节点

i x 011()n n a x x x x b -=<

22331111111()[()()]()()666[,]1,2,,.

i i i i i i i i i i i i i i i

i 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 的计算方法如下：令：

111111116()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 μλ-+++==???-，

对于第一种边界条件下有

???

+-=+---00011

0111

)'],([62]),['(62h f x x M M h x x f f M M n n n n n n

如果令,])

,['(6,1,)'],[(6,11

1000100---==-==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 11011011

μλμλμλ

3.2.3 曲线拟合

由已知的离散数据点选择与实验点误差最小的曲线

)(...)()()(1100x a x a x a x S n n ???+++=

称为曲线拟合的最小二乘法。若记

),()()(),(0

i k i j m

i i k j x x x ??ω??∑==

k i k i m

i i k d x x f x f ≡=∑=)()()(),(0

?ω?

上式可改写为),...,1,0(;),(n k d a k j n

j j k -=∑=??这个方程成为法方程，可写成距阵形式

d Ga =

其中,),...,,(,),...,,(1010T n T n d d d d a a a a ==

????

????????=),(),(),()(),(),(),(),(),(10

1110101000n n n n n n G ?????????????????? 。

它的平方误差为：.)]()([)(||||20

22i i m

i i x f x S x -=

∑=ωδ

3.3 结果比较

图 3.1 三次样条插值

图3.2 拉格朗日插值

图3.1 最小二乘拟合曲线

问题四欧拉法解微分方程

4.1 计算题目

考虑微分方程(i) y x y 2'=,(ii)y x y )1(2'+=, 给定初始条件y(0)=1， (1)求方程的准确解；

(2)在区间[0,1]上，分别取步长为0.1,0.005,0.025，用Euler 方法求出数值解，并作图与准确解进行比较；

(3)在区间[0,1]上，分别取步长为0.1,0.005,0.025，用改进的Euler 方法求出数值解，并作图与准确解进行比较；

4.2.1 方程的准确解

方程的准确解为：(i)33

1x e

y =，(ii)y=2

)1(1+x e e

4.2.2 Euler 方法求解

Euler 单步法求解微分方程的迭代格式为

),(1n n n n y x hf y y +=+ n=0,1，...相应的Matlab 程序在附录程序4.

从对比图可以看出步长取的越小，实际值与近似值的相对误差越小。

图4.2.1 步长为0.1 图4.2.2 步长为0.05

图4.2.3 步长为0.01 图4.2.4 步长0.1

图4.2.5 步长0.05 图4.2 .6 步长0.05

结论：从图片来看，当步长取的越小，计算结果就越接近实际值。

4.2.3改进欧拉方法

如果在计算中，先用Euler 公式求得一个初步的近似值yn+1（称为预测值），再用梯形公式将它校正一次，就称为预测-校正公式，即改进的Euler 方法，其迭代格式为预测：

1),(n n n n y x hf y y +=+

校正：

)](),([2

1,11+-++++=n n n n n n y x f y x f h

y y

微分方程一：

图4.2.7 步长0.05 图4.2.8步长0.1

图4.2.9 步长0.01

微分方程二：

图4.2.10 步长取0.05 图4.2.11步长取0.1

图4.2.12 步长取0.01

问题五四阶龙格-库塔计算常微分方程初值问题

5.1 计算题目

给定初值问题：

≤

)0(

502

用经典的四阶Runge-Kutta 方法求解，步长分别取h=0.1,0.025，0.01计算并打印x=0.1i(i = 0,1,.....,.10)各点的值，并与准确值25031)(x e x y x +=-做比较。

5.2 四阶龙格-库塔方法分析

由拉格朗日微分中值定理

))(,()())(()()(1'1ξξξy hf x y x x y x y x y n n n n n +=-+=++

记))(,(*ξξy hf k =，则得到*1)()(k x y x y n n +=+

这样，只要给出了k 的一种算法，就得到求解微分方程初值问题的一种计算公式。

四阶Runge-Kutta 法就是用4321k k k k 的加权平均值来近似k*，相应的方法称为标准四阶Runge-Kutta 法，最常见为如下公式[1]

)

,()

,21()

,()

22(6

342312143211k y h x hf k k y h x hf k k y h x hf k y x hf k k k k k y y n n n n n n n n n n ++=++=++==++++=+

5.3 程序流程图

是

否

5.4 标准四阶Runge-Kutta 法Matlab 实现

function[x,y]=Runge_Kutta(ydot_fun,x0,y0,h,N) %ydot_fun 为一阶微分方程的函数； % x0,y0为初始条件 %h 为区间步长 %N 为区间的个数 %x 为Xn 构成的向量 %y 为Yn 构成的向量 x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0;y(1)=y0; for n=1:N

输入微分方程，求解的自变量范围，初值，步长等

信息

确定K 的点数是否

到4

)

,()

,21()

,()

22(6

342312143211k y h x hf k k y h x hf k k y h x hf k y x hf k k k k k y y n n n n n n n n n n ++=++=++==++++=+

结束

x(n+1)=x(n)+h;

k1=h*feval(ydot_fun,x(n),y(n));

k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1);

k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2);

k4=h*feval(ydot_fun,x(n)+h,y(n)+k3);

y(n+1)=y(n)+1/6*(k1+2k2+2k3+k4);

End

计算题目初值程序在附录1程序清单5。

5.5 计算结果及比较

运行程序程序附录1-5，并计算了100个值。步长分别取h=0.1,0.025，0.01打印x=0.1i(i = 0,1,.....,.10)各点的值，如表5-1,：