西南交大数值分析上机报告.docx

西南交大数值分析上机报告.docx

《西南交大数值分析上机报告.docx》由会员分享，可在线阅读，更多相关《西南交大数值分析上机报告.docx（22页珍藏版）》请在冰豆网上搜索。

西南交大数值分析上机报告

数值分析实验报告

班级：

桥梁二班

姓名：

学号：

电话：

2013年12月

序言

本次数值分析实验采用MATLAB语言。

MATLAB是当今最流行的科学计算语言之一，与FORTRAN、C语言一样，在科学计算领域被广泛的应用。

本次实验选择MATLAB语言，首先因为MATLAB具有强大的科学计算及数据处理能力，尤其在数值计算与优化方面，具有无可比拟的优越性。

MATLAB拥有600多个工程中要用到的数学运算函数，可以方便地实现用户所需的各种计算功能。

函数中所使用的算法都是科研和工程计算中的最新研究成果，而且经过了各种优化及容错处理，因此使用起来稳定性和可靠性非常高，在通常情况下，可以用它来代替底层编程语言，如C和C++等。

在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。

MATLAB函数所能解决的问题包括矩阵运算、多维数组操作（阵列运算）、复数的各种运算、三角函数和其他初等数学函数运算、非线性方程求根、线性方程组的求解、微分方程及偏微分方程组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、建模和动态仿真等等。

其次 MATLAB具有出色的图形处理功能。

在科学计算中，往往需要用各种图形把数值计算的结果形象地表现出来，以帮助人们更好地理解、认识和发现其中的科学规律。

MATLAB不仅提供数值计算功能和符号运算功能，而且自诞生之日起就具有方便的数据可视化功能，使计算结果的可视化要求得到充分满足。

MATLAB在二维曲线和三维曲面的绘制和处理等方面的功能比一般数据可视化软件更加完善。

在MATLAB中可以对图形对象的属性进行编辑，以提高输出图形的效果。

对一些特殊的可视化要求，例如图形动画等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。

MATLAB最突出的优点是程序语言简单易用。

早期用于科学计算的计算机语言，由于计算机内存容量和运算速度的限制等原因，常常要定义常量、变量、向量和矩阵等的不同的数据类型，结果导致编程过于复杂化。

和这些语言不一样，MATLAB语言对他们进行了高度抽象，实现了数据类型的高度统一，即常量、变量、向量和矩阵等都具有相同的数据类型。

因此，用户不需要事先分别定义常量、变量、向量和矩阵等的数据类型就可以直接使用他们。

MATLAB语言是一种“数学形式的语言”，它的操作和功能函数指令就是用平时计算机和数学书上的英文单词和符号来表达的，比BASIC、FORTRAN和C等语言更接近于人们书写的数学计算公式、更接近于人们进行科学计算的思维方式，因此，MATLAB语言简单自然，学习和使用更容易。

MATLAB程序文件是一个纯文本文件，扩展名为.m，用任何字处理软件都可以对它进行编辑。

MATLAB本身就像一个解释系统，对其中的函数程序的执行以一种解释执行的方式进行，程序不必经过编译就可以直接运行，而且能够及时报告出现的错误，进行出错原因分析。

因此，程序调试容易、编程效率高。

目录

一、解线性方程组1

1.计算与分析1

二、微分方程初值问题2

1.计算与分析2

三、非线性方程根的求解3

1.分析与计算3

四、数值积分4

1.计算与分析4

总结6

附录7

一、Gauss消元法程序7

二、Gauss-Seidel迭代法程序7

三、Runge-Kutta法程序8

四、Newton迭代法程序9

五、复合梯形公式程序9

六、复合Simpson公式程序10

七、Romberg算法程序10

一、

解线性方程组

写出对一般的线性方程组通用的Gauss消元,Gauss-Seidel迭代程序。

并以下面的线性方程组为例进行计算，讨论所得到的计算结果是否与理论一致。

（1）

（2）

（3）

1.计算与分析

（1）利用选列主元的Gauss消元法编写MATLAB程序，解上述5个方程组。

结果如表1.1所示：

表1.1Gauss消元法计算结果

x

（1）

x

（2）

x（3）

x（4）

x（5）

x1

-0.72727273

36.36363636

5.76923077

32.69230769

-0.63636364

x2

0.80808081

-2.07070707

0.76923077

7.69230769

1.54545455

x3

0.25252525

114.04040404

-4.23076923

-42.30769231

Gauss消元法为直接解法，从消元的过程可知其结果为精确解，误差来源于机器误差，即计算机内部的存储精度。

MATLAB用二进制双精度浮点格式来表达一个数，其计算结果至少可以有12位十进制有效数字的可信度。

（2）利用Gauss-Seidel法编写MATLAB程序，用||x（k）-x（k-1）||<ε作为控制迭代的终止条件，ε取10-6。

方程组

（1）系数矩阵A1为行对角占优阵，故Gauss-Seidel法收敛。

方程组

（2）系数矩阵A2为对称正定阵，故Gauss-Seidel法收敛。

方程组（3）的Gauss-Seidel迭代矩阵为

，其特征值为λ1=0，λ2=-21，ρ（B）=21>1，故Gauss-Seidel迭代发散。

将方程组的两个方程顺序对换一下，变为

，则其系数矩阵为行对角占优阵，Gauss-Seidel迭代收敛。

Gauss-Seidel迭代法计算结果如表1.2所示：

表1.2Gauss-Seidel迭代计算结果

x1

x2

x3

迭代次数

x10

1

1

1

15

x1k

-0.72727289

0.80808085

0.25252512

x20

1

1

1

20

x2k

36.36363650

-2.07070701

114.04040413

x30

1

1

1

45

x3k

5.76923210

0.76922890

-4.23076879

x40

1

1

1

51

x4k

32.69230977

7.69230745

-42.30769378

x50

1

1

1

6

x5k

-0.63636362

1.54545454

二、微分方程初值问题

给定初值问题

（精确解为

），用Runge-Kutta4阶算法按步长

求解，分析其中遇到的现象及问题。

1.计算与分析

编写4级4阶的Runge-Kutta法MATLAB程序。

取步长0.2，计算结果列于表2.1中。

表中y（xk）是当x=xk时的准确值，yk是由Runge-Kutta法得到的当x=xk时的值。

表2.1Runge-Kutta法数值计算结果

步长0.2

0

1

1

0.2

5

1.83156E-02

0.4

25

3.35463E-04

0.6

125

6.14421E-06

0.8

625

1.12535E-07

1.0

3125

2.06115E-09

1.2

15625

3.77513E-11

1.4

78125

6.91440E-13

1.6

390625

1.26642E-14

1.8

1953125

2.31952E-16

2.0

9765625

4.24835E-18

从表中可以看出步长取0.2时Runge-Kutta法是发散的。

4级4阶Runge-Kutta法的绝对稳定区间为（-2.78，0），对于实验方程，λh=-20h，当h=0.2时λh=-4，故Runge-Kutta法无法收敛。

步长取0.1，0.05时，λh分别为-2，-1，位于绝对稳定区间，预计Runge-Kutta法收敛。

表2.2列出了步长为0.1和0.05时Runge-Kutta法的计算结果和误差

表2.2Runge-Kutta法数值计算结果

步长0.1

步长0.05

0

1

1

0

0

1

1

0

0.1

3.33333E-01

1.35335E-01

-1.97998E-01

0.05

3.75000E-01

3.67879E-01

-7.12056E-03

0.2

1.11111E-01

1.83156E-02

-9.27955E-02

0.10

1.40625E-01

1.35335E-01

-5.28972E-03

0.3

3.70370E-02

2.47875E-03

-3.45583E-02

0.15

5.27344E-02

4.97871E-02

-2.94731E-03

0.4

1.23457E-02

3.35463E-04

-1.20102E-02

0.20

1.97754E-02

1.83156E-02

-1.45975E-03

0.5

4.11523E-03

4.53999E-05

-4.06983E-03

0.25

7.41577E-03

6.73795E-03

-6.77824E-04

0.6

1.37174E-03

6.14421E-06

-1.36560E-03

0.30

2.78091E-03

2.47875E-03

-3.02162E-04

0.7

4.57247E-04

8.31529E-07

-4.56416E-04

0.35

1.04284E-03

9.11882E-04

-1.30961E-04

0.8

1.52416E-04

1.12535E-07

-1.52303E-04

0.40

3.91066E-04

3.35463E-04

-5.56034E-05

0.9

5.08053E-05

1.52300E-08

-5.07900E-05

0.45

1.46650E-04

1.23410E-04

-2.32400E-05

1.0

1.69351E-05

2.06115E-09

-1.69330E-05

0.50

5.49937E-05

4.53999E-05

-9.59374E-06

1.1

5.64503E-06

2.78947E-10

-5.64475E-06

0.55

2.06226E-05

1.67017E-05

-3.92092E-06

1.2

1.88168E-06

3.77513E-11

-1.88164E-06

0.60

7.73348E-06

6.14421E-06

-1.58927E-06

1.3

6.27225E-07

5.10909E-12

-6.27220E-07

0.65

2.90006E-06

2.26033E-06

-6.39727E-07

1.4

2.09075E-07

6.91440E-13

-2.09074E-07

0.70

1.08752E-06

8.31529E-07

-2.55993E-07

1.5

6.96917E-08

9.35762E-14

-6.96916E-08

0.75

4.07820E-07

3.05902E-07

-1.01918E-07

1.6

2.32306E-08

1.26642E-14

-2.32306E-08

0.80

1.52933E-07

1.12535E-07

-4.03975E-08

1.7

7.74352E-09

1.71391E-15

-7.74352E-09

0.85

5.73498E-08

4.13994E-08

-1.59504E-08

1.8

2.58117E-09

2.31952E-16

-2.58117E-09

0.90

2.15062E-08

1.52300E-08

-6.27618E-09

1.9

8.60390E-10

3.13913E-17

-8.60390E-10

0.95

8.06481E-09

5.60280E-09

-2.46201E-09

2.0

2.86800E-10

4.24835E-18

-2.86800E-10

1.00

3.02430E-09

2.06115E-09

-9.63150E-10

1.05

1.13411E-09

7.58260E-10

-3.75854E-10

1.10

4.25293E-10

2.78947E-10

-1.46346E-10

1.15

1.59485E-10

1.02619E-10

-5.68660E-11

1.20

5.98068E-11

3.77513E-11

-2.20554E-11

1.25

2.24275E-11

1.38879E-11

-8.53960E-12

1.30

8.41033E-12

5.10909E-12

-3.30124E-12

1.35

3.15387E-12

1.87953E-12

-1.27434E-12

1.40

1.18270E-12

6.91440E-13

-4.91262E-13

1.45

4.43513E-13

2.54367E-13

-1.89147E-13

1.50

1.66318E-13

9.35762E-14

-7.27413E-14

1.55

6.23691E-14

3.44248E-14

-2.79443E-14

1.60

2.33884E-14

1.26642E-14

-1.07242E-14

1.65

8.77065E-15

4.65889E-15

-4.11176E-15

1.70

3.28899E-15

1.71391E-15

-1.57509E-15

1.75

1.23337E-15

6.30512E-16

-6.02861E-16

1.80

4.62515E-16

2.31952E-16

-2.30563E-16

1.85

1.73443E-16

8.53305E-17

-8.81126E-17

1.90

6.50411E-17

3.13913E-17

-3.36498E-17

1.95

2.43904E-17

1.15482E-17

-1.28422E-17

2.00

9.14640E-18

4.24835E-18

-4.89805E-18

从表中可以看出，当步长取0.1时，Runge-Kutta法是收敛的，但其最终全局误差尽管绝对在很小，相对误差却很大。

比真实值要大数个量级。

从微分方程的显式解（

）可以看出，该函数是急遽下降的，其导数，即斜率则是急遽增大的。

Runge-Kutta法可由泰勒方法推导而来，其思想即是用xk点处的斜率，xk+1点处的斜率及中点斜率的两个估计值的线性组合来代替

，利用公式

来预测

的值。

如果

变化很剧烈，而步长又取得较大，则可能在中间过程中会出现“大数吃小数”的现象，使结果精度降低。

当步长减半，取0.05时，从表中可以看出其计算精度要好很多，但相对误差仍较大。

对于浮点数而言，越靠近0，精度就越高；越远离0，精度就越低。

MATLAB中的浮点误差值eps为1到下一个能表示的比1大的浮点数之间的距离。

换言之，MATLAB是无法区分1与1+eps之间的数的，它们将被舍入到1或1+eps。

在MATLAB命令窗口键入eps可得到其值为2.220446049250313e-16。

由显示函数

计算的结果从1降到1e-18量级，小于eps，中间过程必然有舍入误差。

三、非线性方程根的求解

利用牛顿法求方程

于区间

的根，考虑不同初值下牛顿法的收敛情况。

1.分析与计算

令

，转化为求

在区间

上的零点。

通过观察方程易知零点在3附近。

取初值x0=3进行牛顿法迭代。

编写MATLAB程序，计算结果如表3.1所示，可见迭代收敛很快。

如果从区间端点开始迭代，初值分别取2和4，计算结果如表3.2所示，同样，迭代收敛很快。

表3.1Newton迭代法计算结果

k

0

3.00000000

-0.098612289

1

3.147918433

0.147918433

0.001177014

2

3.146193441

0.001724992

1.50195E-07

3

3.146193221

2.20177E-07

2.66454E-15

表3.2Newton迭代法计算结果

k

0

2.00000000

-0.69314718

4.00000000

0.613705639

1

3.386294361

1.386294361

0.166558146

3.181725815

0.818274185

0.024302056

2

3.149938394

0.236355967

0.002555499

3.146284845

0.03544097

6.25023E-05

3

3.146194257

0.003744137

7.06991E-07

3.146193221

9.16234E-05

4.24029E-10

4

3.146193221

1.03641E-06

5.41789E-14

对于实验方程可用MATLAB画出其在区间[2，4]上的图形，如图3.1所示。

可见其在区间内只有单根，对于单根，Newton法具有2阶收敛速度。

图3.1f（x）在[2，4]上图形

四、数值积分

已知

，利用复化梯形公式、复化Simpson公式和Romberg算法求

的近似值；并观察实际计算结果，比较它们的收敛速度。

1.计算与分析

（1）复合梯形公式与复合Simpson公式

分别编写复合梯形公式与复合Simpson公式的MATLAB程序，计算结果如表4.1所示。

表4.1复合梯形公式与复合Simpson公式计算结果

子区间数

n

误差

误差

1

3.00000000000000

0.14159265358979

3.133********333

0.00825932025646

2

3.10000000000000

.0415********

3.14156862745098

0.00002402613881

4

3.131********824

.010*********

3.14159250245871

0.00000015113109

8

3.138********109

0.00260415909870

3.14159265122482

0.00000000236497

16

3.14094161204139

0.00065104154840

3.14159265355284

0.00000000003696

可以看到，相对于复合梯形公式，复合Simpson公式的误差更快地收敛到0。

由于梯形公式的误差项为

阶，Simpson公式的误差项为

，因此Simpson公式比梯形公式具有更高的精度，所需计算的次也更少。

（2）Romberg公式

编写Romberg算法的MATLAB程序。

该方法需要逐次计算梯形公式值，Simpson公式和Cotes公式值，通过其线性组合来提高精度。

中间计算值存在矩阵中，如表4.2所示。

表4.2Romberg积分表

k

0

3.00000000000000

1

3.10000000000000

3.13333333333333

2

3.13117647058824

3.14156862745098

3.14211764705882

3

3.13898849449109

3.14159250245871

3.14159409412589

3.14158578376187

4

3.14094161204139

3.14159265122482

3.14159266114256

3.14159263839680

5

3.14142989317497

3.14159265355284

3.14159265370804

3.14159265359003

表4.3为Romberg积分的误差表。

表4.3Romberg积分误差表

k

0

0.14159265358979

1

0.04159265358979

0.00825932025646

2

0.01041618300156

0.00002402613881

0.00052499346903

3

0.00260415909870

0.00000015113109

0.00000144053610

0.00000686982792

4

0.00065104154840

0.00000000236497

0.00000000755277

0.00000001519300

5

0.00016276041482

0.00000000003696

0.00000000011824

0.00000000000024

可以看到Romberg积分表中第二序列正是复合Simpson公式计算结果。

同样容易验证，第三序列是复合Cotes公式计算结果。

第四序列（Romberg序列）则是n=8时的Newton-Cotes公式。

一般而言，右边的列收敛速度要快于左边的列，但并不是说所用公式的阶数越高，结果越精确。

因为高阶Newton-Cotes公式的

很大，导致数值稳定性较差。

从误差表中可以看出，同一行的Cotes序列并不比Simpson序列精度更高。

但对于相同子区间数，右边列比左边列精度要高。

Romberg序列收敛最快，同时精度也最好。

总结

本次数值分析实验采用了MATLAB程序设计语言，分别编写了Gauss消元法、Gauss-Seidel迭代法程序求解线性方程组；编写Newton迭代法程序求解非线性方程；编写Runge-Kutta法程序求解线性微分方程；编写复合梯形公式、复合Simpson法程序和Romberg公式程序求解数值积分问题。

对Gauss消元法采用了选列主元的方法，增加了程序的通用性；对Gauss-Seidel迭代法讨论了应用条件，要求其迭代矩阵满足一定的条件；对Newton迭代比较了不同初值的收敛速度；对Runge-Kutta法讨论了收敛的稳定区间，要求步长满足一定要求。

比较了复合梯形公式、复合Simpson公式和Romberg公式的收敛效率