Problem Locating Event with ode45 - Programmers Heaven

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!


Welcome to the new platform of Programmer's Heaven! We apologize for the inconvenience caused, if you visited us from a broken link of the previous version. The main reason to move to a new platform is to provide more effective and collaborative experience to you all. Please feel free to experience the new platform and use its exciting features. Contact us for any issue that you need to get clarified. We are more than happy to help you.

Problem Locating Event with ode45

mosaranimosarani Posts: 1Member
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;

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.

title('Polar plot for motion in 2D rotating reference frame')

function dxdt = f(t,x);
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

Sign In or Register to comment.