Pendulum without small angle approximation

I'm trying to build a numerical integrator to model a pendulum without using the small angle approximation. Our teacher provided the following code for the small angle approximation, and we've been asked to modify it:
%Script to plot one-dimensional motion (pendulum)
Dt=0.0001; %size of a time step.
Tmax = 10; %Length of time that should be computed.
Nmax=Tmax/Dt; %number of time steps in the entire time interval.
t(1)=0.0; %start at t=0.
theta(1)=0; %starting value for theta
v(1)=2.00; %6.2607; %starting value for v
a(1)=0; %fix the initial acceleration (needed to help loop below).
for n=1:Nmax %the integer n ranges between 1 and Nmax and identifies
%the time step.
theta(n+1)=theta(n)+v(n)*Dt; %advance the position by one time interval
a(n)=-9.8*sin(theta(n)); %use acceleration for real pendulum
v(n+1)=v(n)+a(n)*Dt; %advance the velocity by one time interval
t(n+1)=t(n)+Dt; %advance the time by on time
end
%commands to format the title and axis which work only after
xlabel('t','fontsize', 14, 'fontweight', 'b')
ylabel('\theta','fontsize', 14, 'fontweight', 'b')
title('Pendulum motion', 'fontsize', 14, 'fontweight', 'b')
hold on
plot(t,theta) %makes the plot
So far, I've tried replacing the theta(n+1) line to
theta(n+1)=theta(n)+(sqrt((2*g/L)*(cos(n)-cos(1))))*Dt
but there's an obvious problem because that difference of cosines will just stay at zero and theta will never change. I need to add the differential piece of theta, but how can I still account for the "kick" provided by the initial velocity v(1)=2?

1 Comment

I haven't taken mechanics classes in quite a few years, but
(1) The a(n) and v(n) in your teacher's code template don't look like angular accelerations and velocities, but rather linear accelerations and velocities. Hence, it's hard to see where
theta(n+1)=theta(n)+v(n)*Dt;
came from.
(2) There is no apparent small angle simplification in the original code. If there was supposed to be, wouldn't you have a(n)=-9.8*theta(n), with sin(theta) approximated by theta?
(3) Don't know where your proposed theta update formula came from, but you have cos(n) in there. You don't mean cos(theta(n)) ?

Sign in to comment.

Answers (0)

Asked:

Tim
on 18 Oct 2012

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!