工程数值计算matlab实验报告

1/8《工程数值计算》上机实验报告（第一次）学生姓名***班级********学号*********任课教师***上机时间2019年10月31日，报告完成2019年11月2日1、实验目的：设计一个圆柱凸轮，底圆半径R=200mm，凸轮的上端面不在同一平面上，圆柱的高度要根据动杆位移变化的需要进行设计和加工。工程上，可根据圆周等分9个点的高度值()（见下表），确定圆周上任何一点的圆柱高度y(θ)，用于凸轮的加工和制造。角度04080120160200240280320360高度50.351.432.79.226.2214.723.632.541.450.3任务1.1：分别用以下方法确定角度35°、110°、250°对应的圆柱高度。（1）待定系数法插值（提示：参考程序1.1）（2）拉格朗日插值（参考程序：m2_4；m2_5）任务1.2：分别用3次和5次多项式回归方法，拟合凸轮柱体高度轮廓函数y(θ)（提示：可采用常用函数polyfit，参考程序1.2）。2、计算方法：针对实验任务，结合课堂内容，说明解决方法，如：采用何种理论，列出相关公式，说明计算步骤，写出程序框图等；(1)待定系数法插值计算原理：设多项式为()()将数据代如后得到线性方程组{写成矩阵形式[][][](写成V*a=yi)解得2/8代入x=35、110、250；解得对应y值。程序框图如下：(2)拉格朗日插值计算原理：给定n个插值节点和对应的函数值，利用n次拉格朗日插值多项式公式()∑(),其中()()()()()()()()(),可以得到插值区间内任意x的函数值y为()()。从公式中可以看出，生成的多项式与用来插值的数据密切相关，数据变化则函数要重新计算，所以当插值数据特别多时，计算量会比较大。程序框图：输入𝑉nn，i,j𝑉(𝑖,𝑗)𝑥𝑖(𝑖)(𝑗)i=10?i=i+1j=10?j=j+1a=inv(V)*yi绘制图形计算x=35，110，250时的y输出否是否是3/8(3)3次5次多项式回归计算原理：最小二乘法拟合程序流程图：输入S=0,i=1,j=1i=10?建立矩阵𝑡×j=10?I不等于j?tt×𝑥𝑖−𝑥(𝑗)𝑥(𝑖)−𝑥(𝑗)j=j+1yi=s绘制（xi，yi）图形输出xi=35、110、250时yi值i=i+1S=s+t*y(i)是否否是是否输入数据使用polyfit函数拟合多项式绘制图形4/83、程序设计：根据前面提到的计算方法编写程序；写出程序代码，并结合计算方法对程序中关键步骤进行必要的文字说明；(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阶方阵）fori=1:n%（将已知的0、1、2……n次方填充到矩阵V中）forj=1:nV(i,j)=xi(i)^(j-1);endenda=inv(V)*yi;%（求矩阵V的逆矩阵与矩阵yi的值，即对应待定系数的值）xs=[35;110;250];fori=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.188913.376325.7382%（35、110、250对应高度值输出结果）(2)拉格朗日插值x=[0;40;80;120;160;200;240;280;320;360];5/8y=[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;250]s=0;fori=1:length(x)%（拉格朗日公式插值计算）t=ones(1,length(xi));%（建立一个列向量×）forj=1:length(x)ifj~=i,t=t.*(xi-x(j))/(x(i)-x(j));endends=s+t*y(i);endyi=s;fori=1:length(xi)%（根据插值结果计算35、110、250对应高度值）ifxi(i)==xx(1),yy(1)=yi(i)elseifxi(i)==xx(2),yy(2)=yi(i)elseifxi(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];6/8y=[50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.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输出结果（插值图形）：7/8（2）拉格朗日插值2.1输出结果（图形）：2.2输出结果（35、110、250对应高度值）：待定系数法插值与拉格朗日插值结果相近，从图中明显看出误差较大。（3）3次多项式回归输出结果：8/8（4）5次多项式回归输出结果：三次相比较于五次多项式回归，误差较小。5、其他：说明实验中发现的问题及办法，个人体会、收获、思考等。体会：熟练掌握编成，能通过计算机计算大大减少工作量。