- You have Rposn and Phiposn reversed in the polar plot.
- When you fix that, you'll see that Omega is way too large (the Omega inside f(t,x), not the redundant first Omega
- Your table is accelerating - do you want that? It's hard to keep the step size under control.
Event Location with ode45
3 views (last 30 days)
Show older comments
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
0 Comments
Answers (2)
Andrew Newell
on 4 May 2011
A few things:
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
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!