一.对以下数据分别作二次,三次多项式拟合,并画出图形.
x=1:16;
y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5, 10.55,10.58,10.6];
答：程序如下
（1）x=(1:16);
y=erf(x);
p=polyfit(x,y,2);
f=polyval(p,x);
plot(x,y,x,f);
结果p=
-0.00100.02020.9096
（2）y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5, 10.55,10.58,10.6];
y=erf(x);
p=polyfit(x,y,3)
f=polyval(p,x);
plot(x,y,x,f)
结果
P=
0.0002-0.00710.06280.8404
二.在[0,4pi]画sin(x),cos(x)(在同一个图象中);其中cos(x)图象用红色小圆圈画.并在函数图上标注“y=sin(x)”,“y=cos(x)”,x轴,y轴,标题为“正弦余弦函数图象”.
答：程序如下
x=[0:720]*pi/180;
plot(x,sin(x),x,cos(x),'ro');
x=[2.5;7];
y=[0;0];
s=['y=sin(x)';'y=cos(x)'];
text(x,y,s);
xlabel('正弦余弦函数图象'),ylabel('正弦余弦函数图象')
图形如下
三.选择一个单自由度线性振动系统模型，自定质量、弹簧刚度、阻尼、激振力等一组参数，分别编程（m 文件）计算自由和强迫振动时的响应，并画出振动曲线图。

（要求画出该单自由度线性振动系统模型图）
其中质量为m=1000kg,弹性刚度k=48020N/m,阻尼c=1960N.s/m,激振力f(t)=0.阻尼比ζ的程序p=1960/(2*sqrt(48020*1000))
求得p=0.1414而p为阻尼比ζ
强迫振动时的响应程序
g =tf([-101],[48020048020*1.9848020]);bode(g)
图形
g =tf([001],[0001]);bode(g)
振动曲线图程序：
函数文件
function dx =rigid(t,x)dx =zeros(2,1);dx(1)=x(2);dx(2)=(-48020*x(1)-1960*x(2))/1000;
命令文件
options =odeset('RelTol',1e-4,'AbsTol',[1e-41e-4]);[T,X]=ode45(@rigid,[012],[11],options);plot(T,X(:,1),'-')
其图形如下024681012
-6-5
-4
-3
-2
-1
1
2
3
4
单自由度线性强迫振动系统模型图
其中质量为m=1000kg,弹性刚度k=48020N/m,阻尼c=1960N.s/m，f(t)=cos(3*pi*t)
振动曲线图程序：
函数文件
function dx=rigid(t,x)
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=(-48020*x(1)-1960*x(2))/1000+cos(3*pi*t);
命令文件
options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-4]);
[T,X]=ode45(@rigid,[020],[11],options);
plot(T,X(:,1),'-')
力等一组参数，建立Simulink仿真模型框图进行仿真分析。

答案：参数与第二题一样。

仿真参数设置：simulation/simulation parameters
Time：0，100s，solver options：variable-step，ode45
Simulink仿真模型框图如下：
五．自选一个力学问题（除单自由度线性振动系统外），应用MATLAB求解。

（要求有问题描述、数学模型、编制程序、结果显示）
解：问题描述
梁长L=2m,截面尺寸b*h为0.04m*0.02m,截面惯性矩为2.67*(10e-8)，密度为7920kg/，质量块的质量M为3.0kg弹性模量E=210GPa.求简支梁的模态.
采用集中质量法
力学模型如下:
数学模型：
质量矩阵为
柔度矩阵
矩阵D=FM称作矩阵的动力矩阵MATLAB程序：
m=[100;01.940;001]; h=0.00000195*3.186;
b=[9117;111611;7119]; f=h.*b;
a=f*m;
d=[111]';
for i=1:16;
y=a*d;
end
v=d
s=m*y;
结果：v=
1.0000
1.4271
1.0000
程序：
w1=1/sqrt(s(3,1))
subplot(3,1,1)
plot([02],[00])
hold on
plot([00.511.52],[0v(1,1)v(2,1)v(3,1)0]) p=v'*m*v;
q=v*v'*m;
j=a-(s(3,1)/p)*q;
d=[11-1]';
%x(1)=d;
for i=1:16;
y=j*d;
%l=y(i)
d=(1/y(3,1))*y;
end
v=d
s=m*y;
w2=1/sqrt(s(3,1))
subplot(3,1,2)
plot([02],[00])
hold on
plot([00.511.52],[0v(1,1)v(2,1)v(3,1)0]) q=v*v'*m;
r=j-(s(3,1)/p)*q;
d=[1-11]';
%x(1)=d;
for i=1:16;
y=r*d;
%l=y(i)
d=(1/y(3,1))*y;
end
v=d
s=m*y;
subplot(3,1,3)
plot([02],[00])
hold on
plot([00.511.52],[0v(1,1)v(2,1)v(3,1)0])
结果：
w1=
58.8629
v=
-1.0000
-0.0000
1.0000
w2=
283.6905
v=
1.0018
-0.7230
1.0000
w3=
524.8349。