#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 25 15:00:41 2019

@author: kerkmann
"""

import numpy as np

# Funktion aus (a) (4 Punkte)
def fd_centered(f,m,a,b,btype,order):
    '''Calculates the numerical solution to a given BVP
    using centered finite difference operator of 2nd order'''
    # Input parameters:
    # f       function for right hand side values
    # m       number of inner grid points
    # a       left boundary value (x=0)
    # b       right boundary value (x=1)
    # btype   dictionary with the keys 0 and 1 that indicates if the boundary
    #         conditions are Dirichlet (0) or Neumann (1) boundary conditions
    # order   used order of approximation to Neumann boundary conditions if
    #         they exist
    
    # Output parameters:
    # U   array with size m+2 containing (as output) the numerical solution
    
    x = np.linspace(0,1,m+2)
    h = 1/(m+1)
    
    # Matrix setup
    A = np.zeros((m+2,m+2))
    for i in range(1,m+1):
        A[i,i-1:i+2] = np.array([1, -2, 1])/h**2
       
    # left boundary (x=0)
    neumann_count = 0
    if btype[0] == 0: 
        # Dirichlet
        A[0,0] = 1
    else:
        # Neumann
        neumann_count += 1
        if order == 1:
            A[0,0:2] = np.array([-1,1])/h
        else:
            A[0,0:3] = np.array([-1.5, 2, -0.5])/h

    # right boundary (x=1)
    if btype[1] == 0:
        # Dirichlet
        A[m+1,m+1] = 1
    else:
        # Neumann
        neumann_count += 1
        if order == 1:
            A[m+1,m:m+2] = np.array([1,-1])/h
        else:
            A[m+1,m-1:m+2] = np.array([-0.5, 2, -1.5])/h
            
    if neumann_count == 2:
        raise ValueError('Zwei Neumann Randbedingungen uebergeben. RWP nicht eindeutig loesbar.')
        
    # Right hand side
    F = f(x)
    F[0] = a
    F[-1] = b
    
    U = np.linalg.solve(A,F)
     
    return U


### MAIN PROGRAM ###
    
# (b) (2 Punkte)
    
btype = {0:1, 1:0}
f = lambda x: np.exp(x)
a = 0
b = 0
ms = [20, 40, 80, 160]
n = len(ms)
E = np.zeros(n)
exa = lambda x: np.exp(x) - x + 1 - np.exp(1)

order = 1
for (m,count) in zip(ms,range(n)):
    U = fd_centered(f,m-1,0,0,btype,order)
    x = np.linspace(0,1,m+1)
    E[count] = np.max(abs(U-exa(x)))
    
EOC = np.log(E[:-1]/E[1:])/np.log(2)
print('Order 1:')
print(EOC)

order = 2
for (m,count) in zip(ms,range(n)):
    U = fd_centered(f,m-1,0,0,btype,order)
    x = np.linspace(0,1,m+1)
    E[count] = np.max(abs(U-exa(x)))
    
EOC = np.log(E[:-1]/E[1:])/np.log(2)
print('Order 2:')
print(EOC)