工程数值分析报告实验龙格库塔最小二乘法.doc
实验一：最小二乘法拟合曲线实验一、实验名称：最小二乘法拟合曲线实验实验时间：2015-5-14实验地点： 主楼机房 实验器材： 计算机matlab 二、实验目的：学会用最小二乘法求拟合数据的多项式，并应用算法于实际问题。三、实验要求：1根据最小二乘法和加权最小二乘法的根本理论，编写程序构造拟合曲线的法方程，要求可以方便的调整拟合多项式的次数；2采用列主元法解1中构造的法方程，给出所拟合的多项式表达式；3编写程序计算所拟合多项式的均方误差，并作出离散函数 和拟合函数的图形；4 用MATLAB的内部函数polyfit求解上面最小二乘法曲线拟合多项式的系数与平方误差，并用MATLAB的内部函数plot作出其图形，并与1的结果进展比拟。四、算法描述实验原理与根底理论根本原理：从整体上考虑近似函数 同所给数据点 (i=0,1,m)误差 (i=0,1,m)的大小，常用的方法有以下三种：一是误差 (i=0,1,m)绝对值的最大值 ，即误差 向量 的X数；二是误差绝对值的和 ，即误差向量r的1X数；三是误差平方和 的算术平方根，即误差向量r的2X数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2X数的平方，因此在曲线拟合中常采用误差平方和 来 度量误差 (i=0，1，m)的整体大小。五、实验内容：共有两组给定数据，把给定的数据拟合成多项式。第一组给定数据点如表1所示如下：表1 数据表00506070809101175196219244271300表2 数据表005060708091011751962192442713001236421六、程序流程图开始输入拟合次数N,X,Y的坐标矩阵A计算X,Y平均值获得A的逆矩阵B获得矩阵C计算相关系数r解方程组获得系数a0,a1输出完毕七、实验结果>> zuixiaoerchenfaans =27-May-2015ans = 7.3611e+05ans = 1.0e+03 *>>Figure 1Figure 2Figure 4Figure 3八、实验结果分析实验程序quxiannihe.mclear alldate,now,clockx0=0.0 0.5 0.6 0.7 0.8 0.9 1.0;y0=1 1.75 1.96 2.19 2.44 2.71 3.00;w=ones(size(x0);x=0:0.01:1;%进展五次曲线拟合N=5;for i=1:N a1=LSF(x0,y0,w,i) ; y=polyval(a1,x); figure(i) plot(x0,y0,'ok',x,y,'r') title('最小二乘法'); legend('y0','y'); xlabel('x'); ylabel('y');end实验二：4阶经典龙格库塔法解常微分方程一、实验名称：4阶经典龙格库塔法解常微分方程实验时间：2015-5-14实验地点： 主楼机房 实验器材： 计算机matlab 二、实验目的：学习掌握4阶经典R-K方法，体会参数和步长对问题的影响。三、实验要求：1用4阶经典R-K法编写计算程序，要求用法与ode45一致。并将计算结果画图比拟，并分析步长变化对解的影响。2当激励力幅值分别按0.3，0.33，0.4， 0.43，0.54，0.58，0.75， 0.84，11.21，13.34进展计算。每一个数据画出三幅图，分别为时间位移曲线，时间速度曲线和相图。考察激励力幅值变化引起的系统响应的变化。3请采用MATLAB中的内部库函数ode45求解此常微分方程初值问题的解，并与1中的结果进展比拟。四、算法描述实验原理与根底理论系统方程和表述如下：如此系统的输出按如下求解：其中：这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。五、实验内容：求解常微分方程初值问题，考虑著名的Duffing方程。GDuffing 在1918 年引入了一个带有立方项的非线性振子来描述出现在许多力学问题中的质量、弹簧、阻尼系统。从那时起，Duffing方程在非线性动力学系统的研究中占有重要的地位。Duffing方程的标准形式是其中：是一个含有三次项的非线性函数，是一个周期函数。把，代入上式，可得 1式中：, 步长; 初值向量为:x0=(0, 0.1)。要求考察激励力幅值变化引起的系统响应的变化。积分时间区间为：0, 50。六、程序流程图开始输入a,b,n,x a,y y0输出k,x,yk=n=完毕七、实验结果i x(i) y(i)>>八、实验结果分析实验程序function varargout=saxplaxliu(varargin)clc,clearx0=0;xn=1.2;y0=1;h=0.1;y,x=lgkt4j(x0,xn,y0,h);n=length(x);fprintf('i x(i) y(i)n');for i=1:n fprintf('%2d %4.4f %4.4fn',i,x(i),y(i) );endfunction z=f(x,y)z=-2*x*y2;function y,x=lgkt4j(x0,xn,y0,h)x=x0:h:xn;n=length(x);y1=x;y1(1)=y0;for i=1:n-1 K1=f(x(i),y1(i); K2=f(x(i)+h/2,y1(i)+h/2*K1); K3=f(x(i)+h/2,y1(i)+h/2*K2); K4=f(x(i)+h,y1(i)+h*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);endy=y1;
