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:

• 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 Lösung 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:

• 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:

• 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)) enthält
• ! Methode ist Trapezregel zur Annäherung des Integrals y(t(k)) y(t(k-1)) , da y(t(k)) unbekannt ist, wird es !mit der Eulerregel angenähert 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:

• 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)) enthält
• ! 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:

• 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:

• 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

• 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:

• function erg=f_erg(x,a,y0)
• erg=exp(1.5 .* x.*x).*(y0*exp(-1.5*a^2));

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.