# Why does my code only print straight lines?

1 view (last 30 days)
Colton Norfleet on 23 Sep 2021
Edited: Colton Norfleet on 23 Sep 2021
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 on 23 Sep 2021
Comment the old plot line outside for loop

DGM on 23 Sep 2021
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 on 23 Sep 2021
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!