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

你可能感兴趣的