最小二乘法数值分析实验报告

1 / 38最小二乘法数值分析实验报告数学与信息工程学院 实 课程名称： 实 验 室： 实验台号： 班 级： 姓 名： 实验日期： 验 报 告 数值分析 2016 年 4 月 13 日 数值分析实验报告五 最小二乘法2 / 38 一、 题目 设有如下数据 用三次多项式拟合这组数据，并绘出图形。 二、 方法 最小二乘法 三、 程序 M文件： syms x f; xx=input(请输入插值节点 as x1,x2.n); ff=input(请输入插值节点处对应的函数值 as f1,f2.n); m=input(请输入要求的插值次数 m=);3 / 38 n=length(xx); for i=1:(m+1) syms fai x; H(i,j)=eval(fai); end end A=ff*(H)*inv(H*(H); syms x; f=0; for i=1:(m+1) f=f+A(i) hold on ezplot(f,xx(1),xx(n) 4 / 38四、 结果 save and run 之后： 请输入插值节点 as x1,x2. -3 -2 -1 0 1 2 3 请输入插值节点处对应的函数值 as f1 ,f2. - 请输入要求的插值次数 m=3 f = 133/100+121469856021/35184372088832*x-8042142191733/450359 9627370496*x +1020815915537309/9007199254740992*x 5 / 38五、拓展： 最小二乘法计算方法比较简单，是实际中常用的一种方法，但是必须经计算机来实现，如果要保证精度则需要对大量数据进行拟合，计算量很大。 数值分析实验报告 学院：专业：班级：学号：姓名： _ _ _ _ _ 2016年 1月 5日 实验一 MATLAB 在数值分析中的应用 插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。下面对插值中分段线性插值、拟合中的最为重要的最6 / 38小二乘法拟合加以介绍。 分段线性插值 所谓分段线性插值就是通过插值点用折线段连接起来逼近原曲线，这也是计算机绘制图形的基本原理。实现分段线性插值不需编制函数程序，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 线性7 / 38插值。 spline 三次样条插值。 cubic 三次插值。 所有的插值方法要求 x是单调的。x 也可能并非连续等距的。 正弦曲线的插值示例： x=0:10; y=sin(x); xi=0:10; yi=interp1(x,y,xi); plot(x,y,0,xi,yi) 则可以得到相应的插值曲线。 Matlab也能够完成二维插值的运算，相应的函数为interp2，使用方法与 interpl基本相同，只是输入和输出的参数为矩阵，对应于二维平面上的数据点，详细的用法见 Matlab联机帮助。 最小二乘法拟合 在科学实验的统计方法研究中,往往要从一组实验数据(xi,yi)中寻找出自变量 x 和因变量 y之间的函数关系y=f(x) 。由于观测数据往往不够准确，因此并不要求y=f(x)经过所有的点 (xi,yi),而只要求在给定点 xi上误8 / 38差?i?f(xi)?yi 按照某种标准达到最小，通常采用欧氏范数作为误差量度的标准。这就是所谓的最小二乘法。在MATLAB中实现最小二乘法拟合通常采用 polyfit函数进行。 函数 polyfit是指用一个多项式函数来对已知数据进行拟合，我们以下列数据为例介绍这个函数的用法： x=0:1； 2 y= - 为了使用 polyfit，首先必须指定我们希望以多少阶多项式对以上数据进行拟合，如果我们指定一阶多项式，结果为线性近似，通常称为线性回归。我们选择二阶多项式进行拟合。 P= polyfit (x, y, 2) P= 9 / 38- - 函数返回的是一个多项式系数的行向量，写成多项式形式为： ? 为了比较拟合结果，我们绘制两者的图形： 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现象 10 / 38Lagrange插值 实验目的 掌握 Lagrange插值公式的一般形式和构造方法，并利用Matlab绘出已知函数的 Lagrange插值图像。 实验步骤 已知点?x0,y0?,?x1,y1?xi,yi?xn,yn? Lagrange的基函数： ?x?x0?x?x1?x?xi?1?x?xi?1?x?xn?li?x? xi?x0xi?x1?xi?xi?1xi?xi?1?xi?xn Lagrange插值公式： Ln?x?li?x?yi i?0 11 / 38n 由以上公式可以在 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=k L=L*(z-x0(j)/(x0(k)-x0(j); end end s=s+L*y0(k); end y(i)=s; end y; 贵州师范大学数学与计算机科学学院学生实验报告 课程名称： 数值分析 班级： 09 数本一班 实验日期：学 号： 090704020163 82号 姓名： 指导教师： 12 / 38实验成绩：一、实验名称 实验二: Lagrange 插值与曲线拟合的最小二乘法 二、实验目的及要求 1. 掌握 Lagrange插值的基本思路和步骤. 2. 掌握最小二乘法的基本思路和拟合步骤. 3. 培养 Matlab编程与上机调试能力. 三、实验环境 每人一台计算机，要求安装 Windows XP操作系统，Microsoft officeXX、(或). 四、实验内容 1在区间-5,5内选取 n+1个的等距插值节点,取不同的n,1?x2 13 / 381作 n次 Lagrange插值,将对应的插值多项式图像与被插值函数 f(x)?画在同 21?x 一张图上进行比较. 2. 1. 对函数 f(x)? ,哪一种曲线拟合较好?为什么?能找出更好的拟合曲线吗? 五、算法描述及实验步骤 1. 分别画出原函数图像与 Lagrange多项式插值图像进行比较。 实验步骤： a.编写 Lagrange插值 M文件； 14 / 38b.画出原函数图像； c.调用 Lagrange插值取 n=10画拟合图像； d.比较观察两个图像。 2. 利用最小二乘法对给定数据点分别画一次，二次和三次多项式拟合曲线。 实验步骤： a.输入数据点； b.建立一个划分为四个部分的图像窗口； c.画一次多项式拟合图像在第一部分； d.画二次多项式拟合图像在第二部分； e.画三次多项式拟合图像在第三部分。15 / 38 六、调试过程及实验结果 七、总结 1、从图像可以看出用 lagrange插值函数拟合数据中间拟合的很好，但两边与原函数图象相比波动太大，逼近效果很差，出现所谓的 Runge现象。 2、从图像可以看出用最小二乘法去拟合较少的数据点，曲线拟合比直线拟合得好，高次的会比低次的拟合得好。 3.一般情形高次插值比低次插值精度高,但是插值次数太高也不一定能提高精度. 八、附录 1、 M文件： function cy=Lagrange(x,y,n,cx)16 / 38 m=length(cx);cy=zeros(1,m); for k=1:n+1 t=ones(1,m); for j=1:n+1 if j=k t=t.*(cx-x(j)./(x(k)-x(j); end end cy=cy+y(k).*t; end x=-5:5;17 / 38 y=1./(x. +1); plot(x,y) n=10; x0=-5:10/n:5; y0=1./(1+x0. ); cx=-5:5; cy=Lagrange(x0,y0,n,cx); hold on plot(cx,cy) 2、 x=-3,-1,0,1,3,5;18 / 38 y=-6,-3,-1,0,1,3; subplot(2,2,1) scatter(x,y,filled,r) hold on p1=polyfit(x,y,1); y1=polyval(p1,x); e1=norm(y-y1) e1 = t=-4:6; pt1=polyval(p1,t); 19 / 38 plot(t,pt1) title(一次多项式拟合图象) subplot(2,2,2) scatter(x,y,filled,r) hold on p2=polyfit(x,y,2); y2=polyval(p2,x); e2=norm(y-y2) e2 = t=-4:6; pt2=polyval(p2,t);20 / 38 plot(t,pt2) title(二次多项式拟合图象) subplot(2,2,3) scatter(x,y,filled,r) hold on p3=polyfit(x,y,3); y3=polyval(p3,x) y3 = - - - - y3=polyval(p3,x); e3=norm(y-y2)21 / 38 e3 = t=-4:6; pt3=polyval(p3,t); plot(t,pt3) title(三次多项式拟合图像) xxxx大学 数值分析实验报告 题 目：学 院：专 业：年 级：