> MainFile.m

# 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(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
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