Plots giving me straight lines instead of harmonic oscillation lines
4 views (last 30 days)
Show older comments
Frank Kostecki
on 22 Apr 2021
Answered: Mathieu NOE
on 22 Apr 2021
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/592340/image.png)
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?
0 Comments
Accepted Answer
Mathieu NOE
on 22 Apr 2021
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
0 Comments
More Answers (0)
See Also
Categories
Find more on 2-D and 3-D Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!