MATLAB Answers

Why does my code only print straight lines?

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

Sign in to comment.

Answers (1)

DGM
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)');
  1 Comment
Colton Norfleet
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!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!