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
