#acl AdminGroup:read,write,delete,revert  EditorGroup:read,write,delete,revert StudentGroup:read,write,delete,revert  All:read
#format wiki
#language en
#pragma section-numbers off

= Einführung in die Grundlagen der Numerik und das Programmieren unter Unix/Linux =
Vorlesungsmitschrift und Programmcodes von Alexandra Gronholz
== Iterative Verfahren zur Bestimmung von Nullstellen ==

=== Newton-Verfahren ===
{{{#!python
#bisektion
from pylab import *


def F(T):

    y=(1-alpha)*S-k*T**4
    return y

def Abl(T):
    ydot=-4*k*T**3
    return ydot
    


close()
alpha=0.3
S=1370./4.	#[W/m^2]
k=5.67e-8	#[W/m^2/K^4]
T1=0.		#[K]
T2=350.		#[K]
t=linspace(0,300,100)
durchlauf=30
a,b=0.,350.
x=b
Texakt=((1.-alpha)*S/k)**(1./4.)
diff=zeros(durchlauf)
diffnew=zeros(durchlauf)

figure(1)
plot(t,F(t))
plot(t,zeros(t.size),'r')
title('Strahlungsbilanz')
ylabel('Strahlung/Flaeche [W/m^2]')
xlabel('Temperatur')
show()



#nach bisektion
for i in range(durchlauf):
    c=(a+b)/2

    if F(a)*F(c)<=0:
        b=c
    else:
        a=c

    print c

    diff[i]=c-Texakt
    
    #nach newton:
    x=x-F(x)/Abl(x)
    diffnew[i]=x-Texakt
    print diffnew


figure(2)
semilogy(abs(diff),'.-',label='Bisektionsverfahren')
plot(abs(diffnew),'g',label='Newton-Raphson-Verfahren')
ylabel('Fehlerwert')
xlabel('Iterationsschritt')
legend()
show()
}}}

{{attachment:strahlungbilanz.png}}
{{attachment:bisek_newton.png}}
'''links:''' STrahlungbilanz
'''rechts:''' Vergleich Newton-Raphson-Verfahren - Bisektionsverfahren

=== Euler-Verfahren ===
{{{#!python
from pylab import *

L=100 #raeumlich
M=500 #zeitschritte
#c=0.1  #oder unten als sinus def
dx=1./L
dt=1./10
X=linspace(0,1,L)
u=zeros((M,L))
u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
#u[0,45:55]=1. #rechteck

ion()
figure()
line,=plot(X,u[0,:])

for n in range(M-1):
    
    c=0.1*cos((n*dt)/10*pi)
    CFL=c*dt/dx
    for j in range(L-1):
        u[n+1,j]=u[n,j]-CFL/2*(u[n,j+1]-u[n,j-1])
                                 
    line.set_ydata(u[n+1,:])
    draw()
}}}
{{attachment:euler_start.png}}
{{attachment:euler_t50.png}}
'''links:''' Startsituation
'''rechts:''' Situation nach 50 Zeitschritten

== Finite Differenzen Methode zur Lösung der Advektionsgleichung ==
=== Upstream-Verfahren ===
{{{#!python 

from pylab import *

L=100 #raeumlich
M=500 #zeitschritte
#c=0.1  #oder unten als sinus def
dx=1./L
dt=1./10
X=linspace(0,1,L)
u=zeros((M,L))
u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
#u[0,45:55]=1. #rechteck

ion()
figure()
line,=plot(X,u[0,:])

for n in range(M-1):
    
    c=0.1*cos((n*dt)/10*pi)
    CFL=c*dt/dx
    for j in range(L-1):
        if CFL>0:
            u[n+1,j]=u[n,j]-CFL*(u[n,j]-u[n,j-1])
        else:
            u[n+1,j]=u[n,j]-CFL*(u[n,j+1]-u[n,j])
                                 
    line.set_ydata(u[n+1,:])
    draw()
}}}
{{attachment:upstream_start.png}}
{{attachment:upstream_t50.png}}

'''links:''' Startsituation
'''rechts:''' Situation nach 50 Zeitschritten
=== Upstream-Verfahren mit Werteüberschreibung ===
{{{#!python

from pylab import *

L=100 #raeumlich           #modifiziert mit ueberschr werten/kein feld mehr
M=500 #zeitschritte
#c=0.1  #oder unten als sinus def
dx=1./L
dt=1./10
X=linspace(0,1,L)
u=zeros((2,L))
u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
#u[0,45:55]=1. #rechteck

ion()
figure()
line,=plot(X,u[0,:])

for n in range(M-1):
    
    c=0.1*cos((n*dt)/10*pi)
    CFL=c*dt/dx
    for j in range(L-1):
        if CFL>0:
            u[1,j]=u[0,j]-CFL*(u[0,j]-u[0,j-1])
        else:
            u[1,j]=u[0,j]-CFL*(u[0,j+1]-u[0,j])
                                 
    line.set_ydata(u[1,:])
    draw()
    u[0,:]=u[1,:]

}}}

=== Lax-Wendroff-Verfahren ===
{{{#!python

from pylab import *

L=100 #raeumlich
M=500 #zeitschritte
#c=0.1  #oder unten als sinus def
dx=1./L
dt=1./10
X=linspace(0,1,L)
u=zeros((M,L))
#u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
u[0,45:55]=1. #rechteck

ion()
figure()
line,=plot(X,u[0,:])

for n in range(M-1):
    
    c=0.1*cos((n*dt)/10*pi)
    CFL=c*dt/dx
    for j in range(L-1):
        u[n+1,j]=u[n,j]-CFL/2*(u[n,j+1]-u[n,j-1])+CFL**2/2*(u[n,j+1]-2*u[n,j]+u[n,j-1])
                                 
    line.set_ydata(u[n+1,:])
    draw()
}}}
{{attachment:laxwendroff_start.png}}
{{attachment:laxwendroff_t50.png}}
'''links:''' Startsituation
'''rechts:''' Situation nach 50 Zeitschritten


=== Leapfrog-Verfahren ===
{{{#!python

from pylab import *

L=100 #raeumlich
M=500 #zeitschritte
#c=0.1  #oder unten als sinus def
dx=1./L
dt=1./10
X=linspace(0,1,L)
u=zeros((M,L))
u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
#u[0,45:55]=1. #rechteck

ion()
figure()
line,=plot(X,u[0,:])

c=0.1*cos((0*dt)/10*pi)
CFL=c*dt/dx
   

for j in range(L-1): #j=raum
        u[1,j]=u[0,j]-CFL/2*(u[0,j+1]-u[0,j-1])+CFL**2/2*(u[0,j+1]-2*u[0,j]+u[0,j-1])

for n in range(1,M-1):
    
    c=0.1*cos((n*dt)/10*pi)
    CFL=c*dt/dx
   
    
    for j in range(L-1):
        u[n+1,j]=u[n-1,j]-CFL*(u[n,j+1]-u[n,j-1])
                                 
    line.set_ydata(u[n+1,:])
    draw()
}}}
{{attachment:leapfro_start.png}}
{{attachment:leapfrog_t50.png}}
'''links:'''  Startsituation
'''rechts:''' Situation nach 50 Zeitschritten

----
'''FAZIT'''

Das '''Eulerverfahren''' zeit starke Instabilität.
Während das '''Upstream-Verfahren''' zwar stabiler ist als das '''Euler-Verfahren'''  (sofern CFL-Bedingung erfüllt), zeigt es dennoch starke Diffusion.
Das '''Lax-Wendroff-Verfahren''' ist weniger diffusiv als das '''Upstrem-Verfahren''' und das '''Leapfrog-Verfahren''' zeigt keine Diffusion mehr, dafür jedoch starke Dispersion.

----

=== Ozeanische Grenzschicht ===
{{{#!python

from pylab import *

h=10.           #Deckschichtdicke in  [m]
Cp=4000.        #Waermekapazitaet Ozean in [J/Kg/K]
S=30.           #bei Salzgehalt von...
d=1000          #Anzahl der Durchlaeufe
T=zeros(d)     #bei Temp von...
T[0]=10        #Startwert Temperatur
dt=60*60*24          #Zeitschritteneinheiten [s]-->Tage mit 60*60*24
start=0        #zeitlicher Startwert
#ti=a=r_[start:d:dt] #array mit Zeitpunkten
alph=0.06      #Albedowert Wasser
epsi=0.97      #Emissivitaet
sig=5.67e-8    #Boltzmannkonstante
I=375          #Einstrahlung in [W/m2]-->375: Gleichgewichtswert
rho=1024


for i in range(d-1):
    T[i+1]=T[i]+((((1-alph)*I-epsi*sig*(T[i]+273)**4)*dt)/(h*Cp*rho))

plot(T)
xlabel('Zeit in Tagen')
ylabel(r'Temperatur in $^\circ$C')
title('Energiebilanz der ozeanischen Grenzschicht')
show()
}}}
{{attachment:energiebilanz.png}}

=== Ozeanische Grenzschicht in Fortran ===
{{{

program ozean

implicit none   !damit variablen nicht vordefiniert

integer, parameter :: d=1000  !Anzahl der Durchlaeufe, + als PARAMETER wenn spaeter nochmal benutzt
real :: T(0:d-1)   !automatisch Nullen
real :: h,Cp,S,dt,alph,epsi,sig,I,rho
integer :: j  ! nicht 2mal I, da keine unterscheidung in GROSS U KLEIN!

h=10.           !Deckschichtdicke in  [m]
Cp=4000.        !Waermekapazitaet Ozean in [J/Kg/K]
S=30.           !bei Salzgehalt von...
T(0)=10.        !Startwert Temperatur
dt=60*60*24.          !Zeitschritteneinheiten [s]-->Tage mit 60*60*24
alph=0.06      !Albedowert Wasser
epsi=0.97      !Emissivitaet
sig=5.67e-8    !Boltzmannkonstante
I=375.          !Einstrahlung in [W/m2]-->375: Gleichgewichtswert
rho=1024.


do j=0,d-1
    T(j+1)=T(j)+((((1-alph)*I-epsi*sig*(T(j)+273)**4)*dt)/(h*Cp*rho))
end do

print*, T(d-1)

end program ozean

 !Informationen zum Erzeugen:
 !a) zum kompilieren 'gfortran + name.f90'
 !erzeugt a.out
 !b) ./a.out oder 
 ! (ganzer weg) ohne './' (finde ueber pwd)
}}}