Why does my code only print straight lines?

So I am doing some code for my mechanics class. He typically gives us code to add to or edit. The code that was given was for a spring-mass problem using equations derived from position depedent force, that has the solved analytic equations done for us. Which plots a nice mechanical energy plot with kinetic, potential, force, and total energy. My problem arrises with my portion of the code, I know from my workspace that despite me using the Euler-Cromer numerical approach, my numbers match almost exactly to the analytics, and they should. However, the plotting for my portion of the code just comes up as flat lines. I can't not find any reasoning for it. It is almost as if it's assuming its not a function at all, despite having "y" and "x" values.

%Original code in fofx is for the analytical method

clear;

m=1.0; %mass

k=0.01; %spring constant

w=sqrt(k/m); %natural frequency

x0=0.0; %initial position

v0=0.5; %initial velocity

t=[0:0.05:2*pi/w]; %time array from zero to one oscillation period

E0=0.5*m*v0^2; %total initial energy

x=sqrt(2*E0/k)*sin(w.*t); %position versus time array

v=v0*cos(w.*t); %velocity versus time arry

%a=-k*x/m; %acceleration versus time array if needed

PE=0.5*k*x.^2; %potential energy array

KE=0.5*m*v.^2; %kinetic energy array

E=PE+KE; %total energy array

F=-k*x; %force array

plot(x,PE,'k-',x,KE,'b:',x,E,'r-.',x,F,'m--');

%Below I added the Numerical Approach

NPTS=length(t);

V(1)=v0;

X(1)=0;

%Tmax=2*pi/w;

%dt=Tmax/NPTS;

%format longE %stops the program from rounding up until the 15th decimal place

dt=0.05;

for i=1:NPTS-1

A(i)=(-k/m)*X(i);

V(i+1)=V(i)+A(i)*dt;

X(i+1)=X(i)+V(i+1)*dt;

pe(i)=0.5*k*(X(i))^2;

ke(i)=0.5*m*(V(i))^2;

f(i)=-k*(X(i));

T(i)=pe(i)+ke(i);

end

hold on

plot(X,pe(i),'k>',x,ke(i),'bo',x,T(i),'r+',x,f(i),'ms');

title('Spring-Mass Simple Harmonic Energy-Force Relation')

ylabel('PE, KE, E, F');

xlabel('x(m)');

[VBBV]

Add this line inside for loop and try

plot(X,pe(1:i),'k>',x,ke(1:i),'bo',x,T(1:i),'r+',x,f(1:i),'ms')

[Johan]

Hello,

it looks like your error is in the final plot() function you call. You are ploting the i element of the pe ke T and f arrays, remove the (i) and it should plot the results correctly.

[VBBV]

Comment the old plot line outside for loop

Answer 1

[DGM]

Not sure it's exactly how you want it, but there's this. I put them in subplots just for the sake of clarity. You don't need to do that.

%Original code in fofx is for the analytical method
clear;
m=1.0;                          %mass
k=0.01;                         %spring constant
w=sqrt(k/m);                    %natural frequency
x0=0.0;                         %initial position
v0=0.5;                         %initial velocity
t=[0:0.05:2*pi/w];              %time array from zero to one oscillation period
E0=0.5*m*v0^2;                  %total initial energy
x=sqrt(2*E0/k)*sin(w.*t);       %position versus time array
v=v0*cos(w.*t);                 %velocity versus time arry
%a=-k*x/m;                      %acceleration versus time array if needed
PE=0.5*k*x.^2;                  %potential energy array
KE=0.5*m*v.^2;                  %kinetic energy array
E=PE+KE;                        %total energy array
F=-k*x;                         %force array
subplot(2,1,1)
plot(x,PE,'k-',x,KE,'b:',x,E,'r-.',x,F,'m--');

%Below I added the Numerical Approach
NPTS=length(t);
V(1)=v0;
X(1)=0;
%Tmax=2*pi/w;
%dt=Tmax/NPTS;
%format longE    %stops the program from rounding up until the 15th decimal place
dt=0.05;
for i=1:NPTS-1 % sure this isn't supposed to be 1:NPTS to match vector lengths?
   A(i)=(-k/m)*X(i);
   V(i+1)=V(i)+A(i)*dt;
   X(i+1)=X(i)+V(i+1)*dt;
   pe(i)=0.5*k*(X(i))^2;
   ke(i)=0.5*m*(V(i))^2;
   f(i)=-k*(X(i));
   T(i)=pe(i)+ke(i);
end
X = X(1:end-1); % or do something else to avoid excess entry

subplot(2,1,2);
plot(X,pe,'k>',X,ke,'bo',X,T,'r+',X,f,'ms'); % using X instead of x?
title('Spring-Mass Simple Harmonic Energy-Force Relation')
ylabel('PE, KE, E, F');
xlabel('x(m)');

[Colton Norfleet]

so it appears I did overlook the bottom plot lines x and Xs, the for loop needs to be one less that the total number of points because itll calculate one more than the analytical would. Thank you for making a subplot, but because the professor I am doing this for is very particular he's requiring them to be on the same plot. I'm rather confused on how you managed to get it to work on a subplot. I tried something similar with figures but the same problem was presented. I am extremely perplexed as to why there even is a problem, if simiplying adding a subplot fixed the problem. If there a input limit that a plot can recieve?

Edit: I added your X = X(1:end-1); line and it fixed my problem and kept it on the same plot. Thank you for your help!