Foucault Pendulum

11 次查看（过去 30 天）
Eric
Eric 2011-10-25
编辑： arnold ing 2018-1-21
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;

请先登录，再回答此问题。

采纳的回答

Grzegorz Knor
Grzegorz Knor 2011-10-25
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 个评论
Eric
Eric 2011-10-25
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;

请先登录，再进行评论。

更多回答（2 个）

Grzegorz Knor
Grzegorz Knor 2011-10-26
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 个评论
Eric
Eric 2011-10-29
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 2011-10-29
change to dy = dy(:); or comment it out

请先登录，再进行评论。

arnold ing
arnold ing 2018-1-21
编辑：arnold ing 2018-1-21
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

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by