1、机械振动与噪声控制,江南大学机械工程学院 陈海卫 ,2010-12-28,机械振动与噪声学,第五章 振动方程的数值解法,2010-12-28,机械振动与噪声学,2010-12-28,机械振动与噪声学,第五章 振动方程的数值解法,5.1 MATLAB基础知识 5.2 微分方程的数值解法,2010-12-28,机械振动与噪声学,5.1 MATLAB基础知识,5.1.1 简介 5.1.2 MATLAB基本数学运算 5.1.3 MATLAB中矩阵的表示与运算 5.1.4 MATLAB语言编程 5.1.5 MATLAB图形绘制功能,2010-12-28,机械振动与噪声学,5.1.1 MATLAB简介,M
2、ATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、符号运算、绘制图形、创建用户界面，还可连接其他编程语言的程序等.MATLAB主要应用于工程计算、控制系统设计与仿真、信号分析与处理，图像处理等等。,2010-12-28,机械振动与噪声学,2010-12-28,机械振动与噪声学,5.1 MATLAB基础知识,5.1.1 简介 5.1.2 MATLAB基本数学运算 5.1.3 MATLAB
3、中矩阵的表示与运算 5.1.4 MATLAB语言编程 5.1.5 MATLAB图形绘制功能,2010-12-28,机械振动与噪声学,5.1.2 MATLAB基本数学运算,【例】求 的算术运算结果。(1)用键盘在MATLAB指令窗中输入以下内容 (12+2*(7-4)/32 (2)在上述表达式输入完成后，按【Enter】键MATLAB指令窗中将显示以下结果。ans =2 也可给运算式的结果设定一个变量x： x = (5*2+1.3-0.8)*102/25 x = 42,2010-12-28,机械振动与噪声学,MATLAB基本的算术运算,加 (+)、减 (-)、乘 (*)、除 (/)、幂次方 ()
4、， 例：5+3, 5-3, 5*3, 5/3, 53 例如： 1*2+3*4+5*6+7*8+9*10+11*12+.13*14+15*16ans =744 上例演示了多行指令的收入方法,2010-12-28,机械振动与噪声学,变量命名规则,1.变量名的大小写敏感。2.变量的第一个字符必须为英文字母，而且不能超过31个字符。3.变量名可以包含下连字符、数字，但不能为空格符、标点 。,2010-12-28,机械振动与噪声学,MATLAB常用数学函数,三角函数和双曲函数,2010-12-28,机械振动与噪声学,MATLAB常用数学函数,常用函数,复数函数,2010-12-28,机械振动与噪声学,其
5、他函数,2010-12-28,机械振动与噪声学,5.1 MATLAB基础知识,5.1.1 简介 5.1.2 MATLAB基本数学运算 5.1.3 MATLAB中矩阵的表示与运算 5.1.4 MATLAB语言编程 5.1.5 MATLAB图形绘制功能,2010-12-28,机械振动与噪声学,5.1.3 MATLAB中矩阵的表示与运算,直接输入A=1 2, 3; 4 5 6;7, 8 9A =1 2 34 5 67 8 9 冒号操作符a=0:1:10a =0 1 2 3 4 5 6 7 8 9 10,2010-12-28,机械振动与噪声学,特殊矩阵,eye(n) %生成n维单位方阵 eye(m,n
6、) %生成m行n列的单位阵 ones(n) %生成n维元素全为1的方阵 ones(m,n) %生成m行n列元素全为1的矩阵 zeros(n) %生成n维0方阵 zeros(m,n) %生成m行n列0方阵,2010-12-28,机械振动与噪声学,矩阵的运算,矩阵的代数运算 转置 B=A 加减乘 A+B A-B A*B 左除 AB 即AX=B的解X=A-1B 乘方 An=A*A*A,A 点运算 A.*B A.B A.n,2010-12-28,机械振动与噪声学,点运算的含义, A=ones(2)A =1 11 1 B=1 2;2 1B =1 22 1, C=A.*BC=1 22 1 C=A*BC =
7、3 33 3,注意比较两者运算结果的不同,2010-12-28,机械振动与噪声学,矩阵元素的访问, B =1 2;3 1B =1 23 1a=B(2,1) %第2行第1列的元素a=3,a=B(1,:) %第1行所有元素a =1 2a=B(:,1) %第1列所有元素a =1 3,2010-12-28,机械振动与噪声学,5.1 MATLAB基础知识,5.1.1 简介 5.1.2 MATLAB基本数学运算 5.1.3 MATLAB中矩阵的表示与运算 5.1.4 MATLAB语言编程 5.1.5 MATLAB图形绘制功能,2010-12-28,机械振动与噪声学,5.1.4 MATLAB语言编程,MAT
8、LAB中各种命令可以完成许多单一的任务，对于某些较为复杂的问题，仅靠现有的命令或函数来解决，往往是难以达到目的 。为此，要运用MATLAB编程语言编制程序，形成M-文件。程序是使计算机完成各项运算的命令集，运行一个编制好的程序，计算机会从第一条命令行开始，一行接一行地执行相应的命令，直到终止。程序编写调试完成后，可以存盘为文件。取名时要避开matlab专用命令名或专用变量名。,2010-12-28,机械振动与噪声学,5.1.4 MATLAB语言编程,MATLAB程序类型 for循环 while循环 If语句,2010-12-28,机械振动与噪声学,a.MATLAB程序类型,脚本文件,t=0:0
9、.001:1; x1=0.8*sin(10*t)+0.08*sin(24.5*t); x2=0.4*sin(10*t)-0.16*sin(24.5*t); plot(t,x1,t,x2),2010-12-28,机械振动与噪声学,a.MATLAB程序类型,函数文件,function y=mean(x) % MEAN average or mean Value % For vectors,Mean(x) is the mean value of X % For matrices, Mean(x) is a row vector containing % the mean value of each
10、 column m,n=size(x) if m=1m=n; end y=sum(x)/m;,2010-12-28,机械振动与噪声学,b.For循环,格式： for i=n1:(step):n2commands; end 作用：重复执行命令集commands.i为循环变量；n1为i的起始值；step为每次迭代i变化的量(即i=i+step)，其可正可负，可以是整数也可是小数，如省略后其默认值为1；程序循环的条件为i=n2.,2010-12-28,机械振动与噪声学,例：求奇数和：s=1+3+5+(2k-1)n=101;s=0;for i=1:2:ns=s+i;ends,例：求和：s=1+2+nn
11、=100;s=0;for i=1:ns=s+i;ends %当没有分号;时，将输出s的值。,2010-12-28,机械振动与噪声学,c.while命令,格式： while (condition is true)commands; end 作用：重复执行命令集commands.直到condition变为false.,2010-12-28,机械振动与噪声学,例：求和：s=1+2+ns=0;k=1;while k=3000s=s+k;k=k+1;enddisp(s) %用disp 函数输出计算结果,2010-12-28,机械振动与噪声学,d. If语句,单项选择控制 格式： if (conditio
12、n is true)commands; end 作用：若条件成立，则执行命令集 commands. 否则，不执行。,2010-12-28,机械振动与噪声学,例：求n个实数中最大的数M.n=length(a);M=a(1);for i=2:nif Ma(i)M=a(i);endenddisp(M),2010-12-28,机械振动与噪声学,d. If语句,多项选择控制 格式： if (condition is true)commands; elseif (condition is true)commands; else commands;end 作用：若条件成立，则执行命令集 commands.
13、否则，不执行。,2010-12-28,机械振动与噪声学,例：建立符号函数sign(x) x=10; if x0sn=1; elseif x=0sn=0; elsesn=-1; end disp(sn),2010-12-28,机械振动与噪声学,5.1 MATLAB基础知识,5.1.1 简介 5.1.2 MATLAB基本数学运算 5.1.3 MATLAB中矩阵的表示与运算 5.1.4 MATLAB语言编程 5.1.5 MATLAB图形绘制功能,2010-12-28,机械振动与噪声学,5.1.5 MATLAB图形绘制功能,常用的绘图函数,2010-12-28,机械振动与噪声学,x=0:0.001:1
14、0; % 0到10的1000个点的x座标y=sin(x); % 对应的y座标plot(x,y); % 绘图text(pi,0,leftarrow sin(pi),fontsize,18)hold on %图形保持，即在原基础上继续画图y=cos(x); plot(x,y,r:); % 绘图 r 表示红色, :表示虚线hold onx=0:0.2:10; y=sin(x); plot(x,y,ko); % 绘图 k 表示黑色 o表示圆圈legend(y=sinx,y=cosx,y=sinx) %建立图例xlabel(X)ylabel(Y),2010-12-28,机械振动与噪声学,2010-12-
15、28,机械振动与噪声学,线条和颜色控制符,2010-12-28,机械振动与噪声学,5.2 微分方程的数值解法,5.2.1 Euler法 5.2.2 四阶Runge-Kutta法 5.2.3 MATLAB中的ODE45函数 *5.2.4 固有频率与固有振型的求解,2010-12-28,机械振动与噪声学,Euler法,2010-12-28,机械振动与噪声学,2010-12-28,机械振动与噪声学,2010-12-28,机械振动与噪声学,h=0.001; y=pi/6 0; %赋初值 t=0; for i=0:2000dy=y(2) -9.8/0.2*sin(y(1);y=y+dy*h;t=t+h;
16、 %第i步时间的值 end,2010-12-28,机械振动与噪声学,clear all %清空系统以前所有的变量 close all %关闭所有的图形窗口 h=0.001; ry=zeros(2000,2); rt=zeros(2000,1); %分配空间 y=pi/6 0; %赋初值 t=0; for i=1:2000ry(i,:)=y; rt(i)=t;dy=y(2) -9.8/0.2*sin(y(1);y= y+dy*h;t=t+h; %第i步时间的值 end plot(rt,ry(:,1)*180/pi) %将摆角由弧度转为角度，画图 xlabel(time/s) %横坐标标题为时间
17、ylabel(theta/o) %纵坐标标题为摆角,2010-12-28,机械振动与噪声学,为简化程序的设计，可定义如下函数 function dy=danbai(t,y)dy为导数输出 ， t为时间输入，y为当前系统状态输入，其中y(1)为当前摆角theta的值y(2)为当前z的值，即当前theta导数的值dy=zeros(1,2); dy(1)=y(2); dy(2)=-9.8/0.2*sin(y(1); 此时，系统的源代码可以简化为,2010-12-28,机械振动与噪声学,clear all %清空系统以前所有的变量 close all %关闭所有的图形窗口 h=0.001; ry=ze
18、ros(2000,2); rt=zeros(2000,1); %分配空间 y=pi/6 0; %赋初值 t=0; for i=1:2000ry(i,:)=y; rt(i)=t;y= y+danbai(t,y)*h;t=t+h; %第i步时间的值 end plot(rt,ry(:,1)*180/pi) %将摆角由弧度转为角度，画图 xlabel(time/s) %横坐标标题为时间 ylabel(theta/o) %纵坐标标题为摆角,2010-12-28,机械振动与噪声学,采用Euler法的注意点,Euler法是一种较为简单的求解系统微分方程的方法，同样步长h的情况下，其求解的精度较低。 应用Eu
19、ler法时，迭代步长h不能太大，否则会造成求解过程的发散。,2010-12-28,机械振动与噪声学,5.2 微分方程的数值解法,5.2.1 Euler法 5.2.2 四阶Runge-Kutta法 5.2.3 MATLAB中的ODE45函数 5.2.4 固有频率与固有振型的求解,2010-12-28,机械振动与噪声学,5.2.2四阶Runge-Kutta法(RK4),2010-12-28,机械振动与噪声学,clear all %清空系统以前所有的变量 close all %关闭所有的图形窗口 h=0.01; ry=zeros(200,2); rt=zeros(200,1); %分配空间 y=pi
20、/6 0; t=0; %赋初值 for i=1:200rt(i)=t; ry(i,:)=y; k1= danbai(t,y); k2= danbai(t+0.5*h,y+0.5*k1*h);k3= danbai(t+0.5*h,y+0.5*k2*h); k4= danbai(t+h,y+k3*h); y= y+1/6*(k1+2*k2+2*k3+k4)*h;t=t+h; %第i步时间的值 end plot(rt,ry(:,1)*180/pi) %将摆角由弧度转为角度，画图 xlabel(time/s) %横坐标标题为时间 ylabel(theta/o) %纵坐标标题为摆角,2010-12-28
21、,机械振动与噪声学,四阶Runge-Kutta法的特点,四阶Runge-Kutta法是一种常用的用于求解常微分方程(组)的方法，同样步长h的情况下，其求解的精度和效率都远高于Euler法。以上介绍的Runge-Kutta法是一种定步长的算法，在保证精度的情况下，为提高计算效率，常常对其加以改进：利用高阶结果与当前结果间的误差动态的调整迭代步长h，以形成变步长Runge-Kutta法。常用的如Matlab当中的ODE45函数。,2010-12-28,机械振动与噪声学,5.2 微分方程的数值解法,5.2.1 Euler法 5.2.2 四阶Runge-Kutta法 5.2.3 MATLAB中的ODE
22、45函数 5.2.4 固有频率与固有振型的求解,2010-12-28,机械振动与噪声学,MALTAB中的ODE45函数,ODE45函数是Matlab自带的用于求解常微分方程(组)的变步长Runge-Kutta法。其使用方法如下：t,y=ode45(func,tspan,y0)式中，func为待求解系统的振动微分方程(注意是降阶之后)；tspan为迭代求解的时间范围,如 0 20表示的时间范围为0-20s；y0为系统初值。,2010-12-28,机械振动与噪声学,仍然以单摆系统为例 函数文件为 function dy=danbai(t,y) dy=zeros(2,1); dy(1)=y(2);
23、dy(2)=-9.8/0.2*sin(y(1);求解主文件为 clear all; close all; y0=zeros(2,1);y0(1)=pi/6;y0(2)=0; t,y=ode45(danbai,0,2,y0); plot(t,y(:,1)*180/pi),2010-12-28,机械振动与噪声学,例.二自由度系统数值解,2010-12-28,机械振动与噪声学,2010-12-28,机械振动与噪声学,程序代码如下 y0=zeros(4,1);x10=-0.1; x20=0.1; v10=0; v20=0; %系统初始位置与初始速度 y0(1)=x10; y0(2)= x20; y0(
24、3)=v10; y0(4)=v20 ; t,y=ode45(func1,0 10,y0)；函数文件func1代码如下 function dy=func1(t,y) x=y(1 :2) ; z=y(3 :4) ; M=7 -1 ;-1 7 ; K=12 -4 ;-4 12 ; dy=zeros(4,1) ; dy(1 :2)=z ; dy(3 :4)=inv(M)*K*x;,2010-12-28,机械振动与噪声学,5.2 微分方程的数值解法,5.2.1 Euler法 5.2.2 四阶Runge-Kutta法 5.2.3 MATLAB中的ODE45函数 5.2.4 固有频率与固有振型的求解,201
25、0-12-28,机械振动与噪声学,5.2.4固有频率与固有振型的求解,2010-12-28,机械振动与噪声学,2010-12-28,机械振动与噪声学,程序代码如下 M=7 -1 ;-1 7 ; K=12 -4 ;-4 12 ; D=inv(M)*K; v,e=eig(D); 运行后可得到 v =-0.7071 0.7071 归一化 v(:,1)= v(:,1)/ v(1,1)-0.7071 -0.7071 v(:,2)= v(:,2)/ v(1,2)e =1.3333 00 2.0000 式中e(1,1)与e(2,2)为系统第1阶与第2阶固有频率的平方。v(:,1)为系统的一阶固有振型，v(:,2)为系统的二阶固有振型。分别令v(:,1)/v(1,1), v(:,2)/v(1,2) 可得归1化之后系统的固有振型,2010-12-28,机械振动与噪声学,本章练习,熟悉Matlab编程环境与编程语言 练习本章例题,
