When Lambda reaches a value of approximately 1.33, the code gives an error value.

When Lambda reaches a value of approximately 1.33, the code gives an error value.

Answer 1

[Torsten]

%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