Event Detection on a circle

2 次查看（过去 30 天）
Danny Brett
Danny Brett 2021-12-10
回答： KSSV 2021-12-11
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

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

回答（1 个）

KSSV
KSSV 2021-12-11
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by