数值分析实验报告

《数值分析》实验报告学院：_________________ 专业：__________________班级：___________________学号：___________________姓名：___________________2015年1月5日实验一 MATLAB在数值分析中的应用插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。

随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。

下面对插值中分段线性插值、拟合中的最为重要的最小二乘法拟合加以介绍。

1.1分段线性插值所谓分段线性插值就是通过插值点用折线段连接起来逼近原曲线，这也是计算机绘制图形的基本原理。

实现分段线性插值不需编制函数程序，MATLAB自身提供了内部函数interp1其主要用法如下：interp1(x,y,xi) 一维插值◆yi=interp1(x,y,xi)对一组点(x,y) 进行插值，计算插值点xi的函数值。

x为节点向量值，y为对应的节点函数值。

如果y 为矩阵，则插值对y 的每一列进行，若y 的维数超出x 或xi 的维数，则返回NaN。

◆yi=interp1(y,xi)此格式默认x=1:n ，n为向量y的元素个数值，或等于矩阵y的size(y,1)。

◆yi=interp1(x,y,xi,’method’)method用来指定插值的算法。

默认为线性算法。

其值常用的可以是如下的字符串。

●nearest 线性最近项插值。

●linear 线性插值。

●spline 三次样条插值。

● cubic 三次插值。

所有的插值方法要求x 是单调的。

x 也可能并非连续等距的。

正弦曲线的插值示例： >> x=0:0.1:10; >> y=sin(x); >> xi=0:0.25:10; >> yi=interp1(x,y,xi); >> plot(x,y,’0’,xi,yi) 则可以得到相应的插值曲线。

Matlab 也能够完成二维插值的运算，相应的函数为interp2，使用方法与interpl 基本相同，只是输入和输出的参数为矩阵，对应于二维平面上的数据点，详细的用法见Matlab 联机帮助。

1.2 最小二乘法拟合在科学实验的统计方法研究中,往往要从一组实验数据),(i i y x 中寻找出自变量x 和因变量y 之间的函数关系y=f(x) 。

由于观测数据往往不够准确，因此并不要求y=f(x)经过所有的点 ),(i i y x ,而只要求在给定点i x 上误差i i i y x f -=)(δ按照某种标准达到最小，通常采用欧氏范数2δ作为误差量度的标准。

这就是所谓的最小二乘法。

在MATLAB 中实现最小二乘法拟合通常采用polyfit 函数进行。

函数polyfit 是指用一个多项式函数来对已知数据进行拟合，我们以下列数据为例介绍这个函数的用法：>> x=0:0.1:1；>> y=[ -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.489.30 11.2 ]为了使用polyfit，首先必须指定我们希望以多少阶多项式对以上数据进行拟合，如果我们指定一阶多项式，结果为线性近似，通常称为线性回归。

我们选择二阶多项式进行拟合。

>> P= polyfit (x, y, 2)P=-9.8108 20.1293 -0.0317函数返回的是一个多项式系数的行向量，写成多项式形式为：.92-+x-x8108031720.01293.为了比较拟合结果，我们绘制两者的图形：>> xi=linspace (0, 1, 100); %绘图的X-轴数据。

>> Z=polyval (p, xi); %得到多项式在数据点处的值。

当然，我们也可以选择更高幂次的多项式进行拟合，如10阶：>> p=polyfit (x, y, 10);>> xi=linspace (0, 1,100);>> z=ployval (p, xi);曲线在数据点附近更加接近数据点的测量值了，但从整体上来说，曲线波动比较大，并不一定适合实际使用的需要，所以在进行高阶曲线拟合时，“越高越好”的观点不一定对的。

实验二 Lagrange 插值与Runge 现象2.1 Lagrange 插值2.1.1实验目的掌握Lagrange 插值公式的一般形式和构造方法，并利用Matlab 绘出已知函数的Lagrange 插值图像。

2.1.2实验步骤已知点()()()()n n i i y x y x y x y x ,,,,,1100⋅⋅⋅⋅⋅⋅ Lagrange 的基函数：()()()()()()()()()()()n i i i i i i i n i i i x x x x x x x x x x x x x x x x x x x x x l -⋅⋅⋅--⋅⋅⋅---⋅⋅⋅--⋅⋅⋅--=+-+-11101110Lagrange 插值公式：()()∑=*=ni ii n y x l x L 0由以上公式可以在Matlab 中构造lagrange 函数，构造语法如下：function y=lagrange(x0,y0,x) n=length(x0); m=length(x); for i=1:m z=x(i);s=0; for k=1:n L=1;for j=1:n if j~=kL=L*(z-x0(j))/(x0(k)-x0(j)); end ends=s+L*y0(k); endy(i)=s; end y;由Lagrange 插值函数对x y sin = 做出的图形如下：2.2 Runge 现象2.2.1实验目的：利用Lagrange 插值公式对函数()211x x f +=进行插值，并绘出插值图像与原图像进行比较。

2.2.2实验步骤构建函数()x ffunction y = f(x)y = 1/(1+x*x);构建函数function y=lagrange(x0,y0,x)其构造方法即上面构造的lagrange 函数 执行主程序：a=[-10:0.2:10] ; for i=1:length(a) b(i)=f(a(i)) ;endc=[-5:1:5] ;for i=1:length(c)d(i)=f(c(i)) ;ende=[-5:0.2:5] ;z=lagrange(c,d,e);plot(a,b,'r',e,z)绘出图像如下：2.3 实验结论通过龙格现象可知，并不是插值多项式的次数越高（即插值节点越多）精度就越高在，这个实验中插值函数在两个端点处发生剧烈的波动，造成较大的误差，从数值计算上可解释为高次插值多项式的计算会带来舍入误差的增大，从而引起计算失真，因此，实际应用时，我们一般只用一次、二次最多是三次插值多项式，若想提高插值精度，我们可以采用分段插值。

实验三最小二乘拟合在科学工程试验中，经常需要从试验数据中寻找拟合曲线。

曲线拟合是指函数g(x)拟合给定的节点（xi,yj）,i=1,2,…,n,通常所拟合的节点数n必须大于未知数个数k。

确定函数g(x)参数，使得拟合函数与节点的偏差最小，这种方法称为最小二乘法。

当n=k时，由于拟合曲线通过所用节点，可使问题得到简化。

3.1实验目的明确曲线拟合的含义，用最小二乘法通过matlab进行曲线拟合。

3.2实验步骤编制以函数{}n k k x0=为基的多项式最小二乘拟合程序，并用于对下表中数据作3次多项式最小二乘拟合.x i-1.0 -0.5 0.0 0.5 1.0 1.5 2.0y i-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552取权数1≡iω，求拟合曲线∑=*= *nkkkx aϕ中的参数{}k a、平方误差2δ，并作离散数据{}i i y x,的拟合函数)(*xyϕ=的图形。

编写matlab代码进行三次多项式拟合，程序如下：x0=-1:0.5:2;y0=[-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552];n=3;alph=polyfit(x0,y0,n);y=polyval(alph,x0);r=(y0-y)*(y0-y)';x=-1:0.01:2;y=polyval(alph,x);plot(x,y,'k--');xlabel('x');ylabel('y0 * and polyfit. y-.');hold on;plot(x0,y0,'*');title('离散数据的多项式拟合');grid on;disp(['平方误差:',sprintf('%g',r)]); disp(['参数alph:',sprintf('%g\t',alph)])3.3结果分析输出结果：平方误差:2.17619e-005参数alph:1.99911 -2.99767 -3.96825e-005 0.549119根据给定的7个点的数据，所求的拟合函数的曲线可以基本地反映数据点的 其最小平方误差2δ为：2.17619e-005。

拟合的三次多次多项式为：549119.01096825.399767.299911.1)(5233+⨯--=-x x x x P 离散数据的拟合函数图形为：-1-0.50.511.52-5-4-3-2-1012345xy 0 * a n d p o l y f i t . y -.离散数据的多项式拟合实验四 复化求积公式计算定积分4.1实验目的利用复化求积公式计算下列定积分.（1）dx x ⎰+4=10211π； （2）dx xe e x ⎰=212； 4.2实验要求（1）若用复化梯形公式要求绝对误差限为710-⨯21=ε，利用它们的余项对算法做出步长的事前估计.（2）用复化梯形公式作出计算.（3）将计算结果与精确解做比较，并比较算法的计算量.4.3实验步骤4.3.1事前估计的Matlab 程序如下用复化梯形公式进行事前估计的Matlab 程序 （1）format long g x=0:0.01:1;f=8.*(3*x.^2-1)./(x.^2+1).^3;%二阶导函数x=1; %计算导函数最大值f=8.*(3*x.^2-1)./(x.^2+1).^3; h2=0.5*10^(-7)*12/f;h=sqrt(abs(h2)) %步长 n=1/hn=ceil(1/h)+1 %选取的点数（2）format long g x=1:0.01:2;f=2.*exp(x)+x.*exp(x); %二阶导函数x=2; %计算导函数最大值f=2.*exp(x)+x.*exp(x);h2=0.5*10^(-7)*12/f;h=sqrt(abs(h2)) %步长n=1/hn=ceil(1/h)+1 %选取的点数估计结果步长h及结点数n分别为（1）h = 0.000547722557505166n =1827（2）h = 0.000142479094906909n =70204.3.2积分计算的Matlab程序format long gpromps={'请选择积分公式，若用复化梯形，请输入T'};result=inputdlg(promps,'charpt 4',1,{'T'});Nb=char(result);result=inputdlg({'请输入积分式题号1-4：'},'实验4.1',1,{'1'});Nb_f=str2num(char(result));switch Nb_fcase 1fun=inline('4./(x.^2+1)');a=0;b=1;case 2fun=inline('x.*exp(x)');a=1;b=2;endpromps={'请输入用复化梯形公式应取的步长：'};result=inputdlg(promps,'实验4.2',1,{'0.01'});h=str2num(char(result));tic;N=floor((b-a)/h);detsum=0;for i=1:N-1xk=a+i*h;detsum=detsum+fun(xk);endt=h*(fun(a)+fun(b)+2*detsum)/2;time=toc;switch Nb_fcase 1disp('精确解：pi=3.14159265358979')disp(['绝对误差：',num2str(abs(t-pi))]);disp(['运行时间：',num2str(time)]);case 2disp('精确解：e^2=7.38905609893065')disp(['绝对误差：',num2str(abs(t-7.38905609893065))]);disp(['运行时间：',num2str(time)]);end4.4 结果分析当选用复化梯形公式时：（1）式运行结果为：t =3.14159261385336精确解：pi=3.14159265358979绝对误差：3.9736e-008运行时间：0.005（2）式运行结果为：t =7.38905611970610精确解：e^2=7.38905609893065绝对误差：2.0775e-008运行时间：0.007当选用复化梯形公式时，对步长的事前估计所要求的步长很小，选取的节点很多，误差绝对限要达到710*21-=ε时，对不同的函数n 的取值需达到1000-10000之间，计算量是很大。