Plotting an ODE with an event detection function

2 views (last 30 days)
Danny Brett
Danny Brett on 23 Apr 2022
Answered: Walter Roberson on 23 Apr 2022
I am trying to get a plot of an ODE travelling over this surface (the contour plot and circle aren't important but necessary) and I am trying to find the periodic points.
When I am trying to plot though I just get the contour plot and the circle, I don't get the line the ODE creates, any help would be great.
**You should see the line start on the x-axis go up slightly and then come back down on itself, and stop again once it hits the x-axis
clc
tspan = [0:0.05:80]; %Range for independent var. i.e. time
r = 2;
centreX0 = 0;
centreY0 = 0;
y = 0;
x = -0.31417;
Px = 0;
%Total Energy = Kinetic Energy + Potential Energy
%KE = (1 ./ 2) .* (Px.^2 + Py.^2);
%V = (x.^3 - 3 .* x .* y^2);
%Total Energy = 1
Py = sqrt((1 - x.^3) .* 0.5);
y0 = [x y Px Py];
option = @isYnegative;
[x, y] = ode45(@ODEdbb,tspan,y0,option);
[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
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 = isYnegative(Y)
value = Y(2) == 0;
isterminal = 1;
direction = -1;
end
  1 Comment
Torsten
Torsten on 23 Apr 2022
Since you don't refer to the ODE solution in your plot, it cannot appear therein.

Sign in to comment.

Answers (1)

Walter Roberson
Walter Roberson on 23 Apr 2022
y0 = [x y Px Py];
[x, y] = ode45(@ODEdbb,tspan,y0,option);
The first output from ode45() is the time values. The second output from ode45() is an array that is the value of the boundary conditions corresponding to each time.
Your code appears to assume that because you are thinking of x and y as your first two input coordinates, that the first two outputs are x and y. That is not correct at all. You should be using something like
[t, out] = ode45(@ODEdbb,tspan,y0,option);
x = out(:,1);
y = out(:,2);
after which you can
plot(x, y)

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!