机电系统设计与仿真-系统数值仿真方法

（2）龙格-库塔（Runge-Kutta）法
基本思想是：在泰勒展开法中，泰勒展开式中 f(t,y) 的高 阶导数项的引入可提高数值积分精度。但 f(t,y) 的高阶导数难以 求解，可用函数值 f(t,y) 的线性组合来近似 f(t,y) 的高阶导数， 既可避免计算高阶导数，又可提高数值计算精度。
其中tk满足：a 出发点即为离散化，微分方程方程和时间区 间的离散化。
1、差商法
由初始值向后进行递推求解，求得数值序列。
2、泰勒展开法
当项数n=1时，即与差商法的结果相同。
3、数值积分法 （1）欧拉法 递推求解方法，简单、计算量小、需时少，属于单步法， 自启动法。具有一阶精度，精度较差。
机械振动系统数学模型 对于单自由度振动系统,动力学方程 ¨x+f(.x,x)=P(t) (1) 式中¨x、.x———位移x对时间的二阶和一阶导数; f(.x,x)、P(t)———单位质量物体上作用的恢复力和阻尼力的 合力与合外力。 取x1=x,x2=.x,则式(1)为 .x1=x2 (2) .x2=-f(x1,x2)+P(t) (3) 式(2)和式(3)组成状态方程。可调用相应的微分方程解题器 (Solver)进行时域仿真,并由数值结果以绘图命令plot绘出状态 变量随时间的变化曲线和相平面上的相轨迹。 对于多自由度线性振动系统,其动力学方程是以矩阵方程的形 式表达的。
（一）数值求解的基本概念
一阶常微分方程求解的初值问题
y(t ) f (t , y ) y(t0 ) y0 求该方程的解函数y(t)的数值解，即求函数y(t)在时间区间[a， b]上离散时间点tk（k=0，1，2，… ，N）处的近似值y0， y1，… ，yN，即 yk y(tk ), k 0,1, 2, , N
r为阶次，ki 称为龙格-库塔系数，bi为系数。
精度阶次r越大，计算精度越高，当然计算量会增大。 一阶龙格-库塔法时r=1，即为欧拉法递推公式；r>1时， 均有多个典型算法。如典型的二阶龙格-库塔递推公式有：
龙格-库塔法属于可自启动的单步法。
对于仿真中遇到的大多数工程实际问题，四阶龙格-库塔法 的精度已能满足要求，当步长h取得较小时，误差是很小的。
数值积分算法从原理上就不可避免地存在误差，在计 算机逐点计算时，初始数据的误差、计算过程中的舍入误差 等都会使误差不断积累，若这种误差积累能得到抑制，不会 随计算时间增加而无限增大，则可认为该计算方法是数值稳 定的，反之则是数值不稳定的。 最简单的判断方法（仅适于简单地估计）：取两种显 著不同的步距进行试算，若所得数据基本相同，则一般是稳 定的。 差分方程、离散状态方程稳定性的分析，与计算步距h （相当于采样周期T）关系密切。
机电系统设计与仿真
第3章 系统数值仿真方法
一、常微分方程的数值求解
控制系统数学模型主要是常微分方程形式，且实际遇到 的大部分微分方程难以得到解析解，通常采用计算机进行数值 计算求取数值解。 可采用高级仿真软件（如MATLAB）提供的数值求解函 数和程序段，但也需深入地了解常微分方程的数值求解问题， 从数值算法的底层考虑其编程实现过程，掌握数字仿真技术的 基本技能。 状态方程为一阶微分方程组形式，非常适于用计算机求 解其数值解，很容易直接对它编制相应的求解程序。传递函数 数学模型和高阶微分方程的直接求解很困难，需要转化为状态 方程模型（即实现问题）。
§数值算法中的病态问题 系统中的时间常数差别太大，造成仿真时间步长不好选 取。病态方程，病态系统。 需用多步法、隐式算法、自动变步长的数值算法来处理 病态方程，如隐式吉尔（Gear）法。
MATLAB的数值仿真方法 MATLAB仿真就是运用它所提供的强大的数值分析功能对 微分方程组的数值积分。它提供有多种积分方法,各有着不同的 功能和适用范围以供选择。 对于常微分方程(ODE)的解题器(Solver),可以在给定的初 始时间和条件下,通过数值方法计算每个程序步骤的解 (Solutions),并验证该解是否满足给定的容许误差(Tolerances), 如果满足,该解就是一个正确的解;否则就再试一次,直到求出解为 止。 下面是指令的格式及说明: 指令格式:[t,y]=solver(‘F’,tspan,y0) (1)solver为积分方法,如数值积分方法ode45等; (2)‘F’是描述常微分方程的ODE文件名,便于解题器调用。建立指 令格式如下: function ydot=F(t,y) ydot=[表达式1;表达式2;表达式3;……] 其中的“表达式n”以列向量形式组成为系统的状态方程(一阶微 分方程组);
§多步法
还有多步法，因要用到以前时刻的值，如前r步的数据， 不能从t=0自启动，需要用相同阶次精度的单步法来启动。获得 所需的前r步数据后，方可转入相应多步法。
因算式利用的信息量大，故比单步法更精确。 如梯形数值积分法、阿达姆斯法、吉尔法等为多步法，需 要其它数值积分方法进行启动。
§显式和隐式 显式比较直观，易于计算；隐式则需要数值迭代求解，先 (0) 用另一同阶次的显式公式估计出一个初值 yk 1 ，并求得 f k 1，再 (1) 用隐式公式求校正值 yk 1 ，若未达到所需精度要求，则再次迭 (i ) ( i 1) 代求解，只到相邻两次迭代值 yk 1 、yk 1 间的误差在要求的范 围内为止。故隐式法计算过程复杂，计算速度慢，但计算精度高， 对误差有较强的抑制作用，也利于数值稳定性。 §传递函数的离散化（如双线性变换法）与状态方程离散化 §数值稳定性 数值积分法求解微分方程，实质上为通过差分方程作为递 推公式进行。在连续微分方程离散化为差分方程中，原来的最小 相位稳定系统有可能稳定性变差，甚至变为不稳定，或非最小相 位系统。
(3)tspan为一个行向量,描述运算的起止时间:[t0,te]; (4)y0为初始状态值,用列向量来表示,维数为状态向量数。
MATLAB控制系统工具箱的时域响应函数有: initial(零输 入响应),impulse(单位冲激响应),step(单位阶跃响应),lsim(任意 输入响应)等;频域响应函数有:bode(系统频率响应———波德图 线),nyquist(系统奈奎斯特频率曲线)等。这些函数的参数既可为 状态方程形式,也可为传递函数形式。所以对于线性系统,除了以 上述解题器(Solver)进行仿真,还可根据振动系统与控制系统数 学模型的统一性,研究问题的相似性,直接借用MATLAB控制系统 工具箱函数对它进行仿真,为机械振动系统的仿真分析提供一个 功能强大、使用方便的工具。