Foucault Pendulum

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;

Answer 1

[Grzegorz Knor]

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))

[Eric]

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;

Answer 2

[Grzegorz Knor]

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(:);

[Eric]

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]

change to dy = dy(:); or comment it out

Answer 3

[arnold ing]

 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