Potential Energy Function Identically Zero

Question

I am using ode45 to numerically solve the gravitational system of a sun an some planet. Defined the function:

function eqs = prob4_4(t,in2,in3)
%PROB4_4
%    EQS = PROB4_4(T,IN2,IN3)

%    This function was generated by the Symbolic Math Toolbox version 5.11.
%    05-Mar-2014 15:56:04

G = in3(:,2);
M = in3(:,1);
vx = in2(3,:);
vy = in2(4,:);
x = in2(1,:);
y = in2(2,:);
t2 = x.^2;
t3 = y.^2;
t4 = t2+t3;
t5 = 1.0./t4.^(3.0./2.0)
eqs = [vx;vy;-G.*M.*t5.*x;-G.*M.*t5.*y];

I then used ode45 to generate x and y positions, and velocities in the x and y direction, as a function of time:

[tout,yout] = ode45(@(t,Y) lab6(t,Y,[1,4*pi^2]),[0,1],[1;0;0;5],myoptions);

Using the x positions, yout(:,1), and the y positions, yout(:,2), I calculated the potential energy of this two body system:

PE = 4*pi^2/sqrt((yout(:,1)).^2 + (yout(:,2)).^2);

When I went to plot PE versus tout, I got a horizontal line across the t-axis. Looking at the elements of PE, I found that they were all zero. Why did this happen?

When I went to calculate the kinetic energy, I defined the variable KE as

KE = 1/2*(yout(:,3).^2 + yout(:,4).^2);

I plotted this as a function of time and got positive periodic kinetic energy behavior, which is what I expected; I expected a similar plot for PE, but with negative values, because the gravitational system is bound.

Could anyone help?

Answer 1

I figured it out. I made one little matlab mistake, and one little physics mistake. The former was that I forgot to use the element by element division operator, namely, ./ So, the potential energy function becomes PE = -4*pi^2 ---> ./ <---- sqrt((yout(:,1)).^2 + (yout(:,2)).^2); and the latter mistake was that gravitational potential energy values are measured as negative values, so I had to put in a negative symbol.