Why isn't my graph plotting correctly?

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

[Dyuman Joshi]

What is the result you are getting? And what are you hoping to achieve?

[Danny Brett]

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.

Answer 1

[Mathieu NOE]

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

[Mathieu NOE]

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