工程数值计算matlab实验报告

1、工程数值计算上机实验报告（ 第一 次）学生姓名 * 班级 * 学号 * 任课教师 * 上机时间 2019 年 10月 31日， 报告完成 2019 年 11月 2 日1、实验目的： 设计一个圆柱凸轮，底圆半径 R=200mm，凸轮的上端面不在同一平面上，圆柱的高度要根据动杆位移变化的需要进行设计和加工。工程上，可根据圆周等分9个点的高度值yi（见下表），确定圆周上任何一点的圆柱高度 y()，用于凸轮的加工和制造。角度04080120160200240280320360高度50.351.432.79.226.2214.723.632.541.450.3任务 1.1：分别用以下方法确定角度 35、

2、110、250对应的圆柱高度。 （1）待定系数法插值（提示：参考程序 1.1） （2）拉格朗日插值（参考程序：m2_4；m2_5）任务 1.2：分别用 3 次和 5 次多项式回归方法，拟合凸轮柱体高度轮廓函数 y()（提示：可采用常用函数 polyfit，参考程序 1.2）。2、计算方法：针对实验任务，结合课堂内容，说明解决方法，如：采用何种理论， 列出相关公式， 说明计算步骤， 写出程序框图等；(1)待定系数法插值计算原理：设多项式为Pnx=a0+a1x+a2x2+anxn(nm)将数据代如后得到线性方程组a0+a1x0+a2x02+anx0n=y0a0+a1x1+a2x12+anx1n=y

3、1a0+a1xm+a2xm2+anxmn=ym写成矩阵形式1x0x02x0n1x1x12x1n1x2x22x2n1xmxm2xmna0a1a2an=y0y1y2yn (写成V*a=yi)解得a=V-1*yi代入x=35、110、250；解得对应y值。程序框图如下：输入Vnn，i=1,j=1Vi,j=xi(i)(j-1)i=10?i=i+1j=10?j=j+1a=inv(V)*yi绘制图形计算x=35，110，250时的y输出否是否是(2)拉格朗日插值计算原理：给定n个插值节点和对应的函数值，利用n次拉格朗日插值多项式公式Lnx=k=0nyklkx,其中lkx=x-x0(x-xk-1)(x-xk

4、+1)(x-xn)(xk-x0)(xk-xk-1)(xk-xk+1)(xk-xn),可以得到插值区间内任意x的函数值y为y(x)=Lnx。从公式中可以看出，生成的多项式与用来插值的数据密切相关，数据变化则函数要重新计算，所以当插值数据特别多时，计算量会比较大。程序框图：输入S=0 , i=1 , j=1 i=10?建立矩阵t110j=10?I不等于j?t=txi-x(j)xi-x(j)j=j+1yi=s绘制（xi，yi）图形输出xi=35、110、250时yi值i=i+1S=s+t*y(i)是否否是是否(3)3次5次多项式回归计算原理：最小二乘法拟合程序流程图：输入数据使用polyfit函数拟

5、合多项式绘制图形3、程序设计：根据前面提到的计算方法编写程序；写出程序代码，并结合计算方法对程序中关键步骤进行必要的文字说明；(1)待定系数法插值clear;xi=0;40;80;120;160;200;240;280;320;360;%（角度数据）yi=50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3;%（高度数据）n=length(xi);V=ones(n);%（生成一个全是1的n阶方阵）for i=1:n %（将已知x0xn的0、1、2n次方填充到矩阵V中） for j=1:n V(i,j)=xi(i)(j-1); endenda=inv

6、(V)*yi; %（求矩阵V的逆矩阵与矩阵yi的值，即对应待定系数的值）xs=35;110;250;for i=1:length(xs) %（求35、110、250对应的函数值）ys(i)=a(1)+a(2)*xs(i)+a(3)*xs(i)2+a(4)*xs(i)3+a(5)*xs(i)4+a(6)*xs(i)5+a(7)*xs(i)6+a(8)*xs(i)7+a(9)*xs(i)8+a(10)*xs(i)9;endfigure;%（建立图形）plot(xs,ys,r*,xi,yi,ko-)%（绘制图形）ys%（输出35、110、250对应高度值）ys = 51.1889 13.3763 2

7、5.7382 %（35、110、250对应高度值输出结果）(2)拉格朗日插值x=0;40;80;120;160;200;240;280;320;360;y=50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.7;50.3;xi=1:0.01:360;xx=35;110;250s=0;for i=1:length(x) %（拉格朗日公式插值计算） t=ones(1,length(xi);%（建立一个列向量t110） for j=1:length(x) if j=i, t=t.*(xi-x(j)/(x(i)-x(j); end end s=s+t*y(i);end

8、yi=s;for i=1:length(xi)%（根据插值结果计算35、110、250对应高度值） if xi(i)=xx(1), yy(1)=yi(i) elseif xi(i)=xx(2), yy(2)=yi(i) elseif xi(i)=xx(3), yy(3)=yi(i) endendxx%（输出35、110、250）yy%（输出35、110、250对应高度值）plot(xi,yi,r-);%(绘制插值图形).(3.1)三次多项式回归x=0;40;80;120;160;200;240;280;320;360;y=50.3;51.4;32.7;9.22;6.22;14.7;23.6;3

9、2.5;41.7;50.3;a=polyfit(x,y,3);%（3阶多项式进行最小二乘法拟合）xi=1:0.01:360;yi=polyval(a,xi);%（多项式估值运算）plot(x,y,go,xi,yi,b-);%（绘制拟合曲线）(3.2)五次多项式回归x=0;40;80;120;160;200;240;280;320;360;y=50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.7;50.3;a=polyfit(x,y,5); %（5阶多项式进行最小二乘法拟合）xi=1:0.01:360;yi=polyval(a,xi); %（多项式估值运算）plot(x,y,go,xi,yi,b-); %（绘制拟合曲线）4、结果分析： 给出计算结果（ 可用数值、图表、曲线表示）， 并进行分析；（1）待定系数法插值1.1输出结果（插值图形）：（2）拉格朗日插值2.1输出结果（图形）：2.2输出结果（35、110、25