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

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]程守洙，江之永.普通物理学下册（