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일
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 CenterFile 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