文章目录
计算机如何理解连续系统的动态特性？欧拉法求解微分方程龙格库塔法求解微分方程MATLAB代码编写和仿真效果计算机如何理解连续系统的动态特性？
一般连续系统的动态特性可以由一个微分方程，或者一组微分方程描述。因此，如果要对连续系统进行仿真，就需要对微分方程进行求解。求解微分方程，一般使用数值积分方法，那么计算机中积分怎么求呢？
设一个微分方程如下：
数值积分，就是求出定义域区间内若干个离散点处的近似解，然后相加，这就引出了欧拉法和龙哥库塔法。
欧拉法求解微分方程
欧拉法的思想，就是把积分曲线用折线加以近似。
如果能求出1点的导数，乘以一定的步长，再加上前一个点的值，就得到了后一个点的值。
如此往复迭代，就求出了每一个点的值。
该算法的优点：方法简单，计算量小。
该算法的缺点：精度低，适当减小步长可以改善这个问题。
龙格库塔法求解微分方程
和欧拉法不同，龙格库塔法使用泰勒级数进行展开
如果能将整个泰勒级数公式代入，精度无疑是最高的，但是为了计算方便，龙格和库塔两人提出了间接使用泰勒级数公式的方法，即用泰勒公式确定系数，然后乘上各阶的函数值。
实际工业使用上，四阶龙格库塔法的计算精度已经能够满足仿真要求，其递推公式如下：
简单地来说，两种方法的区别就是迭代时用的公式不同，欧拉法用的是一阶导数，而龙格库塔法用了4阶导数，计算更加精确。
MATLAB代码编写和仿真效果
求解如图所示SISO系统的输出
首先把传递函数转换成状态空间方程，然后求解。
% 龙格库塔法和欧拉法的比较clear all;close all;clc;% 欧拉法的计算步长h = 0.3;% 仿真步数L = 15/h;% SISO对象的零极点型z = [-1 -2];p = [-4 -0.5+j -0.5-j];k = 2.5;% 转换成状态空间[A,B,C,D] = zp2ss(z,p,k)% 输入和初值u = 1*ones(L,1);u0 = 0;% 对象的阶次n = length(p);% 龙格库塔和欧拉的状态初值xrk0 = zeros(n,1); % 3*1xer0 = zeros(n,1);for i = 1:Ltime(i) = i*h;% 欧拉法，更新单步xer = xer0 + h * (A * xer0 + B * u0);yer(i) = C * xer;% 更新数据，将当前值作为初值，便于下一次更新xer0 = xer;u0 = u(i);endfor i = 1:Ltime(i) = i*h;% runge_kutta龙格库塔法迭代公式k1 = A * xrk0 + B * u0;k2 = A * (xrk0 + h * k1/2) + B * u0;k3 = A * (xrk0 + h * k2/2) + B * u0;k4 = A * (xrk0 + h * k3) + B * u(i);xrk = xrk0 + h * (k1 + 2 * k2 + 2 * k3 + k4)/6;yrk(i) = C * xrk;% 更新数据xrk0 = xrk;u0 = u(i);endplot(time,yer,'b',time,yrk,'r');legend('euler','runge-kutta');
蓝色部分为欧拉法的计算输出，红色部分为龙格库塔法的计算输出。
适当调高精度，将h调整为0.1，可以看到欧拉法的精度有了明显的提高，不过在一定的时间之后，两者的值倾向于相同。
