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

2014数值分析上机实习报告

姓名：吴朋朋

学号：2014200328

专业：车辆工程

联系电话：[1**********]

序言

数值分析方法在工程技术领域中的应用越来越广泛，作为连接工程实际与计算机运算处理的基础科学也越来越受到人们的重视，它能够使复杂的工程计算较好地利用计算机进行求解分析. 并得出结论。本上机实习报告以所学的《数值分析》课程为基础，运用MATLAB 编写计算程序得出计算结果，并对结果进行了一些必要的分析。

Matlab 是矩阵实验室(Matrix Laboratory)的简称，是由Mathworks 公司研发推出的一套高性能的数值计算和可视化软件。经过二十多年的不断发展和完善，现已发展成为由开发环境、数学函数库、Matlab 语言、图形处理系统和应用程序接口这五大部分组成的强大数学应用软件，广泛应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。

Matlab 是一个写程序与执行命令同步的交互式系统，其基本的数据元素没有维数的限制，与其他非交互式的编程语言（如C 语言）相比有其独特的特点和优势：

(1)工作平台和编程环境简单，能及时地报告出现的错误及进行出错原因分析；

(2)编程语言简单实用，语言可移植性好、可拓展性极强；

(3)强大的科学计算与处理能力，拥有庞大的数学、统计及工程函数库，且具有运算速度快和结果可靠性高的特点；

(4)多层面的图形图像处理系统，可以实现2D 与3D 数据图示、图像处理、动画生成、图形显示等功能以及开发GUI 应用程序等；

(5)提供了许多专业领域的丰富而高效的Matlab 工具箱；

(6)提供了便利的API 应用程序接口，用户可以编写C 或C++语言程序实现与Matlab 的交互。

因此本实习报告选用了MATLAB 做为编程计算工具进行求解分析。

上机题目及解答 . ...................................................................................................................... 1

1. 插值多项式和多项式拟合 . ........................................................................................... 1

1.1 题目 . ................................................................................................................... 1

1.2 N 次多项式拟合 .............................................................................................. 1

①当N=3时 ...................................................................................................... 1

②当N=4时 ...................................................................................................... 2

③当N=5时 ...................................................................................................... 3

④当N=6时 ...................................................................................................... 4

⑤小结 . .............................................................................................................. 5

1.3 牛顿插值法求解 . ............................................................................................... 5

①调用求解 . ...................................................................................................... 5

②验证评估 . ...................................................................................................... 6

③小结 . .............................................................................................................. 7

2. 雅格比迭代法、高斯—赛德尔迭代法 . ............................................................................... 8

2.1 题目 . ................................................................................................................... 8

2.2 小题(1)的求解分析 ........................................................................................... 8

①Jacobi 法求解Ax1=b1 .................................................................................. 8

②Guass-Seidel 法求解Ax1=b1 . ....................................................................... 8

③Jacobi 法求解Ax2=b2 .................................................................................. 9

④Guass-Seidel 法求解Ax2=b2 . ....................................................................... 9

2.2 小题(2)的求解分析 ......................................................................................... 10

①Jacobi 法求解Ax1=b1 ................................................................................ 10

②Guass-Seidel 法求解Ax1=b1 . ..................................................................... 10

③Jacobi 法求解Ax2=b2 ................................................................................ 11

④Guass-Seidel 法求解Ax2=b2 . ..................................................................... 11

2.2 小题(3)的求解分析 ......................................................................................... 12

①Jacobi 法求解Ax1=b1 ................................................................................ 12

②Guass-Seidel 法求解Ax1=b1 . ..................................................................... 12

3. 三次样条插值 . ........................................................................................................... 13

3.1题目 . .................................................................................................................. 13

3.2三次样条插值多项式的求解 . .......................................................................... 13

3.3比较分析 . .......................................................................................................... 14

总结与体会 . ............................................................................................................................ 15

附录程序清单 . ...................................................................................................................... 16

上机题目及解答

1. 插值多项式和多项式拟合

1.1 题目

某过程涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi 与yi 之间的对应数据如下，xi=1,2,…,10

yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795

103.5743 97.4847 78.2392

(1)请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。

(2)请用插值多项式给出最好近似结果

下列数据为另外的对照记录，它们可以作为近似函数的评价参考数据。

xi =

Columns 1 through 7

1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000

Columns 8 through 14

4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000

Columns 15 through 17

7.1000 7.5000 7.9000

yi =

Columns 1 through 7

42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006

Columns 8 through 14

-18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840

Columns 15 through 17

79.5688 93.7700 102.3677

1.2 N 次多项式拟合

用次数分别为3，4，5，6的多项式进行拟合并验证。

当N=3，4，5，6时，分别调用自定义函数nihe ，并代入给定的x 对照数据进行验证，y 的值与对照记录的误差越小说明拟合结果越精确，执行结果如下：

①当N=3时

>> N=3;

X=1:10;

Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392];

H=nihe(X,Y,N) %调用

H =

- ([**************]3*x^3)/[**************]8 +

([**************]9*x^2)/[**************] -

([**************]9*x)/[1**********]664 + [**************]9/[1**********]416

>> vpa(H) %转化为小数形式

ans =

- 1.[***********][1**********]45*x^3 +

19.[***********][1**********]7*x^2 -

94.[***********][1**********]5*x +

131.[***********][1**********]

>> x=[1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000

5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000];

y=eval(H);

y0=[42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 93.7700 102.3677];

dif=abs(y0-y) %验证

dif =

Columns 1 through 9

12.0573 26.4646 30.9125 27.0645 17.3272 4.5337 8.4418 19.0747 25.5260

Columns 10 through 17

26.8703 23.1790 15.4468 5.3713 4.9764 13.5426 18.6767 19.4555

②当N=4时

>> N=4;

X=1:10;

Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392];

H=nihe(X,Y,N)

H =

- ([**************]5*x^4)/[**************]84 +

([**************]5*x^3)/[**************]4 -

([**************]7*x^2)/[**************] +

([**************]3*x)/[1**********]416 +

[**************]5/[**************]6

>> vpa(H)

ans =

- 0.[***********][1**********]007*x^4 +

7.[***********][1**********]32*x^3 -

42.[***********][1**********]4*x^2 +

73.[***********][1**********]6*x +

0.[***********][1**********]549

>> x=[1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000

5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000];

y=eval(H);

y0=[42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 93.7700 102.3677];

dif=abs(y0-y)

dif =

Columns 1 through 9

2.9932 7.5801 9.2237 7.5935 3.3364 2.2393 7.5498 11.1262 11.9506

Columns 10 through 17

9.6842 4.7502 1.7394 8.2041 12.9249 14.4346 11.9037

5.4647

③当N=5时

>> N=5;

X=1:10;

Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392];

H=nihe(X,Y,N)

H =

([**************]1*x^5)/[**************]68 -

([**************]7*x^4)/[**************] +

([**************]5*x^3)/[**************] -

([1**********]391*x^2)/[1**********] + ([**************]1*x)/[1**********]416 - [**************]9/[1**********]08

>> vpa(H)

ans =

0.[***********][**************]*x^5 -

3.[***********][1**********]94*x^4 +

34.[***********][1**********]8*x^3 -

163.[***********][1**********]*x^2 +

304.[***********][1**********]*x -

139.[***********][1**********]

>> x=[1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000

5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000];

y=eval(H);

y0=[42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 93.7700 102.3677];

dif=abs(y0-y)

dif =

Columns 1 through 9

9.1435 6.6990 1.2208 3.5644 5.8723 5.4228 2.9349 0.3847 3.2987

Columns 10 through 17

4.8868 4.7502 3.0580 0.4478 2.1835 3.9498 4.2416

2.9288

④当N=6时

>> N=6;

X=1:10;

Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392];

H=nihe(X,Y,N)

H =

([**************]1*x^6)/[***********] -

([**************]7*x^5)/[**************]6 +

([**************]7*x^4)/[**************] -

([**************]5*x^3)/[**************] -

([**************]7*x^2)/[**************]2 +

([1**********]037*x)/[1**********]76 - [**************]1/[**************]

>> vpa(H)

ans =

0.[***********][**************]*x^6 -

0.[***********][1**********]029*x^5 +

5.[***********][1**********]18*x^4 -

16.[***********][1**********]7*x^3 -

0.[***********][1**********]041*x^2 +

66.[***********][1**********]5*x -

18.[***********][1**********]3

>> x=[1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000

5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000];

y=eval(H);

y0=[42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 93.7700 102.3677];

dif=abs(y0-y)

dif =

Columns 1 through 9

1.7365 1.0853 0.1701 1.1172 1.3124 0.7720 0.1623 1.0210 1.4061

Columns 10 through 17

1.1400 0.3263 0.6889 1.4448 1.5472 0.8526 0.4092

1.6311

⑤小结

由上述执行结果可知6次多项式拟合的误差最小，即近似结果最好。则可取

f(x)≈H 6=

0.[***********][**************]*x^6 -

0.[***********][1**********]029*x^5 +

5.[***********][1**********]18*x^4 -

16.[***********][1**********]7*x^3 -

0.[***********][1**********]041*x^2 +

66.[***********][1**********]5*x -

18.[***********][1**********]3

1.3 牛顿插值法求解

利用Newton 插值法求解近似多项式，调用自定义函数ndcz 及验证过程如下：

①调用求解

>> X=1:10;

Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392];

N=ndcz(X,Y)

N =

(57131*x)/10000 - ([**************]5*(x - 1)*(x - 2))/[**************] + (11771*(x - 1)*(x - 2)*(x - 3))/2500 + ([**************]*(x - 1)*(x - 2)*(x - 3)*(x - 4))/[**************]6 - ([**************]3*(x - 1)*(x - 2)*(x -

3)*(x - 4)*(x - 5))/[**************]2 + ([**************]5*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6))/[***********] + ([1**********]323*(x -

1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7))/[**************]84 - ([**************]5*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x -

7)*(x - 8))/[***********]8 + ([**************]9*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9))/[***********]64 + 289457/10000

>> N=simplify(N) %化简多项式

N =

([**************]9*x^9)/[***********]64 -

([***********]*x^8)/[***********]64 +

([***********]9*x^7)/[***********]32 -

([***********]97*x^6)/[***********]32 +

([***********]697*x^5)/[***********]64 -

([***********]2061*x^4)/[***********]64 +

([***********]81187*x^3)/[***********]000 -

([***********]432279*x^2)/[***********]60000 +

([***********]35309*x)/[***********]000 -

[***********]2371/[***********]000

>> vpa(N)

ans =

0.[***********][***********]*x^9 -

0.[***********][**************]3*x^8 +

0.[***********][**************]*x^7 -

0.[***********][1**********]976*x^6 +

5.[***********][1**********]67*x^5 -

22.[***********][1**********]8*x^4 +

50.[***********][1**********]7*x^3 -

86.[***********][1**********]*x^2 +

113.[***********][1**********]*x - 25.[***********]0874231791

②验证评估

>> x=[1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000

5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000];

y=eval(H);

y0=[42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 93.7700 102.3677];

dif=abs(y0-y)

dif =

Columns 1 through 9

1.7365 1.0853 0.1701 1.1172 1.3124 0.7720 0.1623 1.0210 1.4061

Columns 10 through 17

1.1400 0.3263 0.6889 1.4448 1.5472 0.8526 0.4092

1.6311

③小结

由牛顿插值法得到的近似多项式为

N =

0.[***********][***********]*x^9 -

0.[***********][**************]3*x^8 +

0.[***********][**************]*x^7 -

0.[***********][1**********]976*x^6 +

5.[***********][1**********]67*x^5 -

22.[***********][1**********]8*x^4 +

50.[***********][1**********]7*x^3 -

86.[***********][1**********]*x^2 +

113.[***********][1**********]*x - 25.[***********]0874231791

由上面的验证结果可知，Newton 插值法的误差与6次多项式拟合几乎相同。

2. 雅格比迭代法、高斯—赛德尔迭代法

2.1 题目

用雅格比法与高斯－赛德尔迭代法解下列方程组Ax =b1或Ax =b2，研究其收敛性。上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。

(1)A 行分别为A 1=[6,2,-1],A 2=[1,4,-2],A 3=[-3,1,4]； b 1=[-3,2,4]T ；b 2=[100,-200,345]T 。 (2) A 行分别为A 1=[1,0,8,0.8],A 2=[0.8,1,0.8],A 3=[0.8,0.8,1]；b 1=[3,2,1] T ； b 2=[5,0,-10]T 。 (3)A 行分别为A 1=[1,3],A 2=[-7,1]；b1=[4,6]T 。

2.2 小题(1)的求解分析

①Jacobi 法求解Ax1=b1

>> A=[6,2,-1;1,4,-2;-3,1,4]; b1=[-3;2;4]; x0=[0;0;0]; eps=1e-4;

x1=Jacobi(A,b1,x0,eps)

n = 16 T =

-0.5000 0.5000 1.0000 -0.5000 1.1250 0.5000 -0.7917 0.8750 0.3438 -0.7344 0.8698 0.1875 -0.7587 0.7773 0.2318 -0.7205 0.8056 0.2367 -0.7291 0.7984 0.2582 -0.7231 0.8114 0.2536 -0.7282 0.8076 0.2548 -0.7267 0.8095 0.2520 -0.7278 0.8077 0.2526 -0.7271 0.8083 0.2522 -0.7274 0.8079 0.2526 -0.7272 0.8081 0.2525 -0.7273 0.8080 0.2526 -0.7273 0.8081 0.2525 x1 = -0.7273 0.8081 0.2525

②Gauss-Seidel 法求解Ax1=b1

>> A=[6,2,-1;1,4,-2;-3,1,4]; b1=[-3;2;4]; x0=[0;0;0]; eps=1e-4;

x1=GaussSeidel(A,b1,x0,eps)

n = 10 T =

-0.5000 0.6250 0.4688 -0.6302 0.8919 0.3044 -0.7466 0.8388 0.2304 -0.7412 0.8005 0.2440 -0.7262 0.8035 0.2545 -0.7254 0.8086 0.2538 -0.7272 0.8087 0.2524 -0.7275 0.8081 0.2524 -0.7273 0.8080 0.2525 -0.7272 0.8081 0.2525 x1 = -0.7272 0.8081 0.2525

③Jacobi 法求解Ax2=b2

>> A=[6,2,-1;1,4,-2;-3,1,4]; b2=[100;-200;345]; x0=[0;0;0]; eps=1e-4;

x2=Jacobi(A,b2,x0,eps)

n = 23 T =

16.6667 -50.0000 86.2500 47.7083 -11.0417 111.2500 38.8889 -6.3021 124.7917 39.5660 2.6736 116.9922 35.2742 -1.3954 115.2561 36.3411 -1.1905 113.0545 35.9059 -2.5581 113.8035 36.4866 -2.0747 113.8189 36.3281 -2.2122 114.1336 36.4263 -2.0152 114.0491 36.3466 -2.0820 114.0735 36.3729 -2.0499 114.0304 36.3550 -2.0780 114.0422 36.3664 -2.0677 114.0358 36.3619 -2.0737 114.0417 36.3648 -2.0696 114.0398 36.3632 -2.0713 114.0410 36.3639 -2.0703 114.0402 36.3635 -2.0709 114.0405 36.3637 -2.0706 114.0403 36.3636 -2.0708 114.0404 36.3637 -2.0707 114.0404 36.3636 -2.0707 114.0404 x2 = 36.3636 -2.0707 114.0404

④Gauss-Seidel 法求解Ax2=b2

>> A=[6,2,-1;1,4,-2;-3,1,4]; b2=[100;-200;345]; x0=[0;0;0]; eps=1e-5;

x2=GaussSeidel(A,b2,x0,eps)

n = 18 T =

16.6667 -54.1667 112.2917 53.4375 -7.2135 128.1315 40.4264 3.9591 115.5800 34.6103 -0.8626 112.4234 35.6914 -2.7112 113.6964 36.5198 -2.2818 114.2103 36.4623 -2.0104 114.0993 36.3534 -2.0387 114.0247 36.3503 -2.0752 114.0316 36.3637 -2.0751 114.0415 36.3653 -2.0706 114.0416 36.3638 -2.0701 114.0404 36.3634 -2.0707 114.0402 36.3636 -2.0708 114.0404 36.3637 -2.0707 114.0404 36.3636 -2.0707 114.0404 36.3636 -2.0707 114.0404 36.3636 -2.0707 114.0404 x2 = 36.3636 -2.0707 114.0404

⑤小结

在本小题中，方程组的Jacobi 迭代和Gauss-Seidel 迭代均收敛，Gauss-Seidel

迭代的收敛速度要比Jacobi 迭代快很多。当右端项从b1增长到b2时，迭代次数增多，收敛速度减慢。

2.2 小题(2)的求解分析

①Jacobi 法求解Ax1=b1

②Guass-Seidel 法求解Ax1=b1

>> A=[1,0.8,0.8;0.8,1,0.8; 0.8,0.8,1]; b1=[3;2;1]; x0=[0;0;0]; eps=1e-4;

x1= GaussSeidel(A,b1,x0,eps) >> A=[1,0.8,0.8;0.8,1,0.8; 0.8,0.8,1]; b1=[3;2;1]; x0=[0;0;0]; eps=1e-4;

x1=Jacobi(A,b1,x0,eps)

方程组的Jacobi 迭代可能不收敛 n = 50 T =

1.0e+10 *

0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 …

……

0.0001 0.0001 0.0001 -0.0001 -0.0001 -0.0001 0.0002 0.0002 0.0002 -0.0003 -0.0003 -0.0003 0.0004 0.0004 0.0004

………

0.0180 0.0180 0.0180 -0.0288 -0.0288 -0.0288 0.0460 0.0460 0.0460 -0.0737 -0.0737 -0.0737 0.1179 0.1179 0.1179 -0.1886 -0.1886 -0.1886 0.3018 0.3018 0.3018 -0.4829 -0.4829 -0.4829 0.7726 0.7726 0.7726 -1.2361 -1.2361 -1.2361 x1 =

1.0e+10 * -1.2361 -1.2361 -1.2361

n = 31 T =

3.0000 -0.4000 -1.0800 4.1840 -0.4832 -1.9606 4.9551 -0.3955 -2.6476 5.4345 -0.2295 -3.1640 5.7148 -0.0407 -3.5393 5.8640 0.1403 -3.8034 5.9305 0.2983 -3.9831 5.9478 0.4282 -4.1008

…

5.7725 0.7728 -4.2362 5.7708 0.7724 -4.2345 5.7697 0.7718 -4.2332 5.7691 0.7713 -4.2323 5.7688 0.7708 -4.2317 5.7687 0.7704 -4.2313 5.7687 0.7700 -4.2310 5.7688 0.7698 -4.2309 5.7688 0.7696 -4.2308 5.7689 0.7695 -4.2307 5.7690 0.7694 -4.2307 5.7691 0.7693 -4.2307 x1 = 5.7691 0.7693 -4.2307

……

③Jacobi 法求解Ax2=b2

>> A=[1,0.8,0.8;0.8,1,0.8; 0.8,0.8,1]; b2=[5;0;-10]; x0=[0;0;0]; eps=1e-4;

x2=Jacobi(A,b2,x0,eps)

方程组的Jacobi 迭代可能不收敛 n = 50 T =

1.0e+10 *

0.0000 0 -0.0000 0.0000 0.0000 -0.0000

④Guass-Seidel 法求解Ax2=b2

>> A=[1,0.8,0.8;0.8,1,0.8; 0.8,0.8,1]; b2=[5;0;-10]; x0=[0;0;0]; eps=1e-4;

x2= GaussSeidel(A,b2,x0,eps)

n = 38 T =

5.0000 -4.0000 -10.8000 16.8400 -4.8320 -19.6064 24.5507 -3.9555 -26.4762 29.3453 -2.2953 -31.6400 32.1483 -0.4066 -35.3933 33.6399 1.4027 -38.0341 34.3051 2.9832 -39.8307 34.4780 4.2822 -41.0081 34.3808 5.3019 -41.7461 34.1554 6.0726 -42.1824 33.8878 6.6356 -42.4188

…

……

-0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0001 -0.0001 -0.0001 0.0001 0.0001 0.0001 -0.0001 -0.0001 -0.0001 0.0002 0.0002 0.0002 -0.0003 -0.0003 -0.0003 0.0006 0.0006 0.0006 -0.0009 -0.0009 -0.0009

…

……

32.6906 7.6931 -42.3070 32.6911 7.6927 -42.3070 32.6915 7.6925 -42.3071 32.6917 7.6923 -42.3072 32.6919 7.6922 -42.3073 32.6921 7.6922 -42.3074 32.6922 7.6922 -42.3075 32.6922 7.6922 -42.3076 x2 = 32.6922 7.6922 -42.3076

…

……

0.1572 0.1572 0.1572 -0.2515 -0.2515 -0.2515 0.4024 0.4024 0.4024 -0.6438 -0.6438 -0.6438 1.0301 1.0301 1.0301 x2 =

1.0e+10 * 1.0301 1.0301 1.0301

⑤小结

在本小题中，对于Ax =b1和Ax =b2，它们的Jacobi 迭代在50次以内均未收敛，而Gauss-Seidel 迭代收敛，且当右端项从b1变成b2后，迭代次数有所增加，收敛速度减慢。

2.2 小题(3)的求解分析

①Jacobi 法求解Ax1=b1

>> A=[1,3;-7,1]; b1=[4;6]; x0=[0;0]; eps=1e-4;

x1=Jacobi(A,b1,x0,eps)

方程组的Jacobi 迭代可能不收敛 n = 50 T =

1.0e+33 *

0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 …

②Guass-Seidel 法求解Ax1=b1

>> A=[1,3;-7,1]; b1=[4;6]; x0=[0;0]; eps=1e-4;

x1= GaussSeidel(A,b1,x0,eps)

方程组的GaussSeidel 迭代可能不收敛！ n = 50 T =

1.0e+66 *

0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000

…

-0.0000 0.0000 -0.0000 -0.0000 0.0001 -0.0002 0.0006 0.0005 -0.0016 0.0040 -0.0119 -0.0115 0.0344 -0.0836 0.2509 0.2410 -0.7231 1.7561 x1 =

1.0e+33 * -0.7231 1.7561

-0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0002 -0.0006 -0.0045 0.0136 0.0950 -0.2851 -1.9954 x1 =

1.0e+66 * -0.2851 -1.9954

③小结

在本小题中，方程的Jacobi 迭代和Gauss-Seidel 迭代在50次以内均未收敛，且迭代结果的绝对值均越来越大。

3. 三次样条插值

3.1题目

给定函数f (x ) =

, -5≤x ≤5，及节点x i =-5+i , i =0,1, 10，求其三1+x 2

次样条插值多项式（可取I 型或II 型边界条件），并画图及与f (x ) 的图形进行比较分析。

注：涉及到线性方程组求解问题需采用适当的求解算法。 3.2三次样条插值多项式的求解

调用spline3生成三弯矩三次样条多项式并绘图，调用格式及结果如下：

>> syms x fx=1/(1+x^2); xx=-5:5;

y=spline3(xx,fx)

S1= 0.00065819 x^3 + 0.014784 x^2 + 0.11326 x + 0.31747 ( -5

3.3比较分析

在执行完上一步的程序后，在同一坐标系下产生f (x ) 的图像，输入代码及图形结果如下：

>> hold on x=-5:0.02:5; f0=eval(fx); plot(x,f0,'r-.')

由图可以看出实验求得的分段三次样条插值函数的图像与原函数f (x ) 的图像较为接近，可见插值结果比较准确。

总结与体会

本报告是由本人严格按照实习要求独立完成。在给定的四道上机题中，第一题（插值多项式和多项式拟合）和第二题（雅格比迭代法、高斯—赛德尔迭代法求解方程组）为必做题，后两道为选做题，我选择的是第三题（对给定函数及节点求三次样条插值多项式）。我采用Matlab 完成程序的编写，从中深深感受到Matlab 与其它编程语言相比具有非常巨大的优势。

本次实习总共经历了六个步骤，分别是：选题；学习相关知识（包括Matlab 编程语言的学习及数值分析相关算法的回顾）；分析解题思路；编写程序代码；验证实验数据；编写实习报告。通过使用Matlab 编写数值计算的程序，我熟练掌握了Matlab 软件的使用方法，同时也使对各种算法的原理有了更深一层的理解。

通过上机实习，我学到了数值分析的相关知识，这种理论与实践相结合的方法，不仅提高了计算与编程能力，还增添了我对这门课程的学习兴趣。此外，我还认识到数值分析作为工程计算和科学计算连接计算机的一门基础学科的重要性，作为一名工科研究生，掌握这些基本技能必能为将来的工作和进一步的学习深造带来优势。

附录程序清单

1. 多项式拟合nihe.m

%求变量x,y 对于数据集X,Y 的N 次多项式拟合 function H=nihe(X,Y,N) m=length(X); A=ones(m,2*N+1); for i=1:2*N A(:,i+1)=X.^i; end ff=sum(A); G=zeros(N+1); for j=1:(N+1) g=ff(j:N+j); G(j,:)=g; end

fy=zeros(m,N+1); for k=1:N+1

fy(:,k)=X.^(k-1).*Y; end

b=sum(fy)'; a=G\b;

aa=a(end:-1:1); H=poly2sym(aa);

2. 牛顿插值法ndcz.m

function N=ndcz(X,Y) n=length(X); A=zeros(n); A(:,1)=Y'; w=1; N=Y(1); syms x ; for j=2:n for i=j:n

A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end

w=w*(x-X(j-1)); N=N+A(j,j)*w; end

3. 雅格比法Jacobi.m

function x=Jacobi(A,b,x0,eps) D=diag(diag(A)); L=tril(A,-1); U=triu(A,1); D=inv(D); B=-D*(L+U); f=D*b; n=0; for i=1:50 n=n+1; x=B*x0+f; T(i,:)=x; if norm(x-x0)>=eps x0=x; elseif norm(x-x0)

disp(' 方程组的Jacobi 迭代可能不收敛！' ) end end if i==50

disp(' 方程组的Jacobi 迭代可能不收敛！' ) end

n %n为迭代次数

T %T为各次迭代结果构成的矩阵

4. 高斯－赛德尔迭代法

%高斯-赛德尔法求解Ax=b

function x=GaussSeidel(A,b,x0,eps) D=diag(diag(A)); L=tril(A,-1); U=triu(A,1); C=inv(D+L); B=-C*U; f=C*b; n=0; for i=1:50 n=n+1; x=B*x0+f; T(i,:)=x;

if norm(x-x0)>=eps x0=x; elseif norm(x-x0)

break ;

end

if i==50

disp(' 方程组的GaussSeidel 迭代可能不收敛！' )

end

n %n为迭代次数

T %T为各次迭代结果构成的矩阵

5. 三次样条插值spline3.m

function y=spline3(xx,fx) % fx为给定函数，自变量为x ，xx 为指定节点集 gfx=diff(fx);

g2fx=diff(gfx);

x=xx;

Fx=eval(fx);

Gfx=eval(gfx);

G2fx=eval(g2fx);

n=length(xx)-1;

lamda(1)=1;

for ii=1:n

h(ii)=xx(ii+1)-xx(ii);

end

for i=2:n

lamda(i)=h(i)/(h(i-1)+h(i));

end

for j=1:n-1

miu(j)=1-lamda(j+1);

end

miu(n)=1;

A=2*eye(n+1)+diag(lamda,1)+diag(miu,-1);

d(1)=6/h(1)*((Fx(2)-Fx(1))/h(1)-Gfx(1));

d(n+1)=6/h(n)*(Gfx(n+1)-(Fx(n+1)-Fx(n))/h(n));

for k=2:n

d(k)=6/(h(k)+h(k-1))*((Fx(k+1)-Fx(k))/h(k)-(Fx(k)-Fx(k-1))/h(k-1)); end

d=d';

x0=zeros(n+1,1);

eps=1e-4;

M=GaussSeidel(A,d,x0,eps); %调用高斯赛德尔函数求解AM=d

syms x ;

for m=1:n

Sx(m)=M(m)*(xx(m+1)-x)^3/(6*h(m))+M(m+1)*(x-xx(m))^3/(6*h(m))+(Fx(m)-M(m)/6*h(m)^2)*(xx(m+1)-x)/h(m)+(Fx(m+1)-M(m+1)/6*h(m)^2)*(x-xx(m))/h(m);

end

y=Sx; % y的各元素依次对应相应区间的三次样条多项式

%输出各区间函数表达式及三次样条函数图像

for I=1:n

x=xx(I):0.02:xx(I+1);

a1=poly2str(sym2poly(y(I)),'x' );

fprintf('S%d=%s ( %d

ff=eval(y(I));

plot(x,ff)

if I

hold on