
import numpy as np
import matplotlib.pyplot as plt


def f(y): return -y**2  # differential equation
def Df(y): return -2*y

def g(f,y,h,ynext):  # Backward Euler objective function: we want to choose ynext to make this zero
    return y + h*f(ynext) - ynext

def Dg(f,y,h,ynext):
    return h*Df(ynext) - 1


t0,t1 = 0,20

tt = np.linspace(t0,t1,200)
yy = 1/(tt+1) # exact solution y = 1/(t+1)
plt.plot(tt,yy)

n = 5
h = (t1-t0)/n

t = 0
y = 1
tol = 1.e-6

ts = [t]
ys = [y]
for j in range(n):
    ynext = y # starting guess for iterative solving of implicit formula
    
    # Newton iteration
    while True:
        step = -g(f,y,h,ynext)/Dg(f,y,h,ynext)
        ynext = ynext + step
        if np.abs(step) < tol:
            y = ynext
            t += h
            ys.append(y)
            ts.append(t)
            break
plt.plot(ts,ys,'o',ms=5,alpha=0.5,label='Backward Euler');
            
# Forward Euler for comparison
t = 0
y = 1
ts = [t]
ys = [y]
for j in range(n):
    y = y + h*f(y)
    t += h
    if np.abs(y)>3: break
    ys.append(y)
    ts.append(t)
    
    
plt.plot(ts,ys,'o',ms=5,alpha=0.5,label='Forward Euler')
plt.legend(); plt.xlabel('t'),plt.ylabel('y',rotation=0)
plt.title(f'h={h}')
plt.savefig(f'backward_euler_h{h}.png',bbox_inches='tight')

Example of using an implicit LMF¶