Plots giving me straight lines instead of harmonic oscillation lines

clc; clear;

%Variables
n = 10000; %number of iterations
m = 1; %attached mass(kg)
k = 10; %spring constant(N/m)
y_0 = 0.2; %initial distance stretched downward (m)
v_0 = 0; %initial velocity (m/s)
b = 0; %damping coefficient
F = 0; %force applied
w = 0; %frequency

%Arrays
t = zeros(1,n);
y = zeros(1,n);
v = zeros(1,n);
a = zeros(1,n);
y_ec = zeros(1,n);
v_ec = zeros(1,n);
a_ec = zeros(1,n);

t(1) = 0; %starting time
y(1) = y_0; %starting vertical position for Euler method
v(1) = v_0; %starting velocity Euler method
a(1) = -k*y(1)/m; %starting acceleration Euler method
y_ec(1) = y_0; %starting vertical position for Euler-Cromer method
v_ec(1) = v_0; %starting velocity for Euler-Cromer method

%Loop
for dt = [0.5 0.1 0.01 0.001] %using 4 different step sizes to find the most accurate one
    for i=2:n
    %Euler Method
    t(i) = t(i-1)*dt;
    y(i) = y(i-1)+v(i-1)*dt;
    v(i) = v(i-1)+a(i-1)*dt;
    a(i) = -k/m*y(i-1)-b/m*v(i-1)+F/m*cos(w*t(i-1));
    
    %Euler-Cromer Method
    t_ec(i) = t(i-1)*dt;
    y_ec(i) = y_ec(i-1)+v_ec(i)*dt;
    v_ec(i) = v_ec(i-1)+a_ec(i-1)*dt;
    a_ec(i) = -k/m*y_ec(i-1)-b/m*v_ec(i-1)+F/m*cos(w*t(i-1));
    
    %exact functions for position and velocity
    y_exact(i) = y(1)*cos(sqrt(k/m)*t(i));
    v_exact(i) = -sqrt(k/m)*y(1)*sin(sqrt(k/m)*t(i));
    end
    hold on
    plot(t,y);
    plot(t,y_ec)
end
    
%exact functions for position and velocity
y_exact(i) = y(1)*cos(sqrt(k/m)*t(i));
v_exact(i) = -sqrt(k/m)*y(1)*sin(sqrt(k/m)*t(i));

%plotting Euler and Euler-Cromer methods vs. exact
t = [0.001:0.001:0.001:10000]
plot(t,y_exact,'r');
xlim([0 3]);
ylim([-0.8 0.8]);
xlabel('Time(s)');
ylabel('Position(m)');
title('Comparing Position vs. Time of Exact, Euler, and Euler-Cromer Methods');
legend 

I keep trying to graph these curves. They should look like a harmonic oscillation wave but for some reason the lines are all straight and not even in the right place (shown in the plot above). Any idea why this is happening?

Answer 1

[Mathieu NOE]

Hello

there were some bugs ... where to start ?

first the way you order the recursive equations can have some importance

usually you update first the acceleration , then the velocity, then the position

for the basic Euler method it's not so critical because all equations are based on past values, but the Euler -Cromer Method update the position based on the updated velocity , so it is important in this case that you update first the velocity and then the position (and not vice versa)

then there are some minor issues related to initialisation

time step should be defined in relation to system's natural frequency : it always should be a small fraction of the system's natural period; of course , the smaller the step , the closer to the "true" response you'll get; there is no need for a for loop to optimize as the optimum in terms of accuracy will always be a smaller time step; or use a higher order solver (RK4);

also , what you call "exact solution" is certainly not the exact solution of this second order differential equation; I leave it in the plot but it's not by any way the "exact" solution

here is the code updated :

clc; clearvars;
%Variables
n = 1000; %number of iterations
m = 1; %attached mass(kg)
k = 10; %spring constant(N/m)
y_0 = 0.2; %initial distance stretched downward (m)
v_0 = 0; %initial velocity (m/s)
b = 1; %damping coefficient
F = 1; %force applied
% w = 0; %frequency
w = sqrt(k/m);
dt_max = 0.01/(w/(2*pi));
dt_min = 0.01*dt_max;

%Arrays
t = zeros(1,n);
y = zeros(1,n);
v = zeros(1,n);
a = zeros(1,n);
y_ec = zeros(1,n);
v_ec = zeros(1,n);
a_ec = zeros(1,n);
t(1) = 0; %starting time
y(1) = y_0; %starting vertical position for Euler method
v(1) = v_0; %starting velocity Euler method
a(1) = -k*y(1)/m; %starting acceleration Euler method
y_ec(1) = y_0; %starting vertical position for Euler-Cromer method
v_ec(1) = v_0; %starting velocity for Euler-Cromer method
a_ec(1) = -k*y(1)/m; %starting acceleration Euler-Cromer  method

y_exact(1) = y(1);

%Loop
dtt = linspace(dt_min,dt_max,4);
for dt = dtt % [0.001] %using 4 different step sizes to find the most accurate one
    for i=2:n
        %Euler Method
        t(i) = t(i-1)+dt;
        a(i) = -k/m*y(i-1)-b/m*v(i-1)+F/m*cos(w*t(i-1));
        v(i) = v(i-1)+a(i-1)*dt;
        y(i) = y(i-1)+v(i-1)*dt;

        %Euler-Cromer Method
        %t_ec(i) = t(i-1)+dt; % not used anyway
        a_ec(i) = -k/m*y_ec(i-1)-b/m*v_ec(i-1)+F/m*cos(w*t(i-1));
        v_ec(i) = v_ec(i-1)+a_ec(i-1)*dt;
        y_ec(i) = y_ec(i-1)+v_ec(i)*dt;    

%         %exact functions for position and velocity
%         y_exact(i) = y(1)*cos(w*t(i));
%         v_exact(i) = -w*y(1)*sin(w*t(i));  
    end
    
    hold on
    plot(t,y,'bd');
    plot(t,y_ec,'r*');
%     plot(t,y_exact,'k');
    legend('Euler','Euler-Cromer','exact');
end
    hold off