Event Detection on a circle

1 view (last 30 days)
Danny Brett
Danny Brett on 10 Dec 2021
Answered: KSSV on 11 Dec 2021
I have an ODE trevelling over the surface area.
The ODE starts on the circle however I want the plot to end once it hits the circle again.
Below is the code I am currently using and the photo shows the plot I am currently getting.
The event detection part of the code is the function right at the bottom but I can't seem to get it working.
Any help would be great thanks!
clc
tspan = [0:0.05:80]; %Range for independent var. i.e. time
r = 2;
centreX0 = 0;
centreY0 = 0;
theta = 1.08;
Ptheta = 3;
Pr = -sqrt(2 - 2 * r.^3 .* cos(theta).^3 + 6 .* r.^3 .* cos(theta) .* sin(theta).^2 - (Ptheta.^2 ./ r.^2));
%options = odeset('Events',@MonkeyPlot);
% y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
a = r .* cos(theta);
b = r .* sin(theta);
c = Pr .* cos(theta) - (Ptheta ./ r) .* sin(theta);
d = Pr .* sin(theta) + (Ptheta ./ r) .* cos(theta);
y0 = [a b c d]; % Initial Values for x, y, Px, Py
[x, y] = ode45(@ODEdbb,tspan,y0);
[m,n] = meshgrid(-3:.05:3,-3:.05:3);
contour(m,n,m.^3-3*n.^2.*m);
hold on
viscircles([centreX0,centreY0,],r);
axis square
% plot solution of ODE
x_ode = y(:,1);
y_ode = y(:,2);
plot(x_ode,y_ode)
axis([-3 3 -3 3])
hold off
function f = ODEdbb(A, B)
x=B(1);
y=B(2);
Px=B(3);
Py=B(4);
%Differential Equations
dPxdt = 3 * y^2 -3 * x^2;
dPydt = 6 * x * y;
dxdt = Px;
dydt = Py;
f = [dxdt; dydt; dPxdt; dPydt];
end
function[value,isterminal,direction] = MonkeyPlot(A,Y)
value = (Y(1)^2 + Y(2)^2) - r^2;
isterminal = 1;
direction = 1;
end

Answers (1)

KSSV
KSSV on 11 Dec 2021
clc
tspan = [0:0.05:80]; %Range for independent var. i.e. time
r = 2;
centreX0 = 0;
centreY0 = 0;
theta = 1.08;
Ptheta = 3;
Pr = -sqrt(2 - 2 * r.^3 .* cos(theta).^3 + 6 .* r.^3 .* cos(theta) .* sin(theta).^2 - (Ptheta.^2 ./ r.^2));
%options = odeset('Events',@MonkeyPlot);
% y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
a = r .* cos(theta);
b = r .* sin(theta);
c = Pr .* cos(theta) - (Ptheta ./ r) .* sin(theta);
d = Pr .* sin(theta) + (Ptheta ./ r) .* cos(theta);
y0 = [a b c d]; % Initial Values for x, y, Px, Py
[x, y] = ode45(@ODEdbb,tspan,y0);
Warning: Failure at t=3.877326e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.
% Circle
viscircles([centreX0,centreY0],r);
hold on
th = linspace(0,2*pi) ;
xc = centreX0+r*cos(th) ;
yc = centreY0+r*sin(th) ;
[m,n] = meshgrid(-3:.05:3,-3:.05:3);
% Get points insode circle
idx = inpolygon(m,n,xc,yc) ;
mn = m.^3-3*n.^2.*m ;
mn(~idx) = NaN ;
contour(m,n,mn);
axis square
% plot solution of ODE
x_ode = y(:,1);
y_ode = y(:,2);
plot(x_ode,y_ode)
axis([-3 3 -3 3])
hold off
function f = ODEdbb(A, B)
x=B(1);
y=B(2);
Px=B(3);
Py=B(4);
%Differential Equations
dPxdt = 3 * y^2 -3 * x^2;
dPydt = 6 * x * y;
dxdt = Px;
dydt = Py;
f = [dxdt; dydt; dPxdt; dPydt];
end

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!