python-3.xjupyter-notebookodeeulers-number

Plotting solution 2nd ODE using Euler


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)

Solution

  • 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

    plot of exact solution

    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

    enter image description here

    with the expected increasing amplitude.