MATLAB求解方程解析解和数值解

辽宁工程技术大学上机实验报告

用MATLAB求解质点振动方程

振动是日常生活和工程技术中常见的一种运动形式。利用常系数线性微分方程的理论来讨论有关自由振动和强迫振动的相关问题。利用MA TLAB数学软件大致可分四类情况：（1）无阻尼自由振动情况；（2）有阻尼自由振动；（3）无阻尼强迫振动；（4）有阻尼强迫振动求其数值解和解析解；

MATLAB软件求解微分方程解析解的命令“dsolve()”

求通解的命令格式：（’微分方程’,’自变量’）

注：微分方程在输入时，一阶导数y’应输入Dy，y’’应输入D2y等，D应大写。

1，无阻尼自由振动情况：常见的数学摆的无阻尼微小振动方程

代码如下：

>> t=0:pi/50:2*pi;

>> y=2*sin(3*t+2);

>>plot(t,y,'b')

2，有阻尼自由振动

由无阻尼振动的通解可以看出，无阻尼振动是按照正弦规律运动的，摆动似乎可以无限期的进行下去，但事实上，空气从在阻力，在运动时，我们必须把空气阻力考虑在内，所以我们得到有阻尼摆动方程为：

记u/m=2n,g/l=w^2,这里n,w是正常数，所以：

y=dsolve('D2y+2*n*Dy+w^2*y=0','t'); （4.43）

解得:

y = C3*exp(-t*(n + ((n + w)*(n - w))^(1/2))) + C2*exp(-t*(n - ((n + w)*(n - w))^(1/2)))

（1）小阻尼情形:n

y=exp(-n*t)*(c1*cos(w1*t)+c2*sin(w1*t))

和前面无阻尼的情形一样，可以把上式的通解改写为一下形式：

y=A*exp(-n*t)*sin(w1*t+Q), (4.45)

这里的A，Q为任意常数。

用matlab 操作得到：

t=0:0.1:10;

y=3*exp(-0.1*t).*sin(5*t+4);

plot(t,y,'k-')

如图：

由（4.45）可见，摆动的运动不是周期的，振动的幅度随着时间的增加而不断减小。 （2）大阻尼情形：你>w 时；r2

（3）临界的情形：即，n=w 的情形，方程（4.43）的通解解为 y=exp(-n*t)*(c1+c2*t), 这里的，c1，c2为任意常数； 由MAYLAB 绘制图像得： t=0:0.1:100;

yy=exp(-0.2*t).*(0.5-0.2*t); plot(t,yy,'r-')

01020

30405060708090100

3.有阻尼自由振动

无阻尼自由振动和有阻尼自由振动都属于自由振动，它对应于一个二阶常数齐次线性微分方程。当一个振动系统还经常受到一个外力作用时，这种振动称为强迫振动。

ψ＝Ａｓｉｎ（ωｔ＋θ）＋Ｈ／（ω２－ｐ２）ｓｉｎ（ｐｔ）

取Ａ＝２； ＝５；ｐ＝３； >> t=0:pi/50:2*pi;

>> y=2*sin(5*t+2)+1/9*sin(3*t); >> plot(t,y,'k')

1

2

3

4

5

6

7

-2.5

-2-1.5-1-0.500.511.52

2.5

4.有阻尼强迫振动

这时摆动的运动方程（1.11）变为：

D2y+2*n*Dy+w.^2*y=H*sin(p*t), （4.52）

根据实际情况，我们只讨论小阻尼的情况，即n

Y=A*exp(-n*t)*sin(w1*t+Q), (4.45) 这里的A ，Q 为任意常数,w1=sqrt(w^2-n^2) 现在 求（4.52）的一个特解，这时可以寻求如：

Y=M*cos(p*t)+N*sin(p*t)

的特解，这里M ，N 是待定系数，代入（4,53），（4,52）化简得：

Y=H*sin(Q)*cos(p*t)+H*cos(Q)*sin(P*t)=H*sin(p*t+Q),

因此，（4.52）的通解为：

Y=A*exp(-n*t)*sin(w1*t+Q)+H/sqrt((w^2-p^2)^2+4*n^2*p^2)*sin(p*t+Q); Matlab操作如下：

A=3;n=0.2;w1=2;Q=0.5;H=5;w=3;p=2;

>> t=0:0.1:50;

Y=A*exp(-n*t).*sin(w1*t+Q)+H/sqrt((w.^2-p^2)^2+4.*n.^2.*p^2).*sin(p.*t+Q); plot(t,Y,'k-')

05101520253035404550

利用MATLAB求解方程的数值解和解析解

用龙格-库塔法求方程数值解：

方程：i(t)=si-0.3i i(0)=0.02

s(t)=-si s(0)=0.98

建立函数M-文件：

function y=ill(t,x)

a=1;b=0.3;y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]';

建立程序命令：

ts=0:50;x0=[0.02,0.98];[t,x]=ode45('ill',ts,x0);[t,x]

plot(t,x(:,1),t,x(:,2));grid,pause