Foucault Pendulum

7 views (last 30 days)
Eric
Eric on 25 Oct 2011
Edited: arnold ing on 21 Jan 2018
I am inexperienced with matlab and have not used the ODE command before. I am trying write a code for the Foucault pendulum. I know there are example on the interent but they are not too helpful.
Here are my two differential equations and initial conditions
ddot(x) + (Omega^2)*x - 2*P*dot(y) = 0
ddot(y) - (Omega^2)*x - 2*P*dot(y) = 0
g = 9.81; % acceleration of gravity (m/sec^2)
l = 60; % pendulum length (m)
b = .2;
initial_x = .2; % initial x coordinate (m)
initial_y = b/2; % initial y coordinate (m)
initial_xdot = 0; % initial x velocity (m/sec)
initial_ydot = 0; % initial y velocity (m/sec)
Omega=2*pi/86400; % Earth's angular velocity of rotation about its axis (rad/sec)
lambda=0/180*pi; % (0, 30, 60, 90)latitude in (rad)
P=Omega*sin(lambda)
C=(g/l)^2;

Sign in to answer this question.

Accepted Answer

Grzegorz Knor
Grzegorz Knor on 25 Oct 2011
First of all you have to reduce the order of an ODEs:
http://www.mathworks.com/support/tech-notes/1500/1510.html#reduce
Next write fuction that describes the problem:
function dy = focaultPendulum(t,y)
% define your constants
dy = zeros(4,1);
dy(1) = y(2);
dy(2) = ... % your equation
dy(3) = y(4);
dy(4) = ... % your equation
And solve and plot:
[T,Y] = ode45(@focaultPendulum,tspan,initial_values);
plot(T,Y(:,1),T,Y(:,3))
  1 Comment
Eric
Eric on 25 Oct 2011
I am having an isssue reducing my 2nd order diff equs. Is there anything that stands out? Thank you for your help.
x_ddot = -C*x + 2*P*y'; %Equation of motion 1
y_ddot = -C*y - 2*P*x'; %Equation of motion 2
x_dot = z; % Equ 3; equals the first derivative of the free variable in Equ 1
z_dot = x_ddot; %Equ 4; taking the derivate of each side
y_dot = zz; % Equ 5; equals the first derivative of the free variable in Equ 2
zz_dot = y_ddot; % Equ 6; taking the derivate of each side
z_dot = -C*x + 2*P*zz; %Equ 7
zz_dot = -C*y - 2*P*z; %Equ 8
y(2)= x_dot;
y(4)= y_dot;
dy = zeros(4,1);
dy(1) = y(2);
dy(2) = z_dot;
dy(3) = y(4);
dy(4) = zz_dot;

Sign in to comment.

More Answers (2)

Grzegorz Knor
Grzegorz Knor on 26 Oct 2011
You have to rewrite your equations:
x'' + C*x - 2*P*y' = 0
y'' + C*y + 2*P*x' = 0
to form:
u = y' --> u' = y''
v = x' --> v' = x''
v' + C*x - 2Pu = 0
u' + C*y + 2Pv = 0
assume:
x(1) = x
x(2) = y
x(3) = u
x(4) = v
then:
function xprime = odetest(t,x)
% define P & C
xprime(1) = x(4);
xprime(2) = x(3);
xprime(3) = -2*P*x(4) - C^2*x(2);
xprime(4) = 2*P*x(3) - C^2*x(1);
xprime = xprime(:);
  2 Comments
Eric
Eric on 29 Oct 2011
I almost have the problem solved but my seoncd to last line keeps giving me an error saying "Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)to change the limit.
function dy = focaultPendulum(t,y)
% define your constants
g = 9.81; % acceleration of gravity (m/sec^2)
l = 60; % pendulum length (m)
b = .2; x = b; % initial x coordinate (m)
y = b/2; % initial y coordinate (m)
xdot = 0; % initial x velocity (m/sec)
ydot = 0; % initial y velocity (m/sec)
Omega = 2*pi/86400; % Earth's angular velocity of rotation about its axis (rad/sec)
lambda = 30/180*pi; % latitude in (rad) P = Omega*sin(lambda); C=(g/l)^2;
xddot = -C*x + 2*P*ydot; %Equation of motion 1
yddot = -C*y - 2*P*xdot; %Equation of motion 2
u = xdot;
udot = xddot;
v = ydot;
vdot = yddot;
udot = -C*x + 2*P*v;
vdot = -C*y - 2*P*u;
t0 = 0; %time initial tf = 8640;
%time final i.e. seconds in a day xy0 = [b b/2 0 0];
%initial conditions dy(1) = xdot;
dy(2) = udot;
dy(3) = ydot;
dy(4) = vdot;
dy(:);
[t,dy] = ode45(@focaultPendulum,[t0,tf], xy0); % plot(t,dy(:,1),t,dy(:,3))
bym
bym on 29 Oct 2011
change to dy = dy(:); or comment it out

Sign in to comment.

arnold ing
arnold ing on 21 Jan 2018
Edited: arnold ing on 21 Jan 2018
Can anyone help me with euler forward to solve foucault pendulum ODEs. i got an error calling the function E=euler_2(@edf1,@edff2,0,30,5/10,0,100)
function dy=edf2(z)
dy = zeros(4,1);
dy(2) = z(4);
dy(4) = -1.4580e-04*sin(pi/4)*z(2) - 0.9085*z(3);
end
function dy=edf1(z)
dy = zeros(4,1);
dy(1) = z(3);
dy(3) = 1.4580e-04*sin(pi/4)*z(4) - 0.9085*z(1);
end
function E=euler_2(f1,f2,x0,b,y01,y02,N)
% euler frw
% in [x0,b]
% step size h=(b-x0)/N
h=(b-x0)/N;
X=zeros(1,N+1);
Y1=zeros(1,N+1);
Y2=zeros(1,N+1);
X=(x0:h:b); % Rrjeti i pikave xi+1-xi=h
Y1(1)=y01; % Kushti fillestar y1(x0)=y01
Y2(1)=y02; % Kushti fillestar y2(x0)=y02
for i=1:N
Y1(i+1)=Y1(i)+h*feval(f1,X(i),Y1(i),Y2(i));
Y2(i+1)=Y2(i)+h*feval(f2,X(i),Y1(i),Y2(i));
end
E=[X' Y1' Y2'];
plot(X,Y1,'*',X,Y2,'o')
end

Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!