801_SteifeDGL.py

"""
Steife Differentialgleichungen
Beispiel 2 aus der Zusammenfassung:
- Konvergenz fuer N = 999
- Oszillation mit konstantem Betrag fuer N = 1000
- Divergenz fuer N = 1001
@author: Kevin Kazuki Huguenin-Dumittan
"""

#Probiere verschiedene Werte fuer N aus!

from numpy import *
from matplotlib.pyplot import *
from scipy.linalg import *
from scipy.optimize import fsolve


# Explizites Eulerverfahren
def eE(f,t0,tEnd,y0,N):
    dim = size(y0)
    y = zeros((N+1,dim))
    t, h= linspace(t0,tEnd,N+1, retstep=True)
    y[0,:] = y0

    for i in range(N):
        y[i+1,:] = y[i,:] + h*f(t[i],y[i,:])
    return t,y

# Implizites Eulerverfahren
def iE(f,t0,tEnd,y0,N):
    dim = size(y0)
    y = zeros((N+1,dim))
    t, h= linspace(t0,tEnd,N+1, retstep=True)
    y[0,:] = y0

    for i in range(N):
        F = lambda x: x - y[i,:] - h*f(t[i+1],x)
        y[i+1,:] = fsolve(F,y[i,:] + h*f(t[i],y[i,:]))
    return t,y


# Beispiel 2 aus der Zusammenfassung:
# (Anderes gutes Beispiel: Gedaempftes Pendel aus den Runge-Kutta Programmen)
f = lambda t,y: array([998*y[0] + 1998*y[1], -999*y[0]-1999*y[1]])
u0 = 1
v0 = 0
y0 = array([u0,v0])
t0 = 0.
tEnd = 2.


N = 999 # Divergenz, da (tEnd-t0)/N > 0.002
N = 1000 # Oszillation, da (tEnd-t0)/N = 0.002
N = 1001 # Konvergenz, da (tEnd-t0)/N < 0.002


t_eE,y_eE = eE(f,t0,tEnd,y0,N)
t_iE,y_iE = iE(f,t0,tEnd,y0,N)

#Plotten
sol = lambda t:array([-exp(-1000*t)+2*exp(-t),exp(-1000*t)-exp(-t)])
tsol = linspace(t0,tEnd,1000)
y_exact = sol(tsol)

figure()
plot(t_eE,y_eE[:,0],'r.-',lw=0.05, label='u(t) eE')
plot(t_iE,y_iE[:,0],'b.-',label='u(t) iE')
plot(tsol,y_exact[0],'k',label='u(t) Exact')
#plot(t_eE,y_eE[:,1],'m.-',lw=0.05, label='v(t) eE')
#plot(t_iE,y_iE[:,1],'c.-',label='v(t) iE')
#plot(tsol,y_exact[1],'k--',label='v(t) Exact')
legend(loc='lower right')
title('Verhalten von steifen DGL')
xlabel('Zeit t')
ylabel('Position y(t)')
xlim(0,1.5)
ylim(-4,5)
grid(True)
show()


# Plotten: Detaillierte Ansicht des letzten Teiles durch Einschraenkung
# beider Achsen mit xlim und ylim, der Rest ist genau gleich
figure()
plot(t_eE,y_eE[:,0],'ro-',lw=0.25, label='u(t) eE')
plot(t_iE,y_iE[:,0],'b.-',label='u(t) iE')
plot(tsol,y_exact[0],'k',label='u(t) Exact')
legend(loc='lower right')
title('Verhalten von steifen DGL: Detaillierte Ansicht')
xlabel('Zeit t')
ylabel('Position y(t)')
xlim(1.9,2)
#ylim(-0.2,0.8)
ylim(-1,1.7)
grid(True)
show()