#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon May 27 13:31:24 2019

@author: kerkmann
"""

import numpy as np
import matplotlib.pyplot as plt

# explicit Euler method (from exercise 27)
def Euler_explicit(x,t,f,k):
    return x + k*f(x,t)

# midpoint rule (1 Punkt)
def midpoint_rule(x,xm1,t,f,k):
    return xm1 + 2*k*f(x,t)

### MAIN PROGRAM ###
    
# Consider u'(t) = 2*u(t), u(0) = 1
f = lambda u,t: 2*u
ex = lambda t: np.exp(2*t)

err = np.zeros(5)
for j in range(5):
    spaces = 100*2**j
    t = np.linspace(0,1,spaces+1)
    k = t[1]-t[0]
    u = np.zeros_like(t)

    # (1 Punkt)
    u[0] = 1
    u[1] = Euler_explicit(u[0],0,f,k)
    
    # (2 Punkte)
    for n in np.arange(2,len(t)):
        u[n] = midpoint_rule(u[n-1],u[n-2],t[n],f,k)
        
    plt.plot(t,u)

    # (1 Punkt)
    err[j] = np.max(abs(ex(t) - u))

# (1 Punkt)
print('Fehler:')
print(err)
EOC = np.log(err[:-1]/err[1:])/np.log(2)
print('Experimentelle Konvergenzordnung:')
print(EOC)