Crazy MATLAB Graph: Solving Simple Pendulum Eqns

[Chairman Lmao]

Crazy MATLAB graph!

Hi!
I had been trying to numerically solve the simple pendulum equations using RK4 method in matlab. I had plotted the energy versus time graph. One would expect it to be just a straight line parallel to time axis. However, we got a weird linearly decreasing graph with periodic fluctuations, though the range in which the values were fluctuating was very small. Could anyone please explain this?
http://yfrog.com/89matlabj<-link for the graph(Energy vs time in s). The initial angle was 90deg and initial angular momentum 0.

If it helps, this was my code:

function z = pro1a(x,y)
q=zeros(1,6000);
p=zeros(1,6000);
t=zeros(1,6000);
E=zeros(1,6000);
q(1)=x;
p(1)=y;
t(1)=0;
E(1)=y.^2/2-9.8*cos(x);
for i=2:6000
k1=p(i-1);
j1=-9.8*sin(q(i-1));
k2=p(i-1)+0.005*j1;
j2=-9.8*sin(q(i-1)+0.005*k1);
k3=p(i-1)+0.005*j2;
j3=-9.8*sin(q(i-1)+0.005*k2);
k4=p(i-1)+0.01*j3;
j4=-9.8*sin(q(i-1)+0.01*k3);
q(i)=mod(q(i-1)+0.01*(k1+2*k2+2*k3+k4)/6,4*pi);
if q(i)>2*pi
q(i)=q(i)-4*pi;
end
p(i)=p(i-1)+0.01*(j1+2*j2+2*j3+j4)/6;
t(i)=t(i-1)+0.01;
E(i)=p(i).*p(i)./2-9.8.*cos(q(i));
end
plot(t,E);
z=max(E)-min(E)
end

x=initial angle
y=initial angular momentum
q=angle as function of time
p=angular momentum as function of time

Also, max(E)-min(E) was 4.893*10^(-7)

Thank you for your consideration! :)

 

[Valerie Park]

It looks like the graph you are seeing is due to numerical error from using the Runge-Kutta method. The Runge-Kutta method is an approximate solution to the differential equations, meaning that it will not give an exact value (which would be a straight line parallel to the time axis). The periodic fluctuations you are seeing are probably due to the small stepsize used in the calculations. As the stepsize decreases, the numerical error decreases, resulting in a more accurate graph.