%Imperfection
t0 = 0.01333;
%Half of a strut
length = 1.5;
%Define the length
a=0;
b=length;
s = linspace(a,b,100);
%Initial value of the load
lambda = 0.01;
%pcr=pi^2*EI/b^2;
%initial guess of the solution for the first step
%solinit = bvpinit(s,@init,lambda);
options = bvpset('RelTol',1e-6,'AbsTol',1e-6,'NMax',10000);
tDrive = 0.01*pi/8;
numberOfStepsOne = 2000;
%Loop in which the theta at the first point of the beam is increased
for i=1:numberOfStepsOne;
if i < 300
drive = t0+i*tDrive;
else
drive=t0+300*tDrive+(i*tDrive/100);
end
driveVector(i) = drive;
if i==1
solinit = bvpinit(s,@init,lambda);
else
solinit = bvpinit(s,@(x)deval(sol,x),lambda(i-1));
end
sol = bvp4c(@ode,@bc,solinit,options,drive,t0);
lambda(i) = sol.parameters;
x = 0:length/100:length;
y = deval(sol,x);
plotdata(i,:) = y(1,:);
end
%Plot the deformed shape of half the strut
figure(1)
x = 0:length/100:length;
plot(x,plotdata(1000,:));
%Plot the load-displacement curve
figure(2)
plot(driveVector,lambda)
function dtds = ode(s,t,lambda,drive,t0)
h=t0*(1-sin(s));
dtds = [t(2)
-lambda*(sin(t(1))+sin(h))];
end
function tinit = init(s)
tinit = [ 0.01333*(1-sin(s))
-0.01333*cos(s)];
end
function res = bc(ta,tb,lambda,drive,t0)
res = [ ta(1)-drive
ta(2)
tb(1) ];
end
