Sie sind hier: Numerik > Numerisches Lösen von Anfangswertaufgaben > Runge-Kutta-Verfahren

Runge-Kutta-Verfahren

Mittwoch 21. November 2007 von
Simon Praetorius
Differentialgleichung von "van der Pol"Differentialgleichung von "van der Pol"
Zu untersuchen ist ein Anfangswertproblem (AWP) der Form y' = f(x,y) mit den Anfangsbedingungen y(x0) = y0 auf einem Intervall [x0,xe]. Runge-Kutta-Verfahren (RKV) bilden eine Klasse von Einschrittverfahren, die versuchen diese AWP numerisch zu lösen. Einschrittverfahren werden sie genannte, weil die Verfahren auf (genau) eine vorher berechnete Stelle zurückgreifen, d.h. will man y(xi+1) näherungsweise berechnen, wird nur xi und y(xi) im Verfahren verwendet. Mehrschrittverfahren nutzen zusätzlich weiter zurückliegende vorher berechnete Stellen in der Verfahrensvorschrift aus.
Die verschiedenen RKV unterscheiden sich unter Anderem in der Anzahl der "Stufen". Diese Stufen haben etwas mit der Anzahl der Funktionsauswertungen der Funktion f pro Schritt zu tun. Die einzelnen Funktionsauswertungen werden in Koeffizienten ki zwischengespeichert und in der Verfahrensfunktion mehrfach wiederverwendet.
Zu der Klasse der RKV gehören z.B. das Eulerverfahren (einstufig) oder Heunverfahren (3-stufig) und es gibt weitere viel höherstufige Verfahren.

Explizite, s-stufige RKV berechnen eine Näherung yi für y(xi) auf einem (eventuell äquidistanten) Gitter xi (i=0,...,N, xN=xe) durch den Ansatz yi+1 = yi + hi⋅Φ(xi,yi,hi) und einer Verfahrensfunktion
Φ(x,y,h) =
s
i=1
bi ki = b*k
mit
ki = f(x + h⋅ci, y + h(
i-1
j=1
aij kj))

Die Koeffizienten ci, bi und aij legen das Verfahren eindeutig fest. Das Besondere an Einschrittverfahren ist, dass bei der Berechnung der ki nur die untere Dreiecksmatrix von A = (aij)i,j=1...s verwendet wird. Bei Mehrschrittverfahren sind auch Einträge in der oberen Dreiecksmatrix von Bedeutung.
A, b und c bestimmen auch die Konsistenzordnung des Verfahrens, d.h. inwieweit das Verfahren das gestellte Problem löst. Damit das explizite RKV mit der Anfangswertaufgabe (AWA) konsistent ist, muss gelten:
s
i=1
bi = 1

Das ist die erste Konsistenzbedingung. Für einige Konsistenzordnung p können zusätzlich weitere Bedingungen an die Koeffizienten aufgestellt werden, z.B. für Konsistenzordnung 2 muss zusätzlich gelten:
s
i=1
bi ci = 1/2

Und für die Konsistenzordnung 3 zusätzlich
s
i=1
bi ci2 = 1/3 und
s
i,j=1
bi aij cj = 1/6

Ein Ziel beim Entwurf von numerischen Verfahren zur Lösung von AWP ist sicher eine hohe Konsistenzordnung (und damit oft auch eine hohe Konvergenzordnung) zu erzeugen. Weitere wichtige Charakteristiken eines "guten" Verfahrens sind z.B. Stabilität (z.B. bei gewissen Testaufgaben) und rechnerischer Aufwand. Bei Runge-Kutta-Verfahren lassen sich Verfahren mit recht hoher Konsistenzordnungen finden, mit höherer Konsistenzordnung steigt allerdings auch der rechnerische Aufwand stark an, so dass man dies nicht allzu weit treibt. Die expliziten RKV haben leider alle nur ein sehr beschränktes Stabilitätsgebiet und sind daher nicht A-stabil. Trotzdem sind sie für das Lösen von AWA sehr wichtig und z.B. in Matlab mit der Funktion ode45 oder ode23 implementiert und gehören zu den Standardlösern für einfache DGLn. Einige dieser RKV möchte ich hier konkret vorstellen:
Die Koeffizienten werden der Übersicht halber auch oft in das sogenannte Butcher-Schema geschrieben:
c A
bT

z.B. für das Euler-Verfahren:
0  
1

oder das Heun-Verfahren:
0      
1/3 1/3    
2/3 0 1/3  
1/4 0 3/4

Weitere (größere) Tableaus können z.B. auf der englischen Wikipediaseite unter Runge-Kutta-Verfahren gefunden werden oder auf einigen anderen Internetseiten oder in diverser Literatur (z.B. Hairer, Wanner: "Solving ordinary differential equations I").

Will man dieses Verfahren nun implementieren, genügt eine Schleife mit N Duchläufen, in der die einzelnen yi der Reihe nach berechnet werden, d.h. in jedem Schleifendurchlauf werden zuerst die ki dann die Φi und zuletzt die yi berechnet. Der Einfachheit halber werden zuerst äquidistante Gitter (d.h. h = 1/N ist konstant) verwendet. Wir schreiben also statt hi einfach nur h. Es ist weiterhin darauf zu achten, dass die Funktion f ein mehrdimensionales Ergebnis zurückliefern kann, wenn man z.B. eine Differentialgleichung 2ter Ordnung betrachtet (wie in der Beispielimplementierung). Dadurch müssen in einigen Termen Matrix-Vektor- bzw. Vektor-Matrixprodukte gebildet werden anstelle von Skalarprodukten.

Zur Implementierung des Verfahrens in Fortran habe ich einige Module geschrieben, die z.B. das Ergebnis des Verfahrens als verkettete Liste zurückliefern. Diese Module werden in einem anderen Artikel beschrieben: Collections
fortran Code
  • module rkv
  • use contmod, only: length, setContent
  • use collections
  • implicit none
  •  
  • interface operator(.x.)
  • module procedure matrixVectorProd
  • module procedure vectorMatrixProd
  • end interface operator(.x.)
  •  
  • contains
  •  
  • subroutine solver(f, erg, x0, xe, y0, h, A, b, c)
  • interface
  • function f(x, y) result(z)
  • use contmod, only: length
  • real(kind=8), intent(in) :: x
  • real(kind=8), dimension(length), intent(in) :: y
  • real(kind=8), dimension(length) :: z
  • end function f
  • end interface
  •  
  • ! Ergebnis als Liste zur?eben
  • type(list), intent(inout) :: erg
  • ! zu berechnendes Intervall
  • real(kind=8), intent(in) :: x0, xe
  • ! Startwerte
  • real(kind=8), dimension(length), intent(in) :: y0
  • ! Schrittweite
  • real(kind=8) :: h
  • ! Parameter des Verfahrens
  • real(kind=8), dimension(:,:) :: A
  • real(kind=8), dimension(size(A,1)) :: b, c
  •  
  • real(kind=8) :: x
  • real(kind=8), dimension(length) :: y
  • real(kind=8), dimension(size(A,1),length) :: k
  • integer :: i,s
  •  
  • if(size(A,1)/=size(A,2)) then
  • write(*,*) 'Fehler (in Subroutine solver): uebergebene Matrix der Verfahrenskeffizienten hat die falsche Dimension.'
  • return
  • endif
  •  
  • s = size(A,1)
  •  
  • x = x0
  • y = y0
  • call add(erg,setContent(x,y))
  •  
  • do
  • if(x>xe) then
  • exit
  • endif
  •  
  • do i=1,s
  • k(i,:) = f(x + c(i)*h, y + h*(A(i,1:i-1) .x. k(1:i-1,:)))
  • end do
  •  
  • x = x + h
  • y = y + h*(b .x. k)
  •  
  • call add(erg, setContent(x,y))
  • end do
  •  
  • end subroutine solver
  •  
  • function matrixVectorProd(A,b) result(erg)
  • real(kind=8),dimension(:,:),intent(in) :: A
  • real(kind=8),dimension(size(A,2)),intent(in) :: b
  • real(kind=8),dimension(size(A,1)) :: erg
  • integer :: i
  • do i=1,size(A,1)
  • erg(i) = dot_product(A(i,:),b)
  • end do
  • end function matrixVectorProd
  •  
  • function vectorMatrixProd(a,B) result(erg)
  • real(kind=8),dimension(:,:),intent(in) :: B
  • real(kind=8),dimension(size(B,1)),intent(in) :: a
  • real(kind=8),dimension(size(B,2)) :: erg
  • integer :: i
  • do i=1,size(B,2)
  • erg(i) = dot_product(a,B(:,i))
  • end do
  • end function vectorMatrixProd
  •  
  • end module

Zum Testen des Verfahrens habe ich noch ein kleines Hauptprogramm geschrieben:
fortran Code
  • module testmod
  • use contmod, only: length
  •  
  • implicit none
  •  
  • contains
  •  
  • function vdpl(x,y) result(f) ! Dimension 2, Intervall [0,20], Startwerte [2 0]
  • real(kind=8),intent(in) :: x
  • real(kind=8),dimension(length),intent(in) :: y
  • real(kind=8),dimension(length) :: f
  • f(1) = y(2)
  • f(2) = (1-y(1)**2)*y(2)-y(1)
  • end function
  •  
  • end module testmod
  •  
  • program runge_kutta
  • use contmod, only: length
  • use rkv
  • use testmod
  • implicit none
  •  
  • type(list),target :: erg
  • real(8), dimension(:,:), allocatable :: A
  • real(8), dimension(:), allocatable :: b, c
  • real(8) :: x0, xe, h
  • real(8), dimension(length) :: y0
  • call init(erg)
  •  
  • write(*,*) 'Bereich (x0, xe]: '
  • read(*,*) x0, xe
  • write(*,*) 'y(x0) = y0: '
  • read(*,*) y0
  • write(*,*) 'Schrittweite h: '
  • read(*,*) h
  •  
  • ! Benutze zum Testen das Heun-Verfahren
  • allocate(A(3,3), b(3), c(3))
  •  
  • A(:,1) = [0.D0, 1.D0/3, 0.D0]
  • A(:,2) = [0.D0, 0.D0, 1.D0/3]
  • A(:,3) = [0.D0, 0.D0, 0.D0]
  •  
  • b(:) = [0.25D0, 0.D0, 0.75D0]
  • c(:) = [0.D0, 1.D0/3, 2.D0/3]
  •  
  • call solver(vdpl, erg, x0, xe, y0, h, A, b, c)
  • call write(erg)
  •  
  • call clear(erg)
  •  
  • end program runge_kutta

Diese Implementierung arbeitet natürlich noch nicht sehr effizient, aber funktioniert und kann viele Anfangswertaufgaben auch für Funktionen beliebiger Dimension (einigermaßen) gut lösen. Eine bessere Implementierung erreicht man z.B. durch variable Schrittbreiten, d.h. wenn sich die Funktion stark ändert, dann wird eine kleine Schrittweite gewählt, sonst eine größere. Dazu werde ich demnächst einen weiteren Artikel inklusive Beispielimplementierung veröffentlichen.

Kommentar hinzufügen

Dieses Feld bitten nicht ausfüllen: