Personal tools

MainFile.m

MainFile.m — Objective-C source code, 2Kb

File contents

%%% This program illustrates the use of the Newton-Raphson method to find roots
% of a set of nonlinear equations.  The example presented here is from:
% "The Simplest Walking Model: Stability, Complexity, and Scaling" by Mariano
% Garcia, Anindya Chatterjee, Andy Ruina, and Michael Coleman published in the
% ASME Journal of Biomechanical Engineering, 1998.


global GAMMA TOE delta TSPAN


%%%%%%% INTEGRATION & LIMIT CYCLE SEARCH PARAMETERS
delta = 1e-6;                       % perturbation for gradient
TOE = 1e-12;                        % TOE = Tolerance of Error
                                    % end condition for fixed point search
TINITIAL = 0;                       % Time to start integration, (NDT)
TFINAL = 10;                        % Time to end integration, (NDT)
TSPAN = [TINITIAL TFINAL];          % Integration time span, (NDT)


%%%%%%% MODEL PARAMETERS
GAMMA = 0.004;                      % ground slope, (rad)

%%%%%%% STATE VECTOR (A guess that is hopefully close to a fixed point)
%%% NDT = Non Dimensional Time [seconds*sqrt(g/L)]
Xk    = zeros(4,1);
Xk(1) =  0.200;                     % theta, (rad)
Xk(2) = 2*Xk(1);                    % phi, (rad)
Xk(3) = -0.199;                     % thetadot, (rad/NDT)
Xk(4) = -0.0160;                    % phidot, (rad/NDT)


tic
%%%%%% EVALUATE STRIDEMAP
[t,F,Fendk] = StrideMap(Xk);       % Evaluates StrideMap
% t & F, are the time and state vectors at each integration step
% Fendk is the state of the system immediately after the impact

% If the walker's 1st step is successful, attempt to find a fixed point
if max(Fendk) ~= 0            
    %%%%%% NEWTON-RAPHSON ALGORITHM TO FIND FIXED POINT(S)
    [t,F,qstar] = NewtonRaphson(Xk,Fendk);
    
    toc
    
    % IF the NR algorithm finds a fixed point, check stability & record results
    if max(qstar) ~= 0
        Jqstar = StabilChk(qstar);                % check stability of FP
        [V,D] = eig(Jqstar);
        EigVals = diag(D);
        MaxAbsEigVals = max(abs(EigVals))         % stable if MaxAbsEigVals < 1

        figure(1)
        plot(t,F(:,1),'-r',t,F(:,2),'--b'),grid
        xlabel('\gamma, (rad)')
        ylabel('leg angles, (rad)')
        title('Leg Angles During the Gait Cycle')
        legend('Stance','Swing','location','northeast')

        figure(2)
        SF_elevation = cos(F(:,1))-cos(F(:,2)-F(:,1));
        plot(t,SF_elevation),grid
        axis([0 4 -0.025 0.025])
        xlabel('time, (none)')
        ylabel('elevation, (none)')
        title('Swing Foot Elevation During the Gait Cycle')
        
        figure(3)
        animate(F)
        
        % IF the NR algorithm, fails to converge, display failure message
    else
        disp('IC did not converge to a fixed point')
    end
    % OTHERWISE, if the walker's 1st step is NOT successful
else                          
    disp('Step not successful')
end

Document Actions