MATLAB R―K方法对线性阻尼振动问题的求解和可视化

【摘要】利用MATLAB软件的龙格-库塔（Runge-Kutta）方法对线性阻尼振动的方程进行数值求解，再利用该软件绘制阻尼振动方程的数值解和解析解的图示、阻尼系数对阻尼振动的关系图示和欠阻尼、过阻尼、临界阻尼的图示。

【关键词】MATLAB;阻尼振动;欠阻尼;过阻尼;临界阻尼

龙格-库塔（Runge-Kutta）方法是目前较常用的一种微分方程求解的方法，特点是精度比较高。这种方法取多点的斜率进行加权平均作为平均斜率，再进行递推计算。

四阶R-K方法的近似计算公式为：

，

其中：

本文应用的是取四个斜率的四阶R-K方法[1]，其命令格式为：

其中，是写入了待解微分方程的命令文件的文件名;，和是自变量的初值和终值;为初始条件;为设定的误差，默认值为10-6;为对应的解向量。

1.阻尼振动方程的数值解和解析解

线性阻尼振动方程可写为：

设物体的质量m=10，弹簧的弹性系数k=10，初速度，初位置x0=10，角频率式，式中为阻尼系数。

微分方程的解析解为：

[2]

其中：

微分方程的数值解。令，x（0）=10，，v（0）=0。建立微分方程文件znzd.m，建立求解和绘图文件demo_znzd.m。

Function y=demo_znzd;t=0：0.1：15;w0=1;x0=1;v0=0;delta=0.5;w1=sqrt（w0^2-delta^2）;a=1;phi=0;x1=a*exp（-delta*t/6）.*cos（w1*t+phi）。

%以下是数值解：

t=0：0.1：15;x10=1;x20=0;x0=[x10，x20];[t，x]=ode45（'znzd'，t，x0）;plot（t，x1，t，x（：，1），'o'）;xlabel（'t'）;ylabel（'x'）;title（‘阻尼振动方程的数值解和解析解’）;legend（‘解析解’，’数值解’）。

以下是命令文件的内容：

function y=znzd（t，x）;y=[x（2）;-x（2）-x（1）]

执行求解文件，结果如图1所示。

2.阻尼系数与阻尼振动的关系

建立两个微分方程文件znzd.m 、znzd22.m和求解文件demo_znzd.m。

function y=demo_znzd22;

clc;t=0：0.05：10;x10=1;x20=0;x0=[x10，x20];

[t，x]=ode45（'znzd'，t，x0）;

[t，x22]=ode45（'znzd22'，t，x0）;

plot（t，x22（：，1），'k*'，t，x（：，1））;xlabel（'t'）;ylabel（'x'）;legend（‘阻尼系数较小’，’阻尼系数较大’）;

title（‘阻尼系数与阻尼振动的关系’）

以下是三个命令文件的内容：

function y=znzd（t，x）;y=[x（2）;-x（2）-x（1）]

function y=znzd22（t，x）;y=[x（2）;-0.02*x（2）-x（1）]

执行求解文件的结果如图2所示。

图1 阻尼振动方程的数值解和解析解

图2 阻尼系数

3.欠阻尼、过阻尼、临界阻尼

建立三个微分方程文件znzd.m、znzd0w.m、znzdw0.m，其中的，分别取、、znzd.m、znzd22.m和求解文件demo_znzdw0.m。

function demo_znzdw0;clc;t=0：0.5：30;x10=1;x20=0;x0=[x10，x20];

[t，x1]=ode45（'znzd'，t，x0）; [t，x2]=ode45（'znzd0w'，t，x0）;[t，x3]=ode45（'znzdw0'，t，x0）;

plot（t，x1（：，1），'r'，t，x2（：，1），'b：'，t，x3（：，1），'k*'）;xlabel（'t'）;ylabel（'x'）;

legend（‘欠阻尼’，’过阻尼’，’临界阻尼’）;title（‘欠阻尼、过阻尼、临界阻尼’）;

以下是三个命令文件的内容：

function y=znzd0w（t，x）;y=[x（2）;-30*x（2）-x（1）]

function y=znzd（t，x）;y=[x（2）;-x（2）-x（1）]

function y=znzdw0（t，x）;y=[x（2）;-2*x（2）-x（1）]

执行求解文件的结果如图3所示。

图3 过阻尼、欠阻尼、临界阻尼

4.结论

本文针对《大学物理》中阻尼振动这一教学难点进行分析，运用MATLAB对线性阻尼振动进行可视化，分别关注了阻尼振动方程的解析解和数值解的图示，阻尼系数与阻尼振动的关系图示，欠阻尼、过阻尼和临界阻尼的图示。

参考文献

[1]刘金远，段萍，鄂鹏.计算物理学[M].北京：科学出版社，2012，6.

[2]程守洙，江之永.普通物理学下册（第六版）[M].北京：高等教育出版社，2006，12.

[3]刘卫国，陈昭平，张颖.MATLAB程序设计与应用[M].北京：高等教育出版社，2002，6.

[4]刘延柱，陈立群.非线性振动[M].北京高等教育出版社，2001年8月.

基金项目：安徽省教育科学规划课题（项目编号：JG12203）;合肥师范学院质量工程项目（项目编号：2011jgsf02;2011jxsf05）。

作者简介：

王明美，女，江苏南京人，合肥师范学院电子信息工程学院副教授，主要从事普通物理、近代物理和计算物理的教学和研究。