弹簧阻尼双足行走模型步态周期极限环的Matlab

弹簧阻尼双足行走模型步态周期极限环的Matlab代码如下：

% 弹簧阻尼双足行走模型步态周期极限环
% 参考文献： https://ieeexplore.ieee.org/document/8953486

clear all;
close all;
clc;

% 参数设置
g = 9.8;
m = 75;
l = 0.9;
c = 1.2;
k = 8400;
b = 480;
h = 1.5;
alpha = 0.1;
beta = 0.3;
gamma = 0.2;
delta = 0.5;
epsilon = 0.8;
f = 1000;

% 状态变量初始化
x = [0.1; 0; 0.1; 0];

% 微分方程
dX = @(t,X) [X(2); ...
    (2*X(4)*k + 2*X(4)*c - k*x(1) - c*x(2) - g*m)/m; ...
    X(4); ...
    (-2*X(2)*k - 2*X(2)*c + k*x(1) + c*x(2) - b*(X(4)-X(2)))/m];

% 时间范围
tspan = [0 20];

% 画图
figure;
hold on;

% 循环步态
for j = 1 : 10

    % 初始条件
    x0 = x;
    t = 0;

    % 求解微分方程
    [T,X] = ode45(dX,tspan,x0);
    n = length(T);

    % 获取步态周期
    Tp = T(n);

    % 绘图
    plot(T,X(:,1),'b');
    plot(T,X(:,3),'r');
    plot([0 Tp],[0 0],'k--');
    plot([Tp Tp],[-h h],'k--');

    % 更新初始条件
    x = [alpha*X(n,1); beta*X(n,2); gamma*X(n,3); delta*X(n,4)];
end

xlabel('时间');
ylabel('状态变量');
title('弹簧阻尼双足行走模型步态周期极限环');
legend('x1','x2','Location','best');

注：代码参考了文献中的参数和初值，但可能与原文略有不同。