Event Location with ode45

I'm trying to set up a simply program that models a billiard ball on a circular, rotating table. I would like to do something like the matlab package file ballode which models a bouncing ball. I'm having trouble locating when the integration reaches a point r=C (polar coordinates). I'd like to have the ball bounce off the edge of the circular table.

What am I doing wrong?

function billiard

%initial conditions
InitRSpeed = 5; 
InitPhiSpeed = 0; 
InitR = 3;
InitPhi = 4; 
y0 = [InitR; InitPhi; InitRSpeed; InitPhiSpeed];

%Rate of rotation of frame
Omega = 1;

%Time
t0 = 0;
tfinal = 20;
Deltat = 0.005*tfinal;
tspan = t0:Deltat:tfinal;

options = odeset('Events',@events,'RelTol',1e-5,'AbsTol',1e-4,'OutputSel',[1 3],'OutputFcn',@odeprint);

%Run DiffEQ
[t,x,te,xe,ie] = ode45(@f,tspan,y0,options);

% Extract the first and third columns of the solution array 
Rposn = x(:,1);
Phiposn = x(:,3);
% Plot the result. 

polar(Rposn,Phiposn)
title('Polar plot for motion in 2D rotating reference frame')

function dxdt = f(t,x);
Omega=1;
dxdt = [x(2);2*Omega*x(2)*x(4)+x(1)*Omega^2;x(1)*x(4);-2*Omega*x(2)];

function [value,isterminal,direction] = events(t,x)
% Locate when r=boundary and dr/dt positive 
value = (x(1)-200);     % Detect r = boundary
isterminal = 1;   % Stop the integration
direction = 1;   % positive direction only

Answer 1

[Andrew Newell]

A few things:

1. You have Rposn and Phiposn reversed in the polar plot.
2. When you fix that, you'll see that Omega is way too large (the Omega inside f(t,x), not the redundant first Omega
3. Your table is accelerating - do you want that? It's hard to keep the step size under control.

EDIT: Answer to section question - After the bounce, the radius (x(:,1)) decreases monotonically, going through zero and then past -15. To fix this, use these lines in events instead of yours:

value = (abs(x(1))-15); % Boundary is |R| = 15
isterminal = 1;   % Stop the integration
direction = 0;   % You're always approaching from inside

Answer 2

[Mosarani]

Thanks for the prompt response. I don't mean to take advantage of friendly advice but I've spent some time working on this and I don't quite understand

function billiard

%initial conditions
InitR = 3;
InitRSpeed = 2; 
InitPhi = 4; 
InitPhiSpeed = 0; 

x0 = [InitR; InitRSpeed; InitPhi; InitPhiSpeed];

%Time
t0 = 0;
tfinal = 22;
Deltat = 0.0005*tfinal;
tspan = t0:Deltat:tfinal;

%Run DiffEQ
options = odeset('Events',@events,'RelTol',1e-5,'AbsTol',1e-4,'OutputSel',[1 3],'OutputFcn',@odeprint);
[t,x,te,xe,ie] = ode45(@f,tspan,x0,options);

%This is my explicit re-running of the integration for test purposes
%It is completely continuous with previous integration but the event no  %longer triggers at r=15 and dr/dt>0
nt = length(t);
x0(1) = 15;
x0(2) = -1*x(nt,2);
x0(3) = x(nt,3);
x0(4) = x(nt,4);

[t,x,te,ye,ie] = ode45(@f,[tspan(nt) tspan(length(tspan))],x0,options);

% Extract the first and third columns of the solution array 
Rposn = x(:,1);
Phiposn = x(:,3);

% Plot the result. 
polar(Phiposn,Rposn)
title('Polar plot for motion in 2D rotating reference frame')

function dxdt = f(t,x);

Omega=.005;
dxdt = [x(2);2*Omega*x(4)+x(1)*Omega^2;x(1)*x(4);-2*Omega*x(2)];

function [value,isterminal,direction] = events(t,x)
% Locate when r=boundary and dr/dt positive 
value = (x(1)-15);     % Detect r = boundary
isterminal = 1;   % Stop the integration
direction = 1;   % positive direction only

When I call ode the second time the function no longer stops at r=15 with dr/dt positive as it should. Also, how can I bin the outputs of ode45 so that I can plot the entire function after a collision with the boundary?

If this could be explained to me I'm certain I could clean it up and implement loops for arbitrary collisions. I'm grateful for any guidance on this one.

[Andrew Newell]

See the edit to my answer.