I am trying to calculate the response of a vibrating cantilever beam using ode45 function. There is no error in the code but when i run the code it takes to long to produce the result. I am not sure where i made a mistake. Can someone help me

MATLAB code
width = 5;
height = 0.1;
% The material is steel
E = 3.0e7;
nu = 0.3;
rho = 0.3/386;
% Calculate the coefficient matrix from the material properties
%G = E/(2.*(1+nu));
%mu = 2.0*G*nu/(1-nu);
I = height^3/12;
A = height*width;
% From beam theory, there is a simple expression for the lowest vibration
% frequency of the cantilever beam.
omegasq = 7.885^4*E*I/(rho*A*width^4);
beam=@(t,y)[y(2);(-(2*1*omegasq*y(2))-(omegasq^2*y(1)))];
  ic= [0;5];
  tspan = [0:5];
  [t,y]=ode45(beam,tspan,ic);
  plot(t,y), title ('Time vs Displacement')
  grid on
  hold on

Answer 1

[Steven Lord]

Your problem is stiff. Let's substitute in symbolic values for y and see what your ODE function will return.

>> vpa(beam(1, [sym('y1'); sym('y2')]))
ans =
                                                             y2
 - 1583163086609297.0*y1 - 79577963.950060874223709106445312*y2

A small change in either y1 or y2 results in a huge change in dy2/dt and thus a huge change in y2 at the next time step. So the ODE solver needs to take very, very small steps. You can see this if you plot your ODE:

[t,y]=ode45(beam,tspan,ic, odeset('OutputFcn', @odeplot));

Let that run for a minute or two and see that progress is being made, but it is being made very slowly. You may want to zoom in around y = 0, x between 0 and 0.01. Now let's try this same experiment but with a stiffer solver like ode15s:

[t,y]=ode15s(beam,tspan,ic, odeset('OutputFcn', @odeplot));

The plot (and the ODE solver) finished too quickly for me to see the individual time steps being plotted except as a bit of flicker.

[John D'Errico]

As this same question has been asked about a million times, I wonder if it might be a good idea to put a warning into ODE45. If the time step gets small, a warning message could pop up asking if the user is sure they should not be using a stiff solver instead? It seems like every time someone has a problem with ODE45 running slowly, the odds are very good it is because of a stiff system.

[Steven Lord]

I've entered this idea into the enhancement database.

[Jan]

@John: Or there is a pole, or a discontinuity.

[Aravind Krishnamoorthi]

@StevenLord Thank you for the reply. It worked out fine..

[Mia Finit]

@Steven Lord Hi, could you please help me with my problem:

ode45 to get time response (FFT) - MATLAB Answers - MATLAB Central (mathworks.com)