数值分析上机报告

数值分析上机报告

班级：20级学隧2班 姓名：000000000 学号：00000000000

1

目录

1 序言 ................................................................................................................... 6 2 题目 ................................................................................................................... 7

2.1 题2 ......................................................................................................... 7

2.1.1 题目内容 ..................................................................................... 7 2.1.2 MATLAB程序 ................................................................................ 8 2.1.3 计算结果 ..................................................................................... 8 2.1.4 图形 ............................................................................................. 9 2.1.5 分析 ................................................................................................... 14 2.2 题3 ....................................................................................................... 14

2.2.1 题目内容 ................................................................................... 14 2.2.2 程序 ........................................................................................... 14 2.2.3 计算结果 ................................................................................... 14 2.2.4 图形 ........................................................................................... 15 2.2.5 分析 ........................................................................................... 16 2.3 选做题5 ............................................................................................... 16

2.3.1方法介绍 .................................................................................... 17 2.3.2计算结果及分析 ........................................................................ 17

3总结 .................................................................................................................. 18 4．附录 ............................................................................................................... 19

4.1 题1程序代码 ...................................................................................... 19 4.2 题2程序代码 ...................................................................................... 22 4.3 题3程序代码 ...................................................................................... 26

2

数值分析2015上机实习报告要求

1．应提交一份完整的实习报告。具体要求如下：

（1）报告要排版，美观漂亮（若是纸质要有封面，封面上）要标明姓名、

学号、专业和联系电话；

（2）要有序言，说明所用语言及简要优、特点，说明选用的考量； （3）要有目录，指明题目、程序、计算结果，图标和分析等内容所在位置，

作到信息简明而完全；

（4）要有总结，全方位总结机编程计算的心得体会；

（5）尽量使报告清晰明了，一般可将计算结果、图表及对比分析放在前面，

程序清单作为附录放在后面，程序中关键部分要有中文说明或标注，指明该部分的功能和作用。

2．程序需完好保存到期末考试后的一个星期，以便老师索取用于验证、询问或

质疑部分内容。

3．认真完成实验内容，可以达到既学习计算方法又提高计算能力的目的，还可

以切身体会书本内容之精妙所在，期间可以得到很多乐趣。 4．拷贝或抄袭他人结果是不良行为，将视为不合格。

5．请按任课老师要求的时间和载体（电子或纸质）提交给任课老师。

3

数值分析上机试题

（请在1-4题中选择两个题目，5-6中选择一个题目）

1． 分别用牛顿法，及基于牛顿算法下的Steffensen加速法

(1) 求ln(x+sinx)=0的根。初值x0分别取0.1, 1,1.5, 2, 4进行计算。 (2) 求sinx=0的根。初值x0分别取1，1.4，1.6, 1.8，3进行计算。 分析其中遇到的现象与问题。

2. 某过程测涉及两变量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

3．用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分

4

析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。

(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T, b2=[100,-200,345]T, (2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T, b2=[5,0,-10]T, (3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T, 4. 松弛因子对SOR法收敛速度的影响。

用SOR法求解方程组Ax=b,其中

??41??-3??????1?41??-2????-2?...??, B??? ???.?...????1?41???-2???-3?1?4?????要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代，记录近似解x(k)达到||x(k)-x(k-1)||<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当w?0或w?2会有什么影响？

5. 用Runge-Kutta 4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区

间的作用，推荐两种步长h=0.1,0.2。 注：此方程的精确解为：y=e-20x 6. 实验内容

(1) 实际验证梯形求积公式、Simpson求积公式、Newton-Cotes求积公式的代数精度。 (2) 针对上述三个函数和积分区间[a,b]，实验观察梯形求积公式、Simpson求积公式和

Newton-Cotes求积公式的复化求积公式的实际计算效果。 y=exp(-x.^2).*sin(10*x)+4; a=1; b=3; y=sin(5*x)./x.^3;a=2*pi;b=4*pi; y=sin(5*x)./x.^3;a=2*pi;b=9.4248;

5

1 序言

MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。

总的来说，该软件有三大特点。一是功能强大。具有数值计算和符号计算、计算结果和编程可视化、数学和文字统一处理、离线和在线计算等功能；二是界面友善、语言自然。MATLAB以复数处理作为计算单元，指令表达与标准教科书的数学表达式相近；三是开放性强。主要有以下优点：

1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；

2) 具有完备的图形处理功能，实现计算结果和编程的可视化；

3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；

4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。

MATLAB擅长于矩阵运算，在各种系统仿真方面应用广泛。他不同于普通的编程语言，集成有许多领域专家为各自领域开发的工具箱，直接调用即可。面向具体应用，使用MATLAB更有针对性。这种程序在数值分析方面应用广泛，基于它的众多优点，这次我应用MATLAB编程，把学习这个软件作为我数值分析课程的实践环节。

6

2 题目

2.1 题2

2.1.1 题目内容

2. 某过程测涉及两变量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

7

2.1.2 MATLAB程序

见附录4.1

2.1.3 计算结果 2.1.3.1 插值多项式

(1)线性插值误差 Error = 10.8098 (2)三次函数插值误差 Error = 7.2441

(3)分段三次Hermit插值误差 Error = 7.2441

(4)三次样条函数插值误差 Error = 1.0338 (5)最终结果

三次样条函数插值效果最好，误差最小。

2.1.3.2 多项式拟合

(1) 拟合次数 n = 3时 拟合误差 Error = 57.2759 3次拟合多项式

F(x) = -1.0326 x^3 + 19.3339 x^2 - 94.4787 x + 131.7944 (2) 拟合次数n = 4时 拟合误差 Error = 29.6408 4次拟合多项式

F(x) = -0.38185 x^4 + 7.368 x^3 - 42.1433 x^2 + 73.5334 x + 0.74498 (3) 拟合次数n = 5时 拟合误差

8

Error = 11.3278 5次拟合多项式

F(x) =0.098075 x^5 - 3.0789 x^4 + 34.502 x^3 - 163.5107 x^2 +

304.7282 x - 139.5019

(4)拟合次数 n = 6时 拟合误差 Error = 3.2509 6次拟合多项式

F(x) = 0.019359 x^6 - 0.54079 x^5 + 5.1137 x^4 - 16.8973 x^3 0.86696 x^2 + 66.375 x - 18.6991

(5)最终结果

n = 6时拟合效果最佳，误差error = 3.2509。 最佳拟合多项式为：

F(x) = 0.019359 x^6 - 0.54079 x^5 + 5.1137 x^4 - 16.8973 x^3 0.86696 x^2 + 66.375 x - 18.6991

2.1.4图形

2.1.4.1 插值多项式

9

––线性插值1201008060y40200-2012345x678910 图2-1 线性插值图像

三次函数插值1201008060y40200-2012345x678910 图2-2 三次函数插值图像

10

分段三次Hermite插值1201008060y40200-2012345x678910 图2-3 分段三次Hermit插值图像

三次样条函数插值1201008060y40200-2012345x678910 图2-4 三次样条函数插值图像

11

2.1.4.2 多项式拟合

3次拟合图像1201008060y40200-2012345x678910 图2-5次拟合图像

4次拟合图像1201008060y40200-2012345x678910 图2-6次拟合图像

12

5次拟合图像1201008060y40200-2012345x678910 图2-7次拟合图像

6次拟合图像1201008060y40200-2012345x678910 图2-8次拟合图像

13

2.1.5分析

由图2-1～2-4可以发现，插值曲线逐步变得光滑，插值误差越来越小。 由图2-5～2-8可以发现，随着拟合次数（2

2.2.1 题目内容

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

(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,

b2=[100,-200,345]T,

(2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,

b2=[5,0,-10]T,

(3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T,

2.2.2程序

见附录2

2.2.3计算结果

(1)雅格比法

取迭代初值矩阵=[0 0 0]，经过61步迭代得出稳定结果x1 = -0.7273 x2 = 0.8081 x3 = 0.2525

(2) 高斯－赛德尔迭代法

取迭代初值矩阵=[0 0 0]，经过36步迭代得出稳定结果x1 = -0.7273 x2

14

= 0.8081 x3 = 0.2525

2.2.4图形

雅格比迭代法收敛图像4.543.53相邻值误差2.521.510.500102030迭代步数40506070 图2-9 雅格比迭代法收敛图像

15

高斯－赛德尔迭代法收敛图像1.81.61.41.2相邻值误差10.80.60.40.2005101520迭代步数25303540 图2-10 高斯－赛德尔迭代法

2.2.5分析

运行程序后发现，高斯－赛德尔迭代法比起雅格比法要快得多，在取同样的初值条件下，高斯－赛德尔迭代法经过36步则迭代出稳定的结果；而雅格比法则需要61步。

另外，在右端项取b2的情况下，经过程序运行，发现经过1000次迭代仍然无法收敛，可见右端项对迭代收敛有很大影响。 2.3 选做题5

用Runge-Kutta 4阶算法对初值问题y’=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1、0.2。

注：此方程的精确解为：y=e-20x

16

2.3.1方法介绍

通过用一些点上函数值f(x ,y)的适当线性组合，来替代Euler法中的f(xk ,yk)，从而使方法的阶数更高，这就是Runge-Kutta法的基本思想。在本次上机实习中采用的是标准4阶Runge-Kutta法，其计算形式如下：

2.3.2计算结果及分析

利用FORTRAN语言编写程序program R_K（附录3），计算结果如下。

表2不同步长计算结果对比

y(x) x=0 x=0.1 x=0.2 x=0.3 h=0.1 1.000000 0.333333 0.111111 3.7037037E-02 1.2345679E-02 4.1152262E-03 1.3717421E-03 4.5724737E-h=0.2 1.000000 5.000000 准确解 1.00000000000000 0.13533528323661 0.01831563888873 0.00247875217667 x=0.4 25.000000 0.00033546262790 x=0.5 0.00004539992976 x=0.6 x=0.7 125.000000 17

0.00000614421235 0.00000083152872

04 x=0.8 1.5241579E-04 5.0805262E-05 1.6935088E-05 5.6450294E-06 1.8816764E-06 6.2722546E-07 2.0907515E-07 收敛 625.000000 0.00000011253517 x=0.9 0.00000001522998 x=1.0 3125.000000 0.00000000206115 x=1.1 0.00000000027895 x=1.2 15625.000000 0.00000000003775 x=1.3 0.00000000000511 x=1.4 结论 78125.000000 发散 0.00000000000069 从上面结果中可以看出取步长h=0.1时候结果可以很好的收敛，与准确解之间的误差也越来越小，但是取h=0.2时候结果就发生了发散。说明该方法对步长的选择比较敏感。

3总结

通过这次实《matlab》实验使我学习掌握了许多知识。首先是对matlab有了一个全新的认识，其次是对matlab的更多操作和命令的使用有了更高的掌握，最重要的事对matlab的处理能力有了一个更高的飞跃尤其是对相关函数的使用及相关问题的处理。

通过这次上机练习，让我对数值分析所介绍的迭代求解方法及其理论有了更深层次的理解，了解了各种方法之间的优缺点，并且认识到了自己在以前的学习中所存在的问题，及时的修补了自己知识上的漏洞。同时也提高了我在编写程序上的熟练程度，所以，我认为这次上机实习是非常有收获的，给予我学习数值分

18

析的帮助也是非常大的。

我也希望通过这篇总结来表达自己对知道老师的感谢之情，谢谢您的不懈努力和耐心指导，才使得我再这次的实验过程中收获的这么多，也正式您的不吝教诲才使得我们在这次实验中学习和收获了许多的有用的知识和技巧，我相信在以后的学习或者工作中一定有其用武之地。过多的感谢无以言表，万分感激，百口不胜言表，至此敬礼！

4．附录

4.1 题1程序代码 （1）插值多项式源程序 clc; clear; %给定数据对 x0 = [1:1:10];

y0 = [34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 103.5743 97.4847 78.2392]; x1 = [1.5:0.4:7.9];

ycheck = [42.1498 41.4620 35.1182 24.3852 11.2732 -12.3006 -18.1566 -17.9069 ...

-11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 102.3677]%检查点矩阵

y1 = interp1(x0,y0,x1,'linear'); %线性插值

y2 = interp1(x0,y0,x1,'spline'); %三次样条函数插值 y3 = interp1(x0,y0,x1,'pchip'); %分段三次Hermite插值 y4 = interp1(x0,y0,x1,'cubic'); %三次函数插值

%插值函数计算值与检算值差距计算

19

75.2795 -1.7813 93.7700 data1 = y1 - ycheck; data2 = y2 - ycheck; data3 = y3 - ycheck; data4 = y4 - ycheck;

%求2范数作为误差指标 error1 = norm(data1,2) error2 = norm(data2,2) error3 = norm(data3,2) error4 = norm(data4,2)

% 计算最小误差

error = [error1,error2,error3,error4]; min_error = error(1); for i = 1:3

if error(i+1) < min_error; min_error = error(i+1); end end min_error figure(1)

plot(x1,y1,'-b',x0,y0,'*r','markersize',5) grid on title('线性插值') figure(2)

plot(x1,y2,'-b',x0,y0,'*r','markersize',5) grid on

20

title('三次样条函数插值') figure(3)

plot(x1,y3,'-b',x0,y0,'.r','markersize',15) grid on

title('分段三次Hermite插值') figure(4)

plot(x1,y4,'-b',x0,y0,'.r','markersize',15) grid on

title('三次函数插值') （2）多项式拟合 clc; clear;

n = input('请输入拟合次数：n=') x0 = [1:1:10];

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

ycheck = [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]; %检查点矩阵

set(0,'defaultfigurecolor','w')

plot(x0,y0,'b',x0,y0,'r*','markersize',10),xlabel('x'),ylabel('y') %原始点图像 p = polyfit(x0,y0,n); %求拟合多项式的系数

fx1 = poly2str(p,'x'); %以显式方式写出拟合得到的f(x) x1 = [1.5:0.4:7.9]

21

y1 = polyval(p,x1) %由拟合多项式拟合得到的值

y2 = polyval(p,x0)

%用MATLAB库函数求2范数 data = (y2 - y0).*(y2 - y0); error = norm((y2-y0),2);

%编程求2范数 mistake1 = 0; for i = 1:10

mistake1 = mistake1 + abs(data(i)); end

mistake = sqrt(mistake1) ; figure(2)

set(0,'defaultfigurecolor','w')

plot(x1,y1,'-b',x0,y0,'.r','markersize',20) %绘出拟合图像 figure(3)

set(0,'defaultfigurecolor','w')

plot(x0,y2,'-b',x0,y0,'r--')%拟合图像与原图像的比较

error %拟合误差 mistake %拟合误差

fx = sprintf(fx1) %输出拟合得到的f(x)到屏幕上

4.2 题2程序代码 （1）雅可比法 clc;

22

clear;

orignal = input('请输入迭代初值矩阵：oiginal = ') n = input('请输入迭代次数：n = ')

m = input('请选择右端项(1-b1;2-b2):m = ') A = [6,2,-1;1,4,-2;-3,1,4];

b1 = [-3,2,4]'; %转置为列向量 b2 = [100,-200,345]';

d = diag(A); for i = 1:3 for j = 1:3 if i == j dj(i,j) = d(i); else dj(i,j) = 0; end end end dj

switch m case 1 t = b1; case 2 t = b2; end

B = eye(length(A)) - inv(dj)*A; g = inv(dj)*t;

x{1} = orignal';

%从矩阵A中提取对角元素 %写成对角阵 %选择右端项 %计算迭代矩阵B %计算矩阵g 23

for k = 1:n

X{1} = B*x{1}+g; x{1} = X{1};

diedai = A*X{1} %使用迭代得到的值计算左端项

error = diedai - b1 %左端项与右端项的差距，作为迭代误差指标

if error == 0

step = k; %记录迭代步数 break; end end step

x1 = X{1}(1) %方程组的根x1 x2 = X{1}(2) %方程组的根x2 x3 = X{1}(3) %方程组的根x3 （2）高斯－赛德尔迭代法 clc; clear;

orignal = input('请输入迭代初值矩阵：oiginal = ') n = input('请输入迭代次数：n = ') m = input('请选择右端项(1-b1;2-b2):m = ') A = [6,2,-1;1,4,-2;-3,1,4];

b1 = [-3,2,4]'; %转置为列向量 b2 = [100,-200,345]';

d = diag(A); %从矩阵A中提取对角元素 for i = 1:3 %写成对角阵 for j = 1:3 if i == j D(i,j) = d(i);

24

else D(i,j) = 0; end end end D

switch m case 1 t = b1; case 2 t = b2; end

%计算L矩阵 l = tril(A); L = -(l - D);

%计算U矩阵 u = triu(A); U = -(u - D);

B = inv(D - L)*U g = inv(D - L)*t

x{1} = orignal'; for k = 1:n

X{1} = B*x{1}+g;

%选择右端项 %计算迭代矩阵B %计算矩阵g 25

x{1} = X{1};

diedai = A*X{1} %使用迭代得到的值计算左端项

error = diedai - b1 %左端项与右端项的差距，作为迭代误差指标

if error == 0

step = k; %记录迭代步数 break; end end step

x1 = X{1}(1) %方程组的根x1 x2 = X{1}(2) %方程组的根x2 x3 = X{1}(3) %方程组的根x3 4.3 题3程序代码 a=1; b=3; h=0.5; y(1)=-2; x(1)=a; n=(b-a)/h+1; yy(1)=-2; for i=2:n

k1=(y(i-1)^2+y(i-1))/x(i-1);

k2=((y(i-1)+h*k1/2)^2+(y(i-1)+h*k1/2))/(x(i-1)+h/2); k3=((y(i-1)+h*k2/2)^2+(y(i-1)+h*k2/2))/(x(i-1)+h/2); k4=((y(i-1)+h*k3)^2+(y(i-1)+h*k3))/(x(i-1)+h);

y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6; % 四阶Runge-Kutta公式解 x(i)=x(i-1)+h; %有解区间的值

26

yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i)); %误差项 end [x' y' yy' s']

27

本文来源：https://www.bwwdw.com/article/62wx.html