Sie sind hier: Numerik > Numerisches Lösen von Anfangswertaufgaben

Numerisches Lösen von Anfangswertaufgaben

Samstag 03. Februar 2007 von Simon Praetorius

Übersicht:

  1. Runge-Kutta-Verfahren
Approximation mit 5 StützstellenApproximation mit 5 Stützstellen
Die einfachsten Verfahren zur Lösung von gewöhnlichen Anfangswertproblem, d.h. Differentialgleichungen der Form y'=f(x,y) für x∈[x0,xN] mit einer Anfangsbedingung y(x0)=y0, ersetzen Ableitungsterme durch entsprechende Approximationen, wie Differenzenquotienten, an einer Stelle xi, z.B. y'(xi)≈(y(xi)-y(xi-h))/h, wobei xi=x0+i*h und h=(x0-xN)/N wenn wir h konstant halten.
Damit kann man eine einfache Iteration formulieren, die ausgehend von x0 jedes xi bis hin zu xN berechnet. Die oben genannte Approximation der Ableitung wird mit dem Eulerverfahren umgesetzt. Auch andere Approximationen und Ansätze sind denkbar. Hier werden einige Verfahren verglichen.

Alexander Voigt hat mir einige Fortran-Quelltexte zur Verfügung gestellt und ich habe passende Matlab-Varianten implementiert, um die Kurven auch sichtbar zu machen:

Eulers Polygonzugverfahren:
fortran Code
  • program eulerdiff
  •  
  • use dglmod
  • implicit none
  •  
  • real, dimension(:), allocatable :: t, y
  • real :: a, b, h, y0
  • integer :: n,i
  •  
  • write(*,*) 'y0,a,b,n'
  • read(*,*) y0,a,b,n
  •  
  • allocate(t(1:n))
  • allocate(y(1:n))
  •  
  • ! Initialisierung
  • h=(b-a)/n
  • y(0)=y0
  • t(0)=a
  •  
  • ! einfachste Version der numerischen Lsung einer DGL, sehr ungenau
  • !
  • ! nach Eulers Polygonzug-Verfahren
  • do i=1,n
  • y(i)=y(i-1) + h*f(t(i-1), y(i-1))
  •  
  • t(i)=t(i-1) + h
  • end do
  •  
  • ! Ergebnisausgabe
  • do i=0,n
  • write(*,*) '(',t(i),' , ',y(i),')'
  • end do
  •  
  • end program


modifiziertes Eulerverfahren:
fortran Code
  • program verbesseuler
  •  
  • use dglmod
  • implicit none
  •  
  •  
  • real, dimension(:), allocatable :: t, y
  • real :: a, b, h, y0, k1, k2
  • integer :: n,i
  •  
  • write(*,*) 'y0,a,b,n'
  • read(*,*) y0,a,b,n
  •  
  • allocate(t(1:n))
  • allocate(y(1:n))
  •  
  • h=(b-a)/n
  • t(0)=a
  • y(0)=y0
  •  
  • do i=1,n
  • k1=f(t(i-1), y(i-1))
  • k2=f(t(i-1) + (h/2), y(i-1) + (h/2)*k1)
  • y(i)=y(i-1) + h*k2
  •  
  • t(i)=t(i-1) + h
  • end do
  •  
  • ! Ergebnisausgabe
  • do i=0,n
  • write(*,*) '(',t(i),' , ',y(i),')'
  • end do
  •  
  • end program


Verfahren von Karl Heun:
fortran Code
  • program verbesseuler
  •  
  • use dglmod
  • implicit none
  •  
  •  
  • real, dimension(:), allocatable :: t, y
  • real :: a, b, h, y0, k1, k2, k3
  • integer :: n,i
  •  
  • write(*,*) 'y0,a,b,n'
  • read(*,*) y0,a,b,n
  •  
  • allocate(t(1:n))
  • allocate(y(1:n))
  •  
  • h=(b-a)/n
  • t(0)=a
  • y(0)=y0
  •  
  • do i=1,n
  • t(i)=t(i-1) + h
  •  
  • k1=y(i-1)+h*f(t(i-1), y(i-1))
  • k2=f(t(i),k1)
  • k3=f(t(i-1), y(i-1))
  • y(i)=y(i-1) + (h/2)*(k2+k3)
  •  
  • ! f ist funktion die die jeweilige spezifische f(t,y(t)) enthlt
  • ! Methode ist Trapezregel zur Annherung des Integrals y(t(k)) y(t(k-1)) , da y(t(k)) unbekannt ist, wird es !mit der Eulerregel angenhert berechnet
  • !
  • ! Verfahren hat die Konvergenzordnung 2
  • end do
  •  
  • ! Ergebnisausgabe
  • do i=0,n
  • write(*,*) '(',t(i),' , ',y(i),')'
  • end do
  •  
  • end program


klassisches Runke-Kutta-Verfahren:
fortran Code
  • program runge_kutta
  • use dglmod
  • implicit none
  •  
  •  
  • real, dimension(:), allocatable :: t, y
  • real :: a, b, h, y0, k1, k2, k3, k4
  • integer :: n,i
  •  
  • write(*,*) 'y0,a,b,n'
  • read(*,*) y0,a,b,n
  •  
  • allocate(t(1:n))
  • allocate(y(1:n))
  •  
  • h=(b-a)/n
  • t(0)=a
  • y(0)=y0
  •  
  • do i=1,n
  • k1=f(t(i-1), y(i-1))
  • k2=f(t(i-1) + h/2, y(i-1) + (h/2)*k1)
  • k3=f(t(i-1) + h/2, y(i-1) + (h/2)*k2)
  • k4=f(t(i-1) + h, y(i-1) + h*k3)
  •  
  • y(i)=y(i-1) + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
  •  
  • t(i)=t(i-1) + h
  •  
  • ! f ist funktion die die jeweilige spezifische f(t,y(t)) enthlt
  • ! Runge Kutta aehnelt der Simpsonschen Fassregel, und liefert einen gewichteten Durchschnitt der !Anstiege verschiedener Punkte zwischen y(t(k)) , y(t(k+1))
  • ! Es wird Modul dglmod mit der funktion f (bergebene Argumente (t,y) benoetigt )
  • end
  • end do
  •  
  • ! Ergebnisausgabe
  • do i=0,n
  • write(*,*) '(',t(i),' , ',y(i),')'
  • end do
  •  
  • end program


Un dazu die notwendige Bibliothek, die eine Beispielfunktion enthält:
fortran Code
  • module dglmod
  • use config
  • implicit none
  •  
  • public
  • real(kind=lang),parameter :: tol = 10D0**(-5)
  •  
  • interface operator(.x.)
  • module procedure matrixVectorProd
  • module procedure vectorMatrixProd
  • end interface operator(.x.)
  •  
  • interface operator(.rund.)
  • module procedure equalsApprox
  • end interface operator(.rund.)
  •  
  • interface norm
  • module procedure vectorNorm
  • module procedure matrixNorm
  • end interface norm
  •  
  • contains
  •  
  • function matrixVectorProd(A,b) result(erg)
  • real(kind=lang),dimension(:,:),intent(in) :: A
  • real(kind=lang),dimension(size(A,2)),intent(in) :: b
  • real(kind=lang),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=lang),dimension(:,:),intent(in) :: B
  • real(kind=lang),dimension(size(B,1)),intent(in) :: a
  • real(kind=lang),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
  •  
  • ! Vektornormen
  • function vectorNorm(a,typ)
  • real(kind=lang),dimension(:),intent(in) :: a
  • integer,optional :: typ
  • real(kind=lang) :: vectorNorm
  • if(.not.present(typ)) then; typ=2; endif
  •  
  • if(typ == 0) then
  • ! Maximumsnorm
  • vectorNorm = maxval(abs(a))
  • elseif(typ == 1) then
  • ! Summennorm (l1-Norm)
  • vectorNorm = sum(abs(a))
  • elseif(typ == 2) then
  • ! Euklidische Norm
  • vectorNorm = sqrt(dot_product(a,a))
  • elseif(typ == 3) then
  • ! l3-Norm
  • vectorNorm = sum(abs(a)**3)**(1D0/3)
  • endif
  • end function vectorNorm
  •  
  • ! Matrixnormen
  • function matrixNorm(A,typ)
  • real(kind=lang),dimension(:,:),intent(in) :: A
  • integer,optional :: typ
  • real(kind=lang) :: matrixNorm
  • if(.not.present(typ)) then; typ=1; endif
  •  
  • if(typ == 0) then
  • ! Zeilensumennorm
  • matrixNorm = maxval(sum(abs(A),1))
  • elseif(typ == 1) then
  • ! Spaltensummennorm
  • matrixNorm = maxval(sum(abs(A),2))
  • endif
  • end function matrixNorm
  •  
  • function equalsApprox(a,b)
  • real(kind=lang),intent(in) :: a,b
  • logical :: equalsApprox
  • equalsApprox = (abs(a-b) < tol)
  • end function equalsApprox
  •  
  • end module


Um die Genauigkeit zu untersuchen, habe ich in Matlab die einzelnen Verfahren implementiert und an eine Grafikausgabe "gekoppelt". Allerdings habe ich alle Verfahren in eine Datei gepackt und nur die Beispielfunktion ausgelagert:

matlab Code
  • a = input('Intervall von a: ');
  • b = input('Intervall bis b: ');
  • y0 = input('y(a) als Anfangsbedingung: ');
  • n = input('Anzahl der Teilintervall n: ');
  •  
  • t=linspace(a,b,n);
  • y=zeros(n,4);
  •  
  • y(1,:) = y0;
  •  
  • h = (b-a)/(n-1);
  •  
  • % Euler-Verfahren
  • for i=2:n
  • y(i,1)=y(i-1,1) + h*f(t(i-1), y(i-1,1));
  • end
  •  
  • % modifiziertes Eulerverfahren
  • for i=2:n
  • k1=f(t(i-1), y(i-1,2));
  • k2=f(t(i-1) + (h/2), y(i-1,2) + (h/2)*k1);
  •  
  • y(i,2)=y(i-1,2) + h*k2;
  • end
  •  
  • % Heun-Verfahren
  • for i=2:n
  • k1=y(i-1,3)+h*f(t(i-1), y(i-1,3));
  • k2=f(t(i),k1);
  • k3=f(t(i-1), y(i-1,3));
  •  
  • y(i,3)=y(i-1,3) + (h/2)*(k2+k3);
  • end
  •  
  • % Runge-Kutta-Verfahren
  • for i=2:n
  • k1=f(t(i-1), y(i-1,4));
  • k2=f(t(i-1) + h/2, y(i-1,4) + (h/2)*k1);
  • k3=f(t(i-1) + h/2, y(i-1,4) + (h/2)*k2);
  • k4=f(t(i-1) + h, y(i-1,4) + h*k3);
  •  
  • y(i,4)=y(i-1,4) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
  • end
  •  
  • y
  •  
  • x=linspace(a,b,500);
  • y_erg=f_erg(x,a,y0);
  •  
  • plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),x,y_erg)
  • legend('Euler','Mod-Euler','Heun','RK','Original')

Beliebige Differentialgleichung, z.B. y' = 3xy
matlab Code
  • function erg = f(t,u)
  • erg = 3*t*u;

Die exakte Lösung für ein konkretes Anfangswertproblem mit y(a)=y0 ist folgende:
y = e3/2x2*y0*e-3/2*a2

Und in Matlab-Funktion implementiert so, dass auch vektorwertige Eingaben akzeptiert werden:
matlab Code
  • function erg=f_erg(x,a,y0)
  • erg=exp(1.5 .* x.*x).*(y0*exp(-1.5*a^2));

Approximation mit 10 StützstellenApproximation mit 10 Stützstellen

Wie man vielleicht auf den Bildern erkennen kann, ist das einfache Eulerverfahren sehr ungeeignet für eine gute Approximation und das Runge-Kutta-Verfahren trifft die Funktion schon recht ordentlich. Das einfache Euler hat eine Konsistenzordnung 1 und das Runge-Kutta-Verfahren schon eine Ordnung 4.
Die Runge-Kutta-Verfahren bilden eine allgemeine Klasse von Einschrittverfahren, die in einem anderen Artikel genauer beleuchtet werden.

2 Kommentare

  1. Donnerstag 01.02.2007 12:15 von
    Alexander

    Hallo Simon. Ich hoffe mal jemand kanns gebrauchen. Zum überprüfen eigener Lösungen reicht es auf jedenfall, auch wenns nur diskrete Angaben der Funktion liefert...

  2. Donnerstag 15.02.2007 12:42 von
    Simon

    Mir ging es darum, einige Verfahren zu veranschaulichen. Vielleicht könnte ich auch mal einen Geschwindigkeitsvergleich zwischen Fortran und Matlab machen...

Kommentar hinzufügen

Dieses Feld bitten nicht ausfüllen: