Euler's method in python

Question

I'm trying to implement euler's method to approximate the value of e in python. This is what I have so far:

def Euler(f, t0, y0, h, N):
    t = t0 + arange(N+1)*h
    y = zeros(N+1)
    y[0] = y0
    for n in range(N):
        y[n+1] = y[n] + h*f(t[n], y[n])
        f = (1+(1/N))^N
    return y

However, when I try to call the function, I get the error "ValueError: shape <= 0". I suspect this has something to do with how I defined f? I tried inputting f directly when euler is called, but gave me errors related to variables not being defined. I also tried defining f as its own function, which gave me a division by 0 error.

def f(N):
    for n in range(N): 
        return (1+(1/n))^n

(not sure if N was the appropriate variable to use here...)

Answer 1

The formula you are trying to use is not Euler's method, but rather the exact value of e as n approaches infinity wiki,

$n = \lim_{n\to\infty} (1 + \frac{1}{n})^n$

---

Euler's method is used to solve first order differential equations.

Here are two guides that show how to implement Euler's method to solve a simple test function: beginner's guide and numerical ODE guide.

To answer the title of this post, rather than the question you are asking, I've used Euler's method to solve usual exponential decay:

$\frac{dN}{dt} = -\lambda N$

Which has the solution,

$N(t) = N_0 e^{-\lambda t}$

Code:

import numpy as np
import matplotlib.pyplot as plt
from __future__ import division

# Concentration over time
N = lambda t: N0 * np.exp(-k * t)
# dN/dt
def dx_dt(x):
    return -k * x

k = .5
h = 0.001
N0 = 100.

t = np.arange(0, 10, h)
y = np.zeros(len(t))

y[0] = N0
for i in range(1, len(t)):
    # Euler's method
    y[i] = y[i-1] + dx_dt(y[i-1]) * h

max_error = abs(y-N(t)).max()
print 'Max difference between the exact solution and Euler's approximation with step size h=0.001:'

print '{0:.15}'.format(max_error)

Output:

Max difference between the exact solution and Euler's approximation with step size h=0.001:
0.00919890254720457

---

Note: I'm not sure how to get LaTeX displaying properly.

Answer 2

Are you sure you are not trying to implement the Newton's method? Because Newton's method is used to approximate the roots.

In case you decide to go with Newton's method, here is a slightly changed version of your code that approximates the square-root of 2. You can change f(x) and fp(x) with the function and its derivative you use in your approximation to the thing you want.

import numpy as np


def f(x):
    return x**2 - 2


def fp(x):
    return 2 * x


def Newton(f, y0, N):
    y = np.zeros(N + 1)
    y[0] = y0
    for n in range(N):
        y[n + 1] = y[n] - f(y[n]) / fp(y[n])
    return y


print(Newton(f, 1, 10))

gives

[ 1. 1.5 1.41666667 1.41421569 1.41421356 1.41421356 1.41421356 1.41421356 1.41421356 1.41421356 1.41421356]

which are the initial value and the first ten iterations to the square-root of two.

Besides this a big problem was the usage of ^ instead of ** for powers which is a legal but a totally different (bitwise) operation in python.