Plotting solution 2nd ODE using Euler

Question

I have used the Equation of Motion (Newtons Law) for a simple spring and mass scenario incorporating it into the given 2nd ODE equation y" + (k/m)x = 0; y(0) = 3; y'(0) = 0.

Using the Euler method and the exact solution to solve the problem, I have been able to run and receive some ok results. However, when I execute a plot of the results I get this diagonal line across the oscillating results that I am after.

Current plot output with diagonal line

Can anyone help point out what is causing this issue, and how I can fix it please?

MY CODE:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols
from sympy.abc import x, i
import math


# Given is y" + (k/m)x = 0; y(0) = 3; y'(0) = 0

# Parameters
h = 0.01;  #Step Size
t = 50.0;  #Time(sec)
k = 1;     #Spring Stiffness
m = 1;     #Mass
x0 = 3;
v0 = 0;

# Exact Analytical Solution
x_exact = x0*cos(math.sqrt(k/m)*t);
v_exact = -x0*math.sqrt(k/m)*sin(math.sqrt(k/m)*t);

# Eulers Method
x = np.zeros( int( t/h ) );
v = np.zeros( int( t/h ) );
x[1] = x0;
v[1] = v0;
x_exact = np.zeros( int( t/h ) );
v_exact = np.zeros( int( t/h ) );
te      = np.zeros( int( t/h ) );
x_exact[1] = x0;
v_exact[1] = v0;


#print(len(x));

for i in range(1, int(t/h) - 1):    #MAIN LOOP
    x[i+1] = x[i] + h*v[i];
    v[i+1] = v[i] - h*k/m*x[i];
    te[i]  = i * h
    x_exact[i] = x0*cos(math.sqrt(k/m)* te[i]);
    v_exact[i] = -x0*math.sqrt(k/m)*sin(math.sqrt(k/m)* te[i]);
#    print(x_exact[i], '\t'*2, x[i]);

#plot
%config InlineBackend.figure_format = 'svg'
plt.plot(te, x_exact, te ,v_exact)
plt.title("DISPLACEMENT")
plt.xlabel("Time (s)")
plt.ylabel("Displacement (m)")
plt.grid(linewidth=0.3)

You did not fill the allocated array to the end, so that you have trailing zeros that connect the end of the graph back to the origin. — Lutz Lehmann, Sep 20 '20 at 18:30

score 0 · Accepted Answer · answered Sep 20 '20 at 19:01

An in some details more direct computation is

te = np.arange(0,t,h)
N = len(te)

w = (k/m)**0.5
x_exact = x0*np.cos(w*te);
v_exact = -x0*w*np.sin(w*te);

plt.plot(te, x_exact, te ,v_exact)

resulting in

Note that arrays in python start at the index zero,

x = np.empty(N)
v = np.empty(N)
x[0] = x0;
v[0] = v0;

for i in range(N - 1):    #MAIN LOOP
    x[i+1] = x[i] + h*v[i];
    v[i+1] = v[i] - h*k/m*x[i];

plt.plot(te, x, te ,v)

then gives the plot

with the expected increasing amplitude.

Thank you for your enlightening description. I'm still getting the hang of utilizing python. You have helped me for sure! — C.L, Sep 21 '20 at 09:32

Plotting solution 2nd ODE using Euler

1 Answers1