matlab仿真动力学方程的几种方法，总结，以范德波振子为例

1、 直接写方程 diff(y, 2) == k*(1 - y^2)*diff(y) - y，借助 odeToVectorField 和matlabFunction 和ode45函数求解。优点是，构造的方程非常直接；

syms y(t);
[V] = odeToVectorField(diff(y, 2) == k*(1 - y^2)*diff(y) - y);%把高阶变成一阶；
M = matlabFunction(V,'vars', {'t','Y'});
[t, y] = ode45(M,[0 30],[xxx yyy]);
n1=length(y);
xxx=y(n1,1);
yyy=y(n1,2);  %

2、把一个高阶微分方程，拆成多个一阶的微分方程组；并形成一个独立的文件。不好的地方是，怎么拆解，需要你有数学功底；

% 范德波振子的方程组；
% 有周期性外力 8*sin(4*t) ；
% y(1) 代表的就是y，振子的幅度；
% y(2) 代表的是y'，振子幅度的一阶导数；

function dydt = vdp1(t,y,params)

    k = params(1);  %外力的幅度和频率也可以是参数；
    dydt = [y(2); k*(1-y(1)^2)*y(2)-y(1)+8*sin(4*t)];

end

[T1,Y1]=ode45(@(t,X) vdp1(t,X,params),[0 5],[xxx;yyy]);

xxx，yyy是初值；

[0 5]是t的范围；

3、自己先把把一个高阶微分方程变成多个一阶的微分方程组，然后自己写  龙格库塔迭代方程。这种方法的优势，是每次迭代你都能参与，并且把想加的耦合项加入进去。前两种方法就做不到这一点。因为他默认耦合力都是已知的，实际上并不是。

%进行动力学迭代
for k=2:n
     %告诉迭代方程，当前神经元的编号
    [dx1, dx2] = dxdt_Vonderpol(k,x(k-1,1), x(k-1,2),1,adjMatix,x);

    x(k,1) = x(k-1,1) + dx1;
    x(k,2) = x(k-1,2) + dx2; 

end

%dx1/dt
function g=f1(t, x1, x2)    
    g = x2;
end

%dx2/dt
function g=f2(t, x1, x2, index, adjMatix, x)
    % t 是当前时间； 
    % x1 就是当前神经元的信号值；
    % x2 当前神经元信号的一阶导数；
    % index 代表的是当前神经元的编号；
    % adjMatix 代表的神经元之间的邻接矩阵；
    % x是过去的所有神经元状态；
    
    k = 3; %范德波振子的衰减系数    
    temp = k*(1-x1*x1)*x2 - x1;  %神经元的自发影响；
    
    % 根据神经元的index和adjMatix，计算其他神经元对当前神经元的影响；
    e = 0.9; %其他神经元对当前神经元的耦合项；   
    
    couple = 0;
    neuNum = size(adjMatix,1);
    
    for i = 1:neuNum
        con = adjMatix(index, i);% 获得两个神经元之间的连通性；
        if con == 1  %如果联通则进行耦合
            couple = couple + e * (x(t-1, 2*i-1) - x1);  %完成了神经元震荡的耦合；
        end
    end
    
    couple = couple/(neuNum-1);%其他神经元对当前神经元的耦合影响；    
    g = temp + couple;
    
end
 

范德波震荡方程的四阶龙格库塔函数

function [dx1, dx2] = dxdt_Vonderpol(t, x1, x2, index, adjMatix, x)

    % t是当前的时间；
    % x1是von der pol的 amp；
    % x2是von der pol的 amp的一阶 导数；
    % index是振子的编号；
    % adjMatix是振子之间的连接矩阵；

    h = 1e-2;   %步长

    K1 = f1(t, x1, x2);
    K2 = f1(t, x1 + h*K1/2, x2 + h*K1/2);
    K3 = f1(t, x1 + h*K2/2, x2 + h*K2/2);
    K4 = f1(t, x1 + h*K3, x2 + h*K3);

    L1 = f2(t, x1, x2, index, adjMatix, x);
    L2 = f2(t, x1 + h*L1/2, x2 + h*L1/2, index, adjMatix, x);
    L3 = f2(t, x1 + h*L2/2, x2 + h*L2/2, index, adjMatix, x);
    L4 = f2(t, x1 + h*L3, x2 + h*L3, index, adjMatix, x);

    dx1 = (K1 + 2*K2 + 3*K3 + K4)*h/6;
    dx2 = (L1 + 2*L2 + 3*L3 + L4)*h/6;
    
end