Why isn't my graph plotting correctly?

1 view (last 30 days)
Danny Brett
Danny Brett on 29 Nov 2021
Commented: Mathieu NOE on 30 Nov 2021
I am trying to get a plot ontop of this contour plot but while staying within the circle (ie. starting on the circle and finishing once it hits the circle again).
As you can see below it is not doing quite what I was hoping for, any help would be greatly appreciated.
function f = ODEms(A, Y)
x=Y(1);
y=Y(2);
Px=Y(3);
Py=Y(4);
%Differential Equations
dPxdt = 3 * y^2 -3 * x^2;
dPydt = 6 * x * y;
dxdt = Px;
dydt = Py;
f = [dPxdt; dPydt; dxdt; dydt];
end
clc
xspan = [0 80] %Range for independent var. i.e. time
y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
options = odeset('Events',@MonkeyPlot);
[x y] = ode45(@ODEms,xspan,y0,options);
[x0,y0] = meshgrid(-2:.05:2,-2:.05:2);
contourf(x0.^3-3*y0.^2.*x0);
radius = 30;
centreX0 = 40;
centreY0 = 40;
viscircles([centreX0,centreY0,],radius);
axis square
hold on
plot(y(:,1),y(:,2));
legend('Px','Py')
ylabel('y');
xlabel('x');
axis([0 80 0 80]);
title('Initial Energy and Kinetic Energy of the Particle')
hold off
function[value,isterminal,direction] = MonkeyPlot(A,Y)
value = (Y(1)^2 + Y(2)^2) - 30^2;
isterminal = 1;
direction = 1;
end
  2 Comments
Danny Brett
Danny Brett on 30 Nov 2021
This is currently the ODE with the initial conditions as stated in the code.
I however want the plot to start on the edge of the circle, plot a path inside the circle and at the moment the plot hits the circle again to stop.
So (please excuse my quick paint job) it would hopefully look something like this.

Sign in to comment.

Answers (1)

Mathieu NOE
Mathieu NOE on 30 Nov 2021
hello
I have a bit modified your code , I had to replace viscircles with a similar function (cercle ) - but you can easily come back to your own viscircles implementation
I had to understand that there are sometimes a factor 1000 between some distance units / arrays to make the data and the plot consistent
I also coded the start point of the inner line based on the radius and center value of the red circle to make the code more robust
now, there is still a bit of work to make the ODE create more (refined) output points so that the inner line does not stop so far away for the circle - but I didn't want to mess up too much here - not sure exactly to understand what MonkeyPlot is here for
so far , this is the result I could obtain :
and the code :
clc
xspan = [0 80]; %Range for independent var. i.e. time
radius = 30;
centreX0 = 40;
centreY0 = 40;
% y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
a = (-radius/sqrt(2)+centreX0)/1000;
b = (-radius/sqrt(2)+centreY0)/1000;
y0 = [a b a b] % Initial Values for x, y, Px, Py
options = odeset('Events',@MonkeyPlot);
[x, y] = ode45(@ODEms,xspan,y0,options);
[x0,y0] = meshgrid(-2:.05:2,-2:.05:2);
contourf(x0.^3-3*y0.^2.*x0);
hold on
radius = 30;
centreX0 = 40;
centreY0 = 40;
% viscircles([centreX0,centreY0,],radius); % replaced by next two lines
[Xc,Yc] = cercle([centreX0,centreY0],radius,360);
plot(Xc,Yc,'r');
axis square
% plot solution of ODE
x_ode = y(:,1)*1000;
y_ode = y(:,2)*1000;
% remove data points outside the circle
ind = find(sqrt((x_ode-centreX0).^2+(y_ode-centreY0).^2) > radius);
x_ode(ind) = [];
y_ode(ind) = [];
plot(x_ode,y_ode);
legend('Px','Py')
ylabel('y');
xlabel('x');
axis([0 80 0 80]);
title('Initial Energy and Kinetic Energy of the Particle')
hold off
function [XData,YData] = cercle(centre,rayon,points)
if(nargin < 2)
points = 50;
end
theta = 0:2*pi/(points-1):2*pi;
XData = centre(1)+rayon.*cos(theta);
YData = centre(2)+rayon.*sin(theta);
if(nargout < 2)
XData = plot(XData,YData);
end;
end
function f = ODEms(A, Y)
x=Y(1);
y=Y(2);
Px=Y(3);
Py=Y(4);
%Differential Equations
dPxdt = 3 * y^2 -3 * x^2;
dPydt = 6 * x * y;
dxdt = Px;
dydt = Py;
f = [dPxdt; dPydt; dxdt; dydt];
end
function[value,isterminal,direction] = MonkeyPlot(A,Y)
value = (Y(1)^2 + Y(2)^2) - 30^2;
isterminal = 1;
direction = 1;
end
  1 Comment
Mathieu NOE
Mathieu NOE on 30 Nov 2021
I finally decided to modify the xspan definition so I get more points for the inner line
result is now :
code :
clc
xspan = [0:0.05:80]; %Range for independent var. i.e. time
radius = 30;
centreX0 = 40;
centreY0 = 40;
% y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
a = (-radius/sqrt(2)+centreX0)/1000;
b = (-radius/sqrt(2)+centreY0)/1000;
y0 = [a b a b]; % Initial Values for x, y, Px, Py
[x, y] = ode45(@ODEms,xspan,y0);
[x0,y0] = meshgrid(-2:.05:2,-2:.05:2);
contourf(x0.^3-3*y0.^2.*x0);
hold on
radius = 30;
centreX0 = 40;
centreY0 = 40;
% viscircles([centreX0,centreY0,],radius); % replaced by next two lines
[Xc,Yc] = cercle([centreX0,centreY0],radius,360);
plot(Xc,Yc,'r');
axis square
% plot solution of ODE
x_ode = y(:,1)*1000;
y_ode = y(:,2)*1000;
% remove data points outside the circle
ind = find(sqrt((x_ode-centreX0).^2+(y_ode-centreY0).^2) > radius);
x_ode(ind) = [];
y_ode(ind) = [];
plot(x_ode,y_ode);
legend('Px','Py')
ylabel('y');
xlabel('x');
axis([0 80 0 80]);
title('Initial Energy and Kinetic Energy of the Particle')
hold off
function [XData,YData] = cercle(centre,rayon,points)
if(nargin < 2)
points = 50;
end
theta = 0:2*pi/(points-1):2*pi;
XData = centre(1)+rayon.*cos(theta);
YData = centre(2)+rayon.*sin(theta);
if(nargout < 2)
XData = plot(XData,YData);
end;
end
function f = ODEms(A, Y)
x=Y(1);
y=Y(2);
Px=Y(3);
Py=Y(4);
%Differential Equations
dPxdt = 3 * y^2 -3 * x^2;
dPydt = 6 * x * y;
dxdt = Px;
dydt = Py;
f = [dPxdt; dPydt; dxdt; dydt];
end
% function[value,isterminal,direction] = MonkeyPlot(A,Y)
%
% value = (Y(1)^2 + Y(2)^2) - 30^2;
% isterminal = 1;
% direction = 1;
%
% end

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!