数值分析上机实验报告

数值分析上机报告

学院：

专业：

学号：

姓名：

指导老师：

联系电话：

目录

1 序言 (3)

1.1 MATLAB语言特点 (3)

1.2 选用原因 (4)

2 第一题 (5)

2.1 题目 (5)

2.2 计算结果分析 (6)

3 第二题 (9)

3.1 题目 (9)

3.2 结果分析 (9)

4 第三题 (11)

4.1 题目 (11)

4.2 计算结果分析 (11)

附录一多项式拟合程序代码 (14)

附录二多项式插值程序代码 (15)

附录三jacobi迭代程序代码 (16)

附录四Gausssidel 迭代程序代码 (17)

附录五四阶Runge-Kutta算法程序代码 (18)

总结与体会 (19)

1序言

1.1 MATLAB语言特点

近年来MATLAB之所以能在各个行业得到广泛的应用，主要得益于MATLAB语言

能强大的工具箱；良好的扩充能力；完善的帮助系统。

1．高效的矩阵运算机制

句才能完成的操作，在MATLAB中直接使用矩阵即可完成计算，因而MATLAB在数据计算分析，特别是对海量数据的处理方面表现出相比其他的编程语言更大的优势。

2．多样化的操作途径

MATLAB语言为用户提供了多种操作方式选择。用户可以编写代码实现各种功能，代码可重复利用，同时，不擅长编程的用户也可以通过MATLAB图形界面操作，完成MATLAB的相应功能。

3．功能强大的工具箱

对于算法的开发，MATLAB提供了大量现成的函数，用户可以直接调用。MATLAB 软件对目前主流的算法都提供了现成的函数，并封装成一个个适用于不同领域的算法工具箱。常用的工具箱包括数学类、经济类、信号处理类工具箱，同时随着MATLAB软件版本的更新，不断有功能更强大的工具箱添加到MATLAB中。MATLAB 语言的这一特点，将大大节省算法开发的时间，用户无须详细了解算法的细节，只要掌握MATLAB中算法相应函数的调用即可。

4．良好的扩展能力

利用MATLAB语言编写的程序具有良好的扩展能力，可以方便地与各种编程语言链接。用户可以方便地在MATLAB中调用其他语言已编写好的程序，同时在其他语言中也可以方便地调用MATLAB的程序。MATLAB语言具有良好的接口编程技术。

5．完善的帮助系统

完善的帮助系统是MATLAB的又一突出特点，MATLAB向用户提供了多种帮助途径。通过MATLAB的帮助系统，用户可以获取MATLAB常用函数的使用方法及应用实例，而且这种帮助可以是实时的、在线的。同时，为了便于用户更好地使用MATLAB软件，在MATLAB中的主要算法都是可以直接看到源代码的。

1.2选用原因

本次数值分析上机实习采用MATLAB语言进行编程，主要是基于以下几方面的考虑：

1．编程效率高

MATLAB语言是一种面向科学与工程计算的高级语言，允许用数学形式的语言编写程序，且比BASIC、FORTRAN和C等语言更加接近我们书写计算公式的思维方式，用MATLAB编写程序犹如在演算纸上排列公式与求解问题。因此，也可通俗地称MATLAB语言为演算纸式科学算法语言。由于它编程简单，所以编程效率高，易学易懂。

2．高效方便的矩阵和数组运算

MATLAB语言像BASIC、FORTRAN和C语言一样规定了矩阵的算术运算符、关系运算符、逻辑运算符、条件运算符及赋值运算符，而且这些运算符大部分可以毫无改变的照搬到数组建的运算中。有些如算术运算符只要增加“.”就可以用于数组建间的运算。另外，它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库函数，使之在求解数字图像处理问题时显得大为简洁、高效、方便，这是其他高级语言所不能相比的。

3．用户使用方便

MATLAB语言是一种解释执行的语言(在没有被专门的工具编译之前)，它灵活、方便，起跳时程序手段丰富，调试速度快，需要学习时间少。人们用任何一种语言编程和调试一般都要经过四个步骤：编辑、编译、连接，以及执行和调试。各个步骤之间是顺序关系，编程的过程就是在他们之间作瀑布型的循环。MATLAB语言与其它语言相比，较好的解决了上述问题，把编辑、编译、连接和执行融为一体。它能在同一画面上进行灵活操作，快速排除输入程序的书写错误、语法错误甚至语义错误，从而加快了用户编写、修改和调试程序的速度，可以说在编程和调试过程中它是一种比VB 还要简单的语言，更不用谈C和C++。

4．扩充性强，交互性好

MATLAB语言有丰富的库函数，在进行复杂的数学运算时可以直接调用，且这些库函数同用户文件在形式上一样，所以用户文件也可以作为MATLAB 的库函数来调用。婴儿，用户可以根据自己的需要方便的建立和扩充新的库函数，提高MATLAB使用效率和扩充它的功能。另外，为了充分利用FORTRAN、C等语言的资源，包括用户自己编好的FORTRAN、C语言程序，通过建立M 文件的形式，混合编程，方便地调用有关的FORTRAN、C语言子程序，还可

以在C语言和FORTRAN语言中方便的使用MATLAB的数值计算功能。良好的交互性增加了代码的重用性。

综上所述，MATLAB最适用的范围是科学计算，对于本实习需要处理矩阵运算来实现图像处理算法来说，MATLAB再合适不过了。

2 第一题

2.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

2.2 计算结果分析

（1）根据实验数据，分别使用3,4,5,6次多项式拟合，可得拟合的多项式如下，并分别作图：

n=3；p=[-1.0326 19.334 -94.479 131.79]

n=4；p=[-0.38185 7.368 -42.143 73.533 0.74498]

n=5；p=[0.098075 -3.0789 34.502 -163.51 304.73 -139.5]

n=6;

p=[0.019359 -0.54079 5.1137 -16.897 -0.86696 66.375 -18.699]

上四幅图是多项式拟合曲线与原散点连线曲线的比较。明显可以看到，随着多项式拟合次数的增加，曲线的拟合度也越来越好，基本与连点曲线一致。

最好近似结果是n=6时多项式拟合，表达式为：

fx=0.019359*x6-0.54079*x5+5.1137*x4-16.897*x3-0.86696*x2+66.375*x1-18.69 9

（2）根据原散点构建拉格朗日插值多项式，对给出的xi的数值进行插值得到yi，在下图可以看到利用拉格朗日插值多项式插值的图像与曲线拟合度很好。

yi =

Columns 1 through 13

42.3840 41.4947 35.0742 24.3601 11.2792 -1.7683 -12.2977 -18.1626 -17.9118 -11.0210 2.0333 19.8565 40.3584

Columns 14 through 17

61.0794 79.5709 93.7788 102.3713

3 第二题

3.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,

3.2 结果分析

在以下的程序中，默认初始值为0，最大迭代次数N为50，精度为0.00001 （1）根据题意可得方程：

对此方程分别采用两种方法，通过运行jacobi.m与gausssidel.m程序进行收敛性分析：

Jacobi：满足误差限的迭代近似解为：-0.7273 0.8081 0.2525

需要迭代18次；

Gausssidel：满足误差限的迭代近似解为：-0.7273 0.8081 0.2525

需要迭代11次。

Jacobi：满足误差限的迭代近似解为：36.3636 -2.0707 114.0404

需要迭代24次；

Gausssidel：满足误差限的迭代近似解为：36.3636 -2.0707 114.0404

需要迭代16次。

通过以上计算分析，此方程的解收敛。根据迭代次数的比较，可以得到Gausssidel 迭代法要比Jacobi收敛性更快，更好。同时，上述两个方程系数阵都相同，但右端不同，用同样方法不同右端比较，得到结论，右端对迭代收敛有影响。

（2）

用jacobi到最大迭代次数50次时仍不能收敛。

用Gausssidel满足误差限的迭代近似解为：5.7692 0.7692 -4.2308需要迭代16次。

用jacobi到最大迭代次数50次时仍不能收敛。

用Gausssidel满足误差限的迭代近似解为：32.6923 7.6923 -42.3077需要迭代43次。

通过上述计算分析，用jacobi迭代时，两个方程都不能收敛，用Gaussssidel 迭代时，可以得到收敛的近似解，既，Gausssidel方法具有良好的收敛性。（3）

对此方程组，无论jacobi还是Gausssidel迭代法都不能得到收敛的近似解。可以推断是给出的行列式无法满足这两种方法的收敛条件。系数阵既不是行对角占优矩阵，也不是列占优阵，因此用这两种方法，得不到收敛的近似解。

4 第三题

4.1 题目

用Ru n ge-Kutt a4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。

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

4.2 计算结果分析

取不同步长得到求解的微分方程的解同精确解的误差分别如表1、表2、表3、表4、表5。

从上表中可以得出，在h=0.0001,0.001,0.01,0.1 时，绝对误差在逐步变小，因此在这种情况下认为此方法是稳定的。当h=0.2 时，绝对误差逐渐扩大，因此在h=0.2 时，算法不稳定。

在不同h下面的精度问题，取不同的步长求y(0.1)的值。从表6中可以看到步长越小，标准四阶Runge-Kutta算法得到的值精度越高，计算所需的步骤也越长。

%多项式拟合

xx= [1 2 3 4 5 6 7 8 9 10];

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

n=input('n=');

%判断x和y向量大小是否一致

if ~isequal(size(xx),size(yy))

error

end

x = xx(:);

y = yy(:);

%构建范德蒙德矩阵

V(:,n+1) = ones(length(x),1,class(x));

for j = n:-1:1

V(:,j) = x.*V(:,j+1)

end

%求解正规方程组的唯一解

%所得的解称为最小二乘多项式

%对矩阵V进行有选择的qr分解.

%当矩阵V为m×n并且m>n,那么只会产生具有前n列的正交矩阵Q

[Q,R] = qr(V,0);

p = R\(Q'*y); % p的保存 p = V\y;

r = y - V*p;

p = p.'; %p的转置一般采用行向量表示多项式系数

%S结构包括cholesky分解的范德蒙德矩阵R，自由度df，标准残差normr S.R = R;

S.df = length(y) - (n+1);

S.normr = norm(r);

%画图

X=xx;

Y=yy;

Y1=polyval(p,xx);

plot(xx,yy,':o',xx,Y1,'-*');

legend('yuanshi','nihe');

xx=[1 2 3 4 5 6 7 8 9 10];

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

xi=[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];

yi=[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]

%xi为标量或向量，被估计函数的自变量；

%yi为xi处的函数估计值。

x=xx(:);

y=yy(:);

xi=xi(:);

n=length(x);

m=length(y);

%插值点与它的函数值应有相同个数

if n~=m

error('The lengths of X and Y must be equal!');

return;

end

yi=zeros(size(xi));

for k=1:n

w=ones (size(xi));

for j=[1:k-1 k+1:n];

%输入的插值节点必须互异

if abs(x(k)-x(j))

error('the data is error');

return;

end

w=(xi-x(j))/(x(k)-x(j)).*w; %计算Lagrange基函数

end

yi=yi+w*y(k); %计算Lagrange插值函数

end

yi

%画图

X=xi;

Y1=load('f:\matlab7\shuzhifenxi\Y1.txt');

Y2=yi(:);

plot(xi,Y1,':o',xi,Y2,'-*');

legend('yuanshi','chazhi');

附录三jacobi迭代程序代码

a=input('请输入系数矩阵a：');

b=input('请输入矩阵b:');

N=input('请输入最大迭代次数N：');

esp=input('请输入近似解的误差限：');

if any(diag(a))==0

error('系数矩阵错误，迭代终止！')

end

D=diag(diag(a));

X0=zeros(size(b));

x1=0;

x2=0;

x3=0;

X1=[x1;x2;x3];

h=inv(D)*b;

B=inv(D)*(D-a);

B1=triu(B);

B2=tril(B);

k=1;

fprintf('雅可比迭代法 ');

t=0;

Y0=zeros(size(b));

while t<=N

Y1=h+B*Y0;

if norm(Y1-Y0,inf)

fprintf('满足误差限 \n')

break ;

end

Y0=Y1;

fprintf('第%2d次迭代得：',t)

disp(Y1');

t=t+1;

end

fprintf('满足误差限的雅可比迭代近似解为：')

disp(Y1');

附录四Gausssidel迭代程序代码

a=input('请输入系数矩阵a：');

b=input('请输入矩阵b:');

N=input('请输入最大迭代次数N：');

esp=input('请输入近似解的误差限：');

if any(diag(a))==0

error('系数矩阵错误，迭代终止！')

end

D=diag(diag(a));

X0=zeros(size(b));

x1=0;

x2=0;

x3=0;

X1=[x1;x2;x3];

h=inv(D)*b;

B=inv(D)*(D-a);

B1=triu(B);

B2=tril(B);

k=1;

fprintf('高斯-赛德尔迭代法 \n');

fprintf('第0次迭代得：')

disp(X1');

while k<=N

x1=h(1,1)+B1(1,:)*X0;

X1=[x1;x2;x3];

x2=h(2,1)+B1(2,:)*X0+B2(2,:)*X1;

X1=[x1;x2;x3];

x3=h(3,1)+B2(3,:)*X1;

X1=[x1;x2;x3];

if norm(X1-X0,inf)

fprintf('已满足误差限。 ')

break ;

end

X0=X1;

fprintf('第%2d次迭代得：',k)

disp(X1');

k=k+1;

end

fprintf('满足误差限的高斯-赛德尔迭代近似解为：') disp(X1');

附录五四阶Runge-Kutta算法程序代码

h=[0.0001,0.001,0.01,0.1]

k=1

for xk=0.1:0.1:1 %取不同的x进行计算

for i=1:1:4 %取不同的步长计算每个x

hj=h(i)

yk=1

for xkp=hj:hj:xk %四阶标准Runge-Kutta迭代算法

k1=hj*(-20*yk);

format long

k1;

k2=hj*(-20*(yk+k1*0.5));

k3=hj*(-20*(yk+k2*0.5));

k4=hj*(-20*(yk+k3));

ykm=yk+(k1+2*k2+2*k3+k4)*(0.16666666666667);

yk=ykm;

end

yk=ykm;

yeal=exp(-20*xkp) ; %求出精确值

ck=abs(yk-yeal); %精确值和计算值之间的绝对误差

ykeal(k)=yeal;

ckk(k)=ck;

yxk(k)=yk;

xkk(k)=xk;

k=k+1;

end

end

D=[yxk' ykeal' ckk']; %将迭代得到的结果及误差存到矩阵D中

for i=1:1:10 %将同一个步长下不同X得到的结果集中到一个矩阵

for j=1:1:3

Ah1(i,j)=D(4*i-3,j);

Ah2(i,j)=D(4*i-2,j);

Ah3(i,j)=D(4*i-1,j);

Ah4(i,j)=D(4*i,j);

end

Axk(i)=xkk(4*i-3);

end

Ch1=[Axk' Ah1]

Ch2=[Axk' Ah2]

Ch3=[Axk' Ah3]

Ch4=[Axk' Ah4]

Cheng=[Ch1(1,[2 3 4]);Ch2(1,[2 3 4]); Ch3(1,[2 3 4]);Ch4(1,[2 3 4])]

总结与体会

本次上机报告选择的是1-4题目中的第二题和第三题，5-6题中的第五题，在实习过程中，通过网上查阅资料，对程序代码进行了研究运行，因为对编程接触的机会太少，在编写过程中也是显得力不从心，也出现了不少问题，同时也很是纠结，最后通过向同学和师兄请教，解决了一些操作问题，通过本次实习，也对C语言及MATLAB进行了一次系统的复习，也有了更深一步的了解，但需要学习的还有很多，需要再接再厉，在克服了种种困难后，收获也是很多的。

总之，通过本次上机实习，学到了数值分析的相关知识，了解了编程运算的重要性，同时为自己专业研究也奠定了基础，编程是必不可少，通过理论与实践相结合的方法，不仅提高了计算与编程能力，还对将来的工作和进一步的学习打好基础。