#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 11 10:06:16 2019

@author: kerkmann
"""

import numpy as np
from sympy import *
import matplotlib.pyplot as plt


x = Symbol('x')
f = x*sin(x)**2

# (a) (1 Punkt)

f1 = diff(f,x)
f2 = diff(f1,x)
f3 = diff(f2,x)

# (b) (2 Punkte)

# erste Ableitung
D11 = lambda x, h: (F(x+h)-F(x))/h
D12 = lambda x, h: (F(x+h)-F(x-h))/2/h
D13 = lambda x, h: (2*F(x+h) + 3*F(x) - 6*F(x-h) + F(x-2*h))/6/h

# zweite Ableitung
D21 = lambda x, h: (2*F(x-h) - 3*F(x) + F(x+2*h))/3/h**2
D22 = lambda x, h: (F(x+h)-2*F(x)+F(x-h))/h**2
D23 = lambda x, h: (11*F(x-h) - 20*F(x) + 6*F(x+h) + 4*F(x+2*h) - F(x+3*h))/12/h**2

# dritte Ableitung
D31 = lambda x, h: (-F(x-h)+3*F(x)-3*F(x+h)+F(x+2*h))/h**3
D32 = lambda x, h: (-F(x-2*h)+2*F(x-h)-2*F(x+h)+F(x+2*h))/2/h**3
D33 = lambda x, h: (-F(x-2*h) - F(x-h) + 10*F(x) - 14*F(x+h) + 7*F(x+2*h) - F(x+3*h))/4/h**3

# (c) (2 Punkte) und (d) (1 Punkt)

F = lambdify(x,f)
F1 = lambdify(x,f1)
F2 = lambdify(x,f2)
F3 = lambdify(x,f3)

def err(d,f):
    return abs(d-f)
        
xbar = 1

hn = np.logspace(-9,-1)
plt.figure(0)
plt.loglog(hn,err(D11(xbar,hn),F1(xbar)),label = 'D11')
plt.loglog(hn,err(D12(xbar,hn),F1(xbar)),label = 'D12')
plt.loglog(hn,err(D13(xbar,hn),F1(xbar)),label = 'D13')
plt.legend()

hn = np.logspace(-6,-1)
plt.figure(1)
plt.loglog(hn,err(D21(xbar,hn),F2(xbar)),label = 'D21')
plt.loglog(hn,err(D22(xbar,hn),F2(xbar)),label = 'D22')
plt.loglog(hn,err(D23(xbar,hn),F2(xbar)),label = 'D23')
plt.legend()

hn = np.logspace(-5,-1)
plt.figure(2)
plt.loglog(hn,err(D31(xbar,hn),F3(xbar)),label = 'D31')
plt.loglog(hn,err(D32(xbar,hn),F3(xbar)),label = 'D32')
plt.loglog(hn,err(D33(xbar,hn),F3(xbar)),label = 'D33')
plt.legend()

# Die Ordnung lässt sich an der Steigung des Graphen ablesen bei Verwendung des log-log Plots (für EOC siehe Übung).