数值分析上机报告

姓名：

学号：

专业：

联系电话：

本次数值分析上机实习采用数学软件。是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。在数值分析应用中可以直接调用软件中已有的函数，同时用户也可以将自己编写的实用导入到函数库中方便自己调用。基于数学软件的各种实用性功能与优点，本次数值分析实习决定采用其作为分析计算工具。

1．语言简洁，编程效率高

因为MATLAB定义了专门用于矩阵运算的运算符，使得矩阵运算就像列出算式执行标量运算一样简单，而且这些运算符本身就能执行向量和标量的多种运算。利用这些运算符可使一般高级语言中的循环结构变成一个简单的MATLAB语句，再结合MATLAB丰富的库函数可使变得相当简短，几条语句即可代替数十行C语言或Fortran语言语句的功能。

2. 交互性好，使用方便

在MATLAB的命令窗口中，输入一条命令，立即就能看到该命令的执行结果，体现了良好的交互性。交互方式减少了编程和调试的工作量，给使用者带来了极大的方便。因为不用像使用C语言和Fortran语言那样，首先编写源，然后对其进行编译、连接，待形成可执行文件后，方可运行得出结果。

3. 强大的绘图能力，便于数据可视化

MATLAB不仅能绘制多种不同坐标系中的二维曲线，还能绘制三维曲面，体现了强大的绘图能力。正是这种能力为数据的图形化表示(即数据可视化)提供了有力工具，使数据的展示更加形象生动，有利于揭示数据间的内在关系

在新版本中也加入了对C、FORTRAN、c++、JAVA的支持，使用时可以直接调用，也可将编写的实用程序导入到matlab函数库中方便以后使用时调用。

本次编程所用的软件为MATLAB，通过这次作业，对它有了初步的认识，以及对数值分析的体会更为深刻，希望为以后的学习和工作奠定一定的基。

不管是从书本中的理论还是实际的拟合分析中，我们知道在多项式拟合中，拟合的次数越高，则拟合函数越逼近原函数，精确度越高，此外，多项式拟合不会出现龙格现象，数值较为稳定，故该题中的最佳拟合为六次拟合。用低次插值多项式去近似被插值函数，即所谓的分段低次插值，其既具有一致收敛性也具有数值稳定性。插值法虽然保证了在节点处函数误差为零，但不一定能反映出被插值函数所对应曲线的总趋势，为此可采用函数逼近与曲线拟合，通过反应其曲线趋势的函数来近似它。

2 必做题二雅格比法迭代与高斯－赛德尔迭代

2.1题目

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

(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]；b1=[4,6]^T。

2.2分析过程

雅格比矩阵迭代格式：

针对

雅格比矩阵迭代格式：

迭代矩阵：

对应的谱半径：，故迭代收敛。

2.3计算结果

用Jacobi法与Gauss--Seidel法迭代法解两个方程组，分析其收敛性。分析结果见表1。为表明的准确性，随后将计算证明判断收敛与否的准确性。

表1.Jacobi 法与Gauss--Seidel法迭代结果对比表

1、对于采用Jacobi法时的问题二、三

问题二，：

根据题意有下列方程组：

计算迭代矩阵的特征值

，特征值；，所以，不收敛，与结论相同。

问题二，：

同理，的值与向量无关，所以，不收敛，与结论相同。

问题三，：

根据题意有下列方程组：

计算迭代矩阵的特征值

，特征值；，不收敛，与结论相同。

2、对于采用Gauss--Seidel法时的问题三

根据题意有下列方程组：

计算迭代矩阵的特征值

，特征值，；，不收敛，与

结论相同。

2.4 结果分析

1.右端项对收敛速度有明显影响；

2.Gauss--Seidel迭代相对Jacobi迭代法在这个问题上收敛速度更显著；

3.Jacobi迭代不收敛的方程用Gauss--Seidel迭代法也可能收敛

3 选做题一

3.1题目三次样条插值

给定函数，及节点，求其三次样条插值多项式（可取I型或II型边界条件），并画图及与的图形进行比较分析。

注：涉及到线性方程组求解问题需采用适当的求解算法。

3.2分析过程

采用三次样条插值多项式（取I型或II型边界条件），并与原函数曲线想比较。

3.3计算结果

从上图可以看出只选取10个样本点，三次样条插值拟合的结果与原函数吻合度很高，集体的数值见表2。

表2. 三次样条插值与原函数值比较表

3.4 结果分析

可以从三次样条插值图像中可以看得出来三次样条拟合的函数具有相当高的光滑度，与原函数吻合度较高。在遇到要求插值多项式拟合原函数的问题尤其对拟合多项式光滑度要求比较高时三次样条插值是理想的数值算法。

总结

通过此次数值分析上机实习，认识到了计算机编程在数值分析课程中的重要性。也让我意识到了学习编程语言的重要性。要想要将数学应用于实际工程中，特别是对于工科的学生来讲，这是一门非常重要的实践课程。

在此次报告中，首次接触了Matlab这门软件，之所以选这门，是因为很多工程中对数据的处理需要使用，这对我本身的专业是非常重要的。虽然之前学过C++，但是对于这门新的计算语言还是很头痛的，例如对函数的定义、调用以及很多形式方面跟之前的比较，还是有了很大的差异，这不得不使我在编程上面再次下功夫了。

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

附录

附录一：必做题一插值法代码

（1）三次多项式拟合：

x=[1,2,3,4,5,6,7,8,9,10]

y=[34.6588,40.3719,14.6448,-14.2721,-13.3570,24.8234,75.2795,103.5743,97.4847,78.2392]

subplot(1,1,1)

axis([0,10,-20,120])

P=polyfit(x,y,3)

Q1=sqrt(sum((y-polyval(p,x)).^2))

xi=0:0.1:10

yi=polyval(p,xi)

plot(x,y,'bo',xi,yi,'r')

（2）四次多项式拟合：

x=[1,2,3,4,5,6,7,8,9,10]

y=[34.6588,40.3719,14.6448,-14.2721,-13.3570,24.8234,75.2795,103.5743,97.4847,78.2392]

subplot(1,1,1)

axis([0,10,-20,120])

p=polyfit(x,y,4)

Q2=sqrt(sum((y-polyval(p,x)).^2))

xi=0:0.1:10

yi=polyval(p,xi)

plot(x,y,'bo',xi,yi,'r')

（3）五次多项式拟合：

x=[1,2,3,4,5,6,7,8,9,10]

y=[34.6588,40.3719,14.6448,-14.2721,-13.3570,24.8234,75.2795,103.5743,97.4847,78.2392]

subplot(1,1,1)

axis([0,10,-20,120])

p=polyfit(x,y,5)

Q3=sqrt(sum((y-polyval(p,x)).^2))

xi=0:0.1:10

yi=polyval(p,xi)

plot(x,y,'bo',xi,yi,'r')

（4）六次多项式拟合：

x=[1,2,3,4,5,6,7,8,9,10]

y=[34.6588,40.3719,14.6448,-14.2721,-13.3570,24.8234,75.2795,103.5743,97.4847,78.2392]

subplot(1,1,1)

axis([0,10,-20,120])

p=polyfit(x,y,6)

Q4=sqrt(sum((y-polyval(p,x)).^2))

xi=0:0.1:10

yi=polyval(p,xi)

plot(x,y,'bo',xi,yi,'r')

附录二：必做题二雅格比法迭代与高斯－赛德尔迭代代码

Jacobi法代码

clc

clear

A=[1,3;-7,1];

b=[4,6]';

[m,n]=size(A);

for j=1:m

for i=1:n

if j<i

L(j,i)=0;

elseif j==i

L(j,i)=0;

else

L(j,i)=-A(j,i);

end

for j=1:m

for i=1:n

if j>i

U(j,i)=0;

elseif j==i

U(j,i)=0;

else

U(j,i)=-A(j,i);

end

for j=1:m

for i=1:n

if j>i

D(j,i)=0;

elseif j<i

D(j,i)=0;

else

D(j,i)=A(j,i);

end

for j=1:m

for i=1:n

B(j,i)=(L(j,i)+U(j,i))./D(j,j);

end

for i=1:n

g(i,1)=b(i,1)./D(i,i);

end

x0=ones(m,1);

x=B*x0+g;

n=1;

while norm(x-x0)>=0.001

x0=x;

x=B*x0+g;

if n>200

disp('Warning:迭代次数太多，可能不收敛！');

break

end

n=n+1;

end

法

clc

clear

A=[1,3;-7,1];

b=[4,6]';

[m,n]=size(A);

for j=1:m

for i=1:n

if j<i

L(j,i)=0;

elseif j==i

L(j,i)=0;

else

L(j,i)=-A(j,i);

end

for j=1:m

for i=1:n

if j>i

U(j,i)=0;

elseif j==i

U(j,i)=0;

else

U(j,i)=-A(j,i);

end

for j=1:m

for i=1:n

if j>i

D(j,i)=0;

elseif j<i

D(j,i)=0;

else

D(j,i)=A(j,i);

end

B=(D-L)\U;

g=(D-L)\b;

x0=ones(m,1);

x=B*x0+g;

n=1;

while norm(x-x0)>=0.001

x0=x;

x=B*x0+g;

if n>200

disp('Warning:迭代次数太多，可能不收敛！');

break

end

n=n+1;

End

附录三：选做题一三次样条插值代码

三次样条差值子

function[yy b c d]=spline3(x,y,xx,flag,vl,vr)

if length(x)~=length(y)

error('输入数据应成对!');

end

n=length(x);

a=zeros(n-1,1);

b=a;d=a;dx=a;dy=a;

A=zeros(n);B=zeros(n,1);

for i=1:n-1

a(i)=y(i);

dx(i)=x(i+1)-x(i);

dy(i)=y(i+1)-y(i);

end

for i=2:n-1

A(i,i-1)=dx(i-1);

A(i,i)=2*(dx(i-1)+dx(i));A(i,i+1)=dx(i);

B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));

end

if flag==0

A(1,1)=1;

A(n,n)=1;

end

if flag==1

A(1,1)=2*dx(1);A(1,2)=dx(1);

A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1);B(1,1)=3*(dy(1)/dx(1)-vl);

B(n,1)=3*(vr-dy(n-1)/dx(n-1));

end

if flag==2

A(1,1)=2;A(n,n)=2;

B(1,1)=vl;B(n,1)=vr;

end

c=A\B;

for i=1:n-1

d(i)=(c(i+1)-c(i))/(3*dx(i));

b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;

end

[mm nn]=size(xx);

yy=zeros(mm,nn);

for i=1:mm*nn

for ii=1:n-1

if xx(i)>=x(ii)&xx(i)<x(ii+1)

j=ii;

break;

elseif xx(i)==x(n)

j=n-1;

end

yy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3;

end

三次样条差值主

clear;clc;

format short g;

x1=[-5 -4 -3 -2 -1 0 1 2 3 4 5]';

y1=[0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.0385]';

u1=0.03846;

un=0.03846;

xx1=[x1(1):0.001:x1(end)]';

[yy1 b1 c1 d1]=spline3(x1,y1,xx1,1,u1,un);

plot(x1,y1,'bo',xx1,yy1,'--k');

grid on;

hold on;

x=-5:0.1:5;

y=1./(1+x.^2);

plot(x,y,'r');

legend('样本点','三次样条插值','y=1/(1+x^2)');

相关推荐