工程数值计算matlab实验报告

《工程数值计算》上机实验报告（第一次）
学生姓名 *** 班级 ******** 学号 *********
任课教师 *** 上机时间 2019 年 10月 31日，报告完成 2019 年 11月 2 日1、实验目的：
设计一个圆柱凸轮，底圆半径 R=200mm，凸轮的上端面不在同一平面上，圆柱的高度要根据动杆位移变化的需要进行设计和加工。

工程上，可根据圆周等分9个点的高度值y(θi)（见下表），确定圆周上任何一点的圆柱高度 y(θ)，用于凸轮的加工和制造。

任务 1.1：分别用以下方法确定角度 35°、110°、250°对应的圆柱高度。

（1）待定系数法插值（提示：参考程序 1.1）
（2）拉格朗日插值（参考程序：m2_4；m2_5）
任务 1.2：分别用 3 次和 5 次多项式回归方法，拟合凸轮柱体高度轮廓函数 y(θ)（提示：可采用常用函数 polyfit，参考程序 1.2）。

2、计算方法：针对实验任务，结合课堂内容，说明解决方法，如：采用何种理论，列出相关公式，说明计算步骤，写出程序框图等；
(1)待定系数法插值
计算原理：
设多项式为P n(x)=a0+a1x+a2x2+⋯a n x n(n≤m)
将数据代如后得到线性方程组
{a0+a1x0+a2x02+⋯+a n x0n=y0 a0+a1x1+a2x12+⋯+a n x1n=y1
⋮
a0+a1x m+a2x m2+⋯+a n x m n=y m
写成矩阵形式
[1x0x02⋯x0n
1x1x12⋯x1n
1x2x22⋯x2n
⋮⋮⋮⋱⋮
1x m x m2⋯x m n][
a0
a1
a2
⋮
a n]
=
[
y0
y1
y2
⋮
y n]
(写成V*a=yi)
解得a=V−1∗yi
代入x=35、110、250；解得对应y 值。

程序框图如下：
(2)拉格朗日插值
计算原理：
给定n 个插值节点和对应的函数值，利用n 次拉格朗日插值多项式公式L n (x )=∑y k l k (x )n k=0,其中l k (x )=
(x−x 0)⋯(x−x k−1)(x−x k+1)⋯(x−x n )(x k −x 0)⋯(x k −x k−1)(x k −x k+1)⋯(x k −x n )
,可以得到插值区间内任意x 的函
数值y 为y(x)=L n (x )。

从公式中可以看出，生成的多项式与用来插值的数据密切相关，数据变化则函数要重新计算，所以当插值数据特别多时，计算量会比较大。

程序框图：
输入
V nn ，i =1,j =1
V (i,j )=xi(i)(j−1)
i<=10?
i=i+1
j<=10?
j=j+1
a=inv(V)*yi
绘制图形
计算x=35，110，250 时的y
输出
否
是
否
是
(3)3次5次多项式回归
计算原理： 最小二乘法拟合
程序流程图：
输入
S=0 , i=1 , j=1
i<=10?
建立矩阵t 1×10
j<=10?
I 不等于j?
t =t ×
xi −x(j)
x (i )−x(j)
j=j+1
yi=s
绘制（xi ，yi ）图形 输出xi=35、110、250时yi 值
i=i+1
S=s+t*y(i)
是
否
否
是
是
否
输入数据
使用polyfit 函数拟合多项式
绘制图形
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 %（将已知x0⋯x n的0、1、2……n次方填充到矩阵V中）
for j=1:n
V(i,j)=xi(i)^(j-1);
end
end
a=inv(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;
end
figure;%（建立图形）
plot(xs,ys,'r*',xi,yi,'ko-')%（绘制图形）
ys%（输出35、110、250对应高度值）
ys =
51.1889 13.3763 25.7382 %（35、110、250对应高度值输出结果）
(2)拉格朗日插值
x=[0;40;80;120;160;200;240;280;320;360];
xi=1:0.01:360;
xx=[35;110;250]
s=0;
for i=1:length(x) %（拉格朗日公式插值计算）
t=ones(1,length(xi));%（建立一个列向量t1×10）
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
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)
end
end
xx%（输出35、110、250）
yy%（输出35、110、250对应高度值）
plot(xi,yi,'r-');%(绘制插值图形) ………………………………………………………………………………………………………………………………………………….
(3.1)三次多项式回归
x=[0;40;80;120;160;200;240;280;320;360];
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、250对应高度值）：
待定系数法插值与拉格朗日插值结果相近，从图中明显看出误差较大。

（3）3次多项式回归输出结果：
（4）5次多项式回归输出结果：
三次相比较于五次多项式回归，误差较小。

5、其他：说明实验中发现的问题及办法，个人体会、收获、思考等。

体会：熟练掌握编成，能通过计算机计算大大减少工作量。