Differences between revisions 3 and 38 (spanning 35 versions)
Revision 3 as of 2010-02-01 09:48:49
Size: 405
Editor: anonymous
Comment:
Revision 38 as of 2010-02-01 10:58:29
Size: 5850
Comment:
Deletions are marked like this. Additions are marked like this.
Line 7: Line 7:
Vorlesungsmitschrift und Programmcodes von Alexandra Gronholz
Line 10: Line 10:
=== 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
Line 11: Line 118:
=== 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.

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

   1 #bisektion
   2 from pylab import *
   3 
   4 
   5 def F(T):
   6 
   7     y=(1-alpha)*S-k*T**4
   8     return y
   9 
  10 def Abl(T):
  11     ydot=-4*k*T**3
  12     return ydot
  13     
  14 
  15 
  16 close()
  17 alpha=0.3
  18 S=1370./4.      #[W/m^2]
  19 k=5.67e-8       #[W/m^2/K^4]
  20 T1=0.           #[K]
  21 T2=350.         #[K]
  22 t=linspace(0,300,100)
  23 durchlauf=30
  24 a,b=0.,350.
  25 x=b
  26 Texakt=((1.-alpha)*S/k)**(1./4.)
  27 diff=zeros(durchlauf)
  28 diffnew=zeros(durchlauf)
  29 
  30 figure(1)
  31 plot(t,F(t))
  32 plot(t,zeros(t.size),'r')
  33 title('Strahlungsbilanz')
  34 ylabel('Strahlung/Flaeche [W/m^2]')
  35 xlabel('Temperatur')
  36 show()
  37 
  38 
  39 
  40 #nach bisektion
  41 for i in range(durchlauf):
  42     c=(a+b)/2
  43 
  44     if F(a)*F(c)<=0:
  45         b=c
  46     else:
  47         a=c
  48 
  49     print c
  50 
  51     diff[i]=c-Texakt
  52     
  53     #nach newton:
  54     x=x-F(x)/Abl(x)
  55     diffnew[i]=x-Texakt
  56     print diffnew
  57 
  58 
  59 figure(2)
  60 semilogy(abs(diff),'.-',label='Bisektionsverfahren')
  61 plot(abs(diffnew),'g',label='Newton-Raphson-Verfahren')
  62 ylabel('Fehlerwert')
  63 xlabel('Iterationsschritt')
  64 legend()
  65 show()

strahlungbilanz.png bisek_newton.png links: STrahlungbilanz rechts: Vergleich Newton-Raphson-Verfahren - Bisektionsverfahren

Euler-Verfahren

   1 from pylab import *
   2 
   3 L=100 #raeumlich
   4 M=500 #zeitschritte
   5 #c=0.1  #oder unten als sinus def
   6 dx=1./L
   7 dt=1./10
   8 X=linspace(0,1,L)
   9 u=zeros((M,L))
  10 u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
  11 #u[0,45:55]=1. #rechteck
  12 
  13 ion()
  14 figure()
  15 line,=plot(X,u[0,:])
  16 
  17 for n in range(M-1):
  18     
  19     c=0.1*cos((n*dt)/10*pi)
  20     CFL=c*dt/dx
  21     for j in range(L-1):
  22         u[n+1,j]=u[n,j]-CFL/2*(u[n,j+1]-u[n,j-1])
  23                                  
  24     line.set_ydata(u[n+1,:])
  25     draw()

euler_start.png euler_t50.png links: Startsituation rechts: Situation nach 50 Zeitschritten

Finite Differenzen Methode zur Lösung der Advektionsgleichung

Upstream-Verfahren

   1 from pylab import *
   2 
   3 L=100 #raeumlich
   4 M=500 #zeitschritte
   5 #c=0.1  #oder unten als sinus def
   6 dx=1./L
   7 dt=1./10
   8 X=linspace(0,1,L)
   9 u=zeros((M,L))
  10 u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
  11 #u[0,45:55]=1. #rechteck
  12 
  13 ion()
  14 figure()
  15 line,=plot(X,u[0,:])
  16 
  17 for n in range(M-1):
  18     
  19     c=0.1*cos((n*dt)/10*pi)
  20     CFL=c*dt/dx
  21     for j in range(L-1):
  22         if CFL>0:
  23             u[n+1,j]=u[n,j]-CFL*(u[n,j]-u[n,j-1])
  24         else:
  25             u[n+1,j]=u[n,j]-CFL*(u[n,j+1]-u[n,j])
  26                                  
  27     line.set_ydata(u[n+1,:])
  28     draw()

upstream_start.png upstream_t50.png

links: Startsituation rechts: Situation nach 50 Zeitschritten

Upstream-Verfahren mit Werteüberschreibung

   1 from pylab import *
   2 
   3 L=100 #raeumlich           #modifiziert mit ueberschr werten/kein feld mehr
   4 M=500 #zeitschritte
   5 #c=0.1  #oder unten als sinus def
   6 dx=1./L
   7 dt=1./10
   8 X=linspace(0,1,L)
   9 u=zeros((2,L))
  10 u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
  11 #u[0,45:55]=1. #rechteck
  12 
  13 ion()
  14 figure()
  15 line,=plot(X,u[0,:])
  16 
  17 for n in range(M-1):
  18     
  19     c=0.1*cos((n*dt)/10*pi)
  20     CFL=c*dt/dx
  21     for j in range(L-1):
  22         if CFL>0:
  23             u[1,j]=u[0,j]-CFL*(u[0,j]-u[0,j-1])
  24         else:
  25             u[1,j]=u[0,j]-CFL*(u[0,j+1]-u[0,j])
  26                                  
  27     line.set_ydata(u[1,:])
  28     draw()
  29     u[0,:]=u[1,:]

Lax-Wendroff-Verfahren

   1 from pylab import *
   2 
   3 L=100 #raeumlich
   4 M=500 #zeitschritte
   5 #c=0.1  #oder unten als sinus def
   6 dx=1./L
   7 dt=1./10
   8 X=linspace(0,1,L)
   9 u=zeros((M,L))
  10 #u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
  11 u[0,45:55]=1. #rechteck
  12 
  13 ion()
  14 figure()
  15 line,=plot(X,u[0,:])
  16 
  17 for n in range(M-1):
  18     
  19     c=0.1*cos((n*dt)/10*pi)
  20     CFL=c*dt/dx
  21     for j in range(L-1):
  22         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])
  23                                  
  24     line.set_ydata(u[n+1,:])
  25     draw()

laxwendroff_start.png laxwendroff_t50.png links: Startsituation rechts: Situation nach 50 Zeitschritten

Leapfrog-Verfahren

   1 from pylab import *
   2 
   3 L=100 #raeumlich
   4 M=500 #zeitschritte
   5 #c=0.1  #oder unten als sinus def
   6 dx=1./L
   7 dt=1./10
   8 X=linspace(0,1,L)
   9 u=zeros((M,L))
  10 u[0,:]=exp(-500*(X-0.5)**2) #gaussfunktion
  11 #u[0,45:55]=1. #rechteck
  12 
  13 ion()
  14 figure()
  15 line,=plot(X,u[0,:])
  16 
  17 c=0.1*cos((0*dt)/10*pi)
  18 CFL=c*dt/dx
  19    
  20 
  21 for j in range(L-1): #j=raum
  22         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])
  23 
  24 for n in range(1,M-1):
  25     
  26     c=0.1*cos((n*dt)/10*pi)
  27     CFL=c*dt/dx
  28    
  29     
  30     for j in range(L-1):
  31         u[n+1,j]=u[n-1,j]-CFL*(u[n,j+1]-u[n,j-1])
  32                                  
  33     line.set_ydata(u[n+1,:])
  34     draw()

leapfro_start.png 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.

LehreWiki: NumerikProgrammierenEinfuehrung2010 (last edited 2010-02-01 14:47:00 by AlexandraGronholz)