详见mathworks

龙格库塔方法

写成矩阵（状态方程）的形式更简洁一点（其实这两种方法结果是一样的，如果C是[1,0,0]的话，就很明显了）

例如：求系统在0-5s内的单位阶跃相应，已知传递函数：

状态方程：

clear
%状态方程
A=[0,1
    -100,-10];
B=[0;
    1];
C=[100,0];
%初值
x=[0;0];y=0;t0=0;tf=5;  
h=0.05;%步长
t=t0;%从t0开始迭代，储存时间
N=round((tf-t0)/h);     %迭代步数
r=1;%单位阶跃输入
for i=1:N   
    k1=A * x+B*r;
    k2=A * (x+h*k1/2)+B*r;
    k3=A * (x+h*k2/2)+B*r;
    k4=A * (x+h*k3)+B*r;
    x=x+h*(k1+2*k2+2*k3+k4)/6;                    %采用四阶龙格库塔法
    y=[y,C*x];                                    %输出值
    t=[t,t(i)+h];
end
plot(t,y);

结果：

• matlab中传函转化状态空间方程的函数

构造传函：

M=100;%分子
N=[1,10,100];%分母
G=tf(M,N);%传函

转换：

    g=ss(G);%状态方程
    A=g.A;%各个系数矩阵
    B=g.B;
    C=g.C;

ODE45

• 函数句柄的创建

1、 直接@

A=@sin(x);
A(pi)

2、创建匿名函数

%定义函数f(x)=x^2
f = @(x) x^2;
f(2)

• ODE45使用
  [t,y] = ode45(odefun,tspan,y0)
  t：自变量
  y(t):因变量
  odefun：函数句柄(自动y’=)
  tspan:t的范围，如[0,10]
  y0：初值

1、一阶方程

a=@(t,y) 2*t;
[t,y] = ode45(a, [0 5], 0);
plot(t,y)

2、二阶方程

先转化为两个一阶微分：

创建m文件：

function dydt = vdp1(t,y)
dydt = [y(2); 
		(1-y(1)^2)*y(2)-y(1)];%括号不可省，不是变量，y是2x1列向量

调用ode45

[t,y] = ode45(@vdp1,[0 20],[2; 0]);
plot(t,y)

输出y两列，分别为y(1)和y(2）

一、计算公式

对于形如以下的常微分方程：

采用四阶龙格库塔法的计算公式：

四阶 龙格库塔法精度为4，属于单步递推法。

单步递推法的基本思想是从（x(i),y(i)）点出发，以某一斜率沿直线达到（x(i+1),y(i+1)）点。

二、实例计算

从上述定义可以看出，龙格库塔实质上是求一阶微分方程，但是如果将一阶导看作变量，则二阶导也不过是这个变量的一阶导而已。

接下来就看实例吧：

对于下述二阶方程：

1、基本思想：

令位移为   

q的一阶导，即位移的一阶导（速度）为  

q的二阶导  

求解位移u(1)的龙格库塔计算方法如下：

KK1=u(2);
KK2=u(2)+h/2*KK1;
KK3=u(2)+h/2*KK2;
KK4=u(2)+h*KK3;
u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);

求解速度u(2)的龙格库塔计算方法如下：

K1=-2*eptheton*u(2)-u(1)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
K2=-2*eptheton*(u(2)+h/2*K1)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
K3=-2*eptheton*(u(2)+h/2*K2)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
K4=-2*eptheton*(u(2)+h*K3)-(u(1)+h)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);

2、编程实现

参数设置：

h=0.001;       % 步长
t0=0:h:400;    % 总时长
w=5;
ep=0.02;
Fm=0.1;
Fah=0.05;
u(1)=0;u(2)=0;  % 初值

总的程序实现

h=0.001;
t0=0:h:400;
w=5;
ep=0.02;
Fm=0.1;
Fah=0.05;
u(1)=0;u(2)=0;
for i=1:length(t0)   % 进行多次迭代
    tao=t0(i);
    u=RK(u,tao,h,ep,w,Fm,Fah);     
    Result(i,:)=u;   % 将每次迭代的位移和速度保存
end
figure(1)
subplot(2,1,1)
plot(t0,Result(:,1))      % 绘制位移图
xlabel('Time')
ylabel('displacement')
subplot(2,1,2)
plot(t0,Result(:,2))      % 绘制速度图
xlabel('Time')
ylabel('velocity')

function u=RK(u,tao,h,ep,w,Fm,Fah)
KK1=u(2);
KK2=u(2)+h/2*KK1;
KK3=u(2)+h/2*KK2;
KK4=u(2)+h*KK3;
u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);
K1=-2*ep*u(2)-u(1)+Fm+Fah*cos(w*tao);
K2=-2*ep*(u(2)+h/2*K1)-u(1)-h/2+Fm+Fah*cos(w*tao);
K3=-2*ep*(u(2)+h/2*K2)-u(1)-h/2+Fm+Fah*cos(w*tao);
K4=-2*ep*(u(2)+h*K3)-u(1)-h+Fm+Fah*cos(w*tao);
u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);
end

结果图如下所示

值得特别注意的是

u=RK(u,tao,h,ep,w,Fm,Fah);

输入的u与输出的u一定要符号一致，从而确保下一次迭代的初始值是上一次的值。确保迭代能一直进行下去。

三、辅助验证

接下来用MATLAB自带的ode45函数来进行验证。

之前已经写过ode45函数的用法，在此不再进行介绍。

源程序如下：

t0=0:0.001:400;
w=5;
ep=0.02;
Fm=0.1;
Fah=0.05;
y0=[0 0];
[t,u]=ode45(@(t,u) odefun(t,u,w,ep,Fm,Fah),t0,y0);
figure(1)
subplot(2,1,1)
plot(t,u(:,1))
xlabel('Time')
ylabel('displacement')
subplot(2,1,2)
plot(t,u(:,2))
xlabel('Time')
ylabel('velocity')
function du=odefun(t,u,w,ep,Fm,Fah)
du=[u(2);
    -2*ep*u(2)-u(1)+Fm+Fah*cos(w*t)]; 
end

运算结果图如下所示

两中方法计算的结果是一样的。

创作不易，如有帮助，记得点个赞呐。