#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 24 15:43:04 2019

@author: kerkmann
"""

import numpy as np
import matplotlib.pyplot as plt

# (a) (2 Punkte)

eps = 1/10
e = np.ones(10)
A = 1/eps*(np.diag(e) - np.diag(e[0:-1],-1))
A[0,-1] = -1/eps

ev = np.linalg.eigvals(-A)
plt.figure(1)
plt.plot(np.real(ev),np.imag(ev),'*')

# Explizite Verfahren eignen sich gut fuer dieses Problem. Der Zeitschritt
# k kann von der Groessenordnung O(eps) gewaehlt werden, da der
# betragsgroesste Eigenwert linear mit epsilon waechst (vergleiche dazu
# mehrere Plots). Z.B. expl. Euler mit k <= eps.

# (b) (2 Punkte)

A = 1/2/eps*(np.diag(e[0:-1],1) - np.diag(e[0:-1],-1))
A[0,-1] = -1/2/eps
A[-1,0] = 1/2/eps

ev = np.linalg.eigvals(-A)
plt.figure(2)
plt.plot(np.real(ev),np.imag(ev),'*')

# Die Eigenwerte liegen alle auf der imaginaeren Achse. Das explizite
# Euler Verfahren laesst sich also nie anwenden, unabhaengig von k.
# Andere explizite Verfahren mit einem etwas groesseren Stabilitats-
# gebiet, das insbesondere einen Teil der imaginaeren Achse beinhaltet,
# koennen mit Zeitschritt k von der Groessenordnugn O(eps) gewaehlt werden.

# (c) (2 Punkte)

A = 1/eps**2*(2*np.diag(e) - np.diag(e[0:-1],1) - np.diag(e[0:-1],-1))
A[0,-1] = -1/eps**2
A[-1,0] = -1/eps**2

ev = np.linalg.eigvals(-A)
plt.figure(3)
plt.plot(np.real(ev),np.imag(ev),'*')

# Die Eigenwerte liegen alle auf der reelen Achse. Der betragsgroesste
# Eigenwert waechst wie eps**2, also muss der Zeitschritt eines expliziten
# Verfahrens von der Ordnung k = O(eps**2) gewaehlt werden. Dies ist fuer
# kleine eps sehr unguenstig, daher sollte hier ein implizites Verfahren
# gewaehlt werden, das A(alpha) stabil ist.