disp(' ')
disp(' sdof_stiff2_ode45.m ')
disp(' ver 1.0  July 25, 2006 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the free vibration of a  ')
disp(' single-degree-of-freedom system with nonlinear stiffness, where ');
disp(' ')
disp('                   Stiffness Force = k x^3/abs(x)   ')
disp(' or equivalently   Stiffness Force = k x^2 sign(x)  ')
disp(' ')
disp(' The solution method is ode45, which is a numerical method.')
disp(' ')
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear t;
clear x;
clear w;
%
disp(' ')
v0 = input(' input initial velocity (meters/sec) ');
disp(' ')
d0 = input(' input initial displacement (meters) ');   
%
disp(' ')
m = input(' input mass kg ');
%
disp(' ')
k = input(' input stiffness coefficient (N/m^2) ');
%
omegan=sqrt(k/m);    %  units = (rad/sec)/sqrt(meters)
%
%   State Space Solution
%
fn=omegan/(2*pi);   % approximation only
period=1/fn;
%
dt=period/64;
sr=1/dt;
dur=40*period;
N=dur/dt;
%
% out6 = sprintf('\n dt = %f sec    dur = %f sec ',dt,dur);
% disp(out6)
%
t=linspace(0,dur,N);
t=t';
%
w(1)=d0;
w(2)=v0;
%
omegan2=omegan^2;
[t,w] = ode45(@sdof_stiff2_free_function,t,w,[],[],omegan2);
%
dmax=max(w(:,1));
vmax=max(w(:,2));
%
out6 = sprintf('\n max disp = %9.3g meters   max vel = %9.3g m/sec',dmax,vmax);
disp(out6)
%
period=6.87/(omegan*sqrt(abs(dmax)));
%
out6 = sprintf('\n period = %8.3g sec ',period);
disp(out6)
%
clear t;
clear w;
%
dur=period*8;
N=1000;
t=linspace(0,dur,N);
t=t';
w(1)=d0;
w(2)=v0;
[t,w] = ode45(@sdof_stiff2_free_function,t,w,[],[],omegan2);
%
disp(' ');
disp(' plot displacement?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(1);
    plot(t,w(:,1));
    xlabel('Time (sec)'); 
    ylabel('Disp (meters)');
    grid on;    
end
%
disp(' ');
disp(' plot velocity?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(2);
    plot(t,w(:,2));
    xlabel('Time (sec)'); 
    ylabel('Vel (meters/sec)');      
    grid on;   
end
%
disp(' ');
disp(' plot phase portrait?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(3);
    plot(w(:,1),w(:,2));
    ylabel('Vel (meters/sec)');
    xlabel('Disp (meters)');     
    grid on;   
end