Search code examples
python-3.xnumerical-methodsdifferential-equationsnumerical-analysisnumerical-stability

Euler Method implementation in Python gives a stable result but it should be unstable

I am trying to solve this differential equation with the Euler Method using Python3:

enter image description here

According to Wolfram Alpha, that's the plot of the correct equation.

enter image description here

Again, according to Wolfram Alpha, in this case, the classic Euler Method should not be stable, as you can see by the end of the interval:

enter image description here

However, on my implementation, Euler method provides a stable result, which is strange. I wonder that my implementation is wrong for some reason. Nonetheless, I can't find the error.

I generated some points and a plot comparing my approximation and the analytic output of the function. In blue, the analytic result as control group. In red, the output of my implementation:

enter image description here

That's my code:

import math
import numpy as np
from matplotlib import pyplot as plt
import pylab

def f(x):

    return (math.e)**(-10*x)

def euler(x):

    y_init = 1
    x_init = 0

    old_dy_dx = -10*y_init

    old_y = y_init 

    new_y = None

    new_dy_dx = None

    delta_x = 0.001

    limite = 0

    while x>limite:

        #for i in range(1,6):

        new_y = delta_x*old_dy_dx + old_y
        #print ("new_y", new_y)

        new_dy_dx = -10*new_y
        #print ("new dy_dx", new_dy_dx)

        old_y = new_y
        #print ("old_y", old_y)

        old_dy_dx = new_dy_dx
        #print ("old delta y_delta x", old_dy_dx)
        #print ("iterada",i)

        limite = limite +delta_x

    return new_y

t = np.linspace(-1,5, 80)

lista_outputs = []

for i in t:
    lista_outputs.append(euler(i))
    print (i)

# red dashes, blue squares and green triangles
plt.plot(t, f(t), 'b-', label='Output resultado analítico')
plt.plot(t , lista_outputs, 'ro', label="Output resultado numérico")
plt.title('Comparação Euler/Analítico - tolerância: 0.3')
pylab.legend(loc='upper left')
plt.show()

Thanks for the help.

============================================================

UPDATE

With @SourabhBhat help, I was able to see that my implementation was, actually, right. It was, indeed, generating an instability. Besides increasing the step size, I needed to do some zoom in to see it happening.

The picture bellow speaks for itself (step size of 0.22):

enter image description here

Solution

  • Depending on the time step size, the Euler integration can be stable or unstable, since it is an explicit method. You have chosen a pretty small time step. If you increase it you will start seeing the oscillations, as shown in figure below.

    Oscillations

    Here is a small testing program which I wrote (try slowly increasing the steps variable [20,30,40,50....]):

    import numpy as np
import matplotlib.pyplot as plt

steps = 20


def exact_solution(t):
    return np.exp(-10.0 * t)


def numerical_solution(y0, dt, num_steps):
    y = np.zeros(num_steps + 1)
    y[0] = y0
    for step in range(num_steps):
        y[step + 1] = y[step] - 10.0 * y[step] * dt

    return y


if __name__ == "__main__":
    t0 = 0
    time = np.linspace(t0, 5, steps + 1)
    num_sol = numerical_solution(exact_solution(t0), time[1] - time[0], steps)
    exact_sol = exact_solution(time)

    plt.plot(time, num_sol, ".--", label="numerical")
    plt.plot(time, exact_sol, label="exact")
    plt.legend(loc="best")
    plt.show()