文档介绍：非线性单摆运动的数值解
() 文章编号 :1004 - 9762 200701 - 0094 - 03
Ξ
非线性单摆运动的数值解
陈向华 ,赵国忠
()包头师范学院 数学系 ,内蒙古 包头 014030
关键词 :阻尼 ; λλ ,.= = 2 12 2
2 非线性方程的数值计算结果 φπ图 1 阻尼较小的情形 , b < 2 k ,= / 2 0 ( ) ( () ) 为了讨论方程 3即方程 1的数值解, 不妨 φπFig. 1 The situation with smaller da mping ,= / 2 , 0 ( ) 假设 g = 9. 8 , l = 1 即 k = 9. 8, 单摆在初始时刻 b < 2 k φ d φ与铅垂线的夹角为 , 初始角速度= 0. 我们对0 d t
( ) 几种不同阻尼的情况分别求方程 3的数值解. 具体计算采用了MATLAB 下求解一阶微分方程组初值
( ) 问题数值解的最常用方法中的 ode45 函数 , 该函 数采用的是变步长的四阶五级 RUN GE KUTTA
[2 , 3 ] Felhberg 算法. 由经典的 RUN GE KUTTA 方法理
论分析可知 , 四阶的 RUN GE KUTTA 方法其收敛阶
5( ) ( 为 O h其中 h 为计算过程中所采用的时间步) 长, 下文计算选取的步长适当的小就可以保证所需 精度的解. 具体的数值化过程如下 :
首先编写一个 MATLAB 函数 dbhs. m 来描述系 φπ图 2 临界阻尼的情形 , b = 2 k ,= / 2 0 统的动态模型. 其内容为 φπFig. 2 The situation with smaller da mping ,= / 2 , 0 ()f unction f = dbhs t ,x b = 2 k () ()() ( dbhs = x 2, - b 3 x 2 - ksinx 1; 其中 b
) 和 k 要根据求解的需要输入不同的值
( ) 这时 , 可以调用微分方程数值解函数 ode45 对 dbhs. m 函数描述系统的动态模型进行数值求解. 其 大致过程为 :
t final = 15 ;
φ d ( ) φ; x= , 定义初值0 0 d t t = 0
( ) t ,x = ode45 ‘dbhs’, 0 ,t final ,x ; 0
() () plot t ,x. 数值结果可视化
( 不同阻尼的计算结果如图1 ,3 所示 在后文所
φ有的图形中 , b 代表阻尼 , f = 代表初始时刻单摆 1 0 π φ 图 3 阻尼较大的情形 , b > 2 k ,= / 2 0 与铅垂线之间的角度. 另外 , 为了比较 , 我们将对应 φπFig. 3 The situation with larger da mping ,= / 2 , 0 ) 不同阻尼的线性化解析解也画于相应的同一图中; b > 2 k
内 蒙 古 科 技 大 学 学 报 2007 年 3 月第 26 卷 第 1 期96
3 结果与讨论
( ) 众所周知 , 通常我们在无阻尼 b = 0的小角度
( φ φ) φ单摆运动中 即若 ?5,? 则