Search code examples
pythonnumerical-methodsnumericalnumerical-integration

Explicit Euler method doesn't behave how I expect


I have implemented the following explicit euler method in python:

def explicit_euler(df, x0, h, N):
    """Solves an ODE IVP using the Explicit Euler method.

    Keyword arguments:
    df  - The derivative of the system you wish to solve.
    x0 - The initial value of the system you wish to solve.
    h  - The step size.
    N  - The number off steps.
    """
    x = np.zeros(N)
    x[0] = x0

    for i in range(0, N-1):
        x[i+1] = x[i] + h * df(x[i])

    return x

Following the article on wikipedia I can plot the function and verify that I get the same plot: Explicit euler. I believe that here the method I have written is working correctly.

Next I tried to use it to solve the last system given on this page and instead of the plot shown there I obtain this:

enter image description here

I am not sure why my plot doesn't match the one shown on the webpage. The explicit euler method seems to work fine when I use it to solve systems where the slope doesn't change, but for an oscillating function it never seems to mimic it at all. Not even showing the expected error gain as indicated on the linked webpage. I am not sure what is wrong with the method I have implemented.

Here is the code used for plotting and the derivative:

def g(t):
    return -0.5 * np.exp(t * 0.5) * np.sin(5 * t) + 5 * np.exp(t * 0.5) 
    * np.cos(5 * t)

h = 0.001
x0 = 0
tn = 4
N = int(tn / h)

x = ee.explicit_euler(f, x0, h, N)
t = np.arange(0, tn, h)

fig = plt.figure()
plt.plot(t, x, label="Explicit Euler")
plt.plot(t, (np.exp(0.5 * t) * np.sin(5 * t)), label="Analytical 
solution")
#plt.plot(t, np.exp(0.5 * t), label="Analytical solution")
plt.xlabel('Timesteps t')
plt.ylabel('x(t)=e^(0.5*t) * sin(5*t)')
plt.legend()
plt.grid()
plt.show()

Edit:

As requested here is the current equation I am applying the method to:

y'-y=-0.5*e^(t/2)*sin(5t)+5e^(t/2)*cos(5t)

Where y(0)=0.

I would like to make clear however that this behaviour doesn't occur just for this equation but all equations where the slope has a change in sign, or oscillating behaviour.

Edit 2: Ok thanks. Yes the code below does indeed work. But I have one further question. In the simple example I had for the exponential function, I had defined a method:

def f(x): 
    return x

for the system f'(x)=x. This gave the output of my first graph which looks correct. I then defined another function:

def k(x): 
    return cos(x)

for the system f'(x)=cos(x), this does not give expected output. But when I change the function definition to

def k(t, x): 
    return cos(t) 

I get the expected output. If I change my function

def f(t, x): 
    return t 

I get an incorrect output. Am I always actually evaluating the function at a time step and is it just by chance for the system x'=x that at each time step the value is just the value of x?

I had understood that the Euler method used the value of the previously calculated value in order to get the next value. But if I run code for my function k(x)=cos(x), I get output pictured below, which must be incorrect. This now uses the updated code you provided.

enter image description here

def k(t, x):
    return np.cos(x)

h  = 0.1         # Step size
x0 = (0, 0)        # Initial point of iteration
tn = 10        # Time step to iterate to
N  = int(tn / h)   # Number of steps

x = ee.explicit_euler(k, x0, h, N)
t = np.arange(0, tn, h)

Solution

  • The problem is that you have incorrectly raised the function g, you want to solve the equation:

    enter image description here

    From where we observe that:

    y' = y -0.5*e^(t/2)*sin(5t)+5e^(t/2)*cos(5t)
    

    Then we define the function f(t, y) = y -0.5*e^(t/2)*sin(5t)+5e^(t/2)*cos(5t) as:

    def f(t, y):
        return y -0.5 * np.exp(t * 0.5) * np.sin(5 * t) + 5 * np.exp(t * 0.5) * np.cos(5 * t) 
    

    The initial point of iteration is f0=(t(0), y(0)):

    f0 = (0, 0)
    

    Then from Euler's equations:

    def explicit_euler(df, x0, h, N):
        """Solves an ODE IVP using the Explicit Euler method.
    
        Keyword arguments:
        df  - The derivative of the system you wish to solve.
        x0 - The initial value of the system you wish to solve.
        h  - The step size.
        N  - The number off steps.
        """
        x = np.zeros(N)
        t, x[0] = x0
    
        for i in range(0, N-1):
            x[i+1] = x[i] + h * df(t ,x[i])
            t += h
    
        return x
    

    Complete Code:

    def explicit_euler(df, x0, h, N):
        """Solves an ODE IVP using the Explicit Euler method.
    
        Keyword arguments:
        df  - The derivative of the system you wish to solve.
        x0 - The initial value of the system you wish to solve.
        h  - The step size.
        N  - The number off steps.
        """
        x = np.zeros(N)
        t, x[0] = x0
    
        for i in range(0, N-1):
            x[i+1] = x[i] + h * df(t ,x[i])
            t += h
    
        return x
    
    def df(t, y):
        return -0.5 * np.exp(t * 0.5) * np.sin(5 * t) + 5 * np.exp(t * 0.5) * np.cos(5 * t) + y
    
    h = 0.001
    f0 = (0, 0)
    tn = 4
    N = int(tn / h)
    
    x = explicit_euler(df, f0, h, N)
    t = np.arange(0, tn, h)
    
    fig = plt.figure()
    plt.plot(t, x, label="Explicit Euler")
    plt.plot(t, (np.exp(0.5 * t) * np.sin(5 * t)), label="Analytical solution")
    #plt.plot(t, np.exp(0.5 * t), label="Analytical solution")
    plt.xlabel('Timesteps t')
    plt.ylabel('x(t)=e^(0.5*t) * sin(5*t)')
    plt.legend()
    plt.grid()
    plt.show()
    

    Screenshot:

    enter image description here

    Dump y' and what is on the right side is what you should place in the df function.

    We will modify the variables to maintain the same standard for the variables, and will y be the dependent variable, and t the independent variable.

    Equation 2: In this case the equation f'(x)=cos(x) will be rewritten to:

    y'=cos(t)
    

    Then:

    def df(t, y):
        return np.cos(t)
    

    enter image description here

    In conclusion, if we have an equation of the following form:

    y' = f(t, y)
    

    Then:

    def df(t, y):
        return f(t, y)