Differentialgleichungen elegant mit Haskell lösen
Demonstration am Beispiel der Besselgleichung

Die Besselfunktionen J_n und Y_n sind die Basislösungen
für die Besselgleichung
  x^2 * y''(x) + x * y'(x) + (x^2 - m^2) * y(x) = 0  .

Diese Gleichung wollen wir numerisch lösen.

> module Numerics.ODEEuler where
> import Data.List (transpose, zipWith4, )
> import Useful ((^!), )

Fangen wir mit einem einfachen numerischen Integral an.

> integrateStep1 :: Num a => a -> [a] -> [a]
> integrateStep1 = scanl (+)

scanl akkumuliert die Werte einer Folge mit Hilfe einer anderen Funktion, hier (+).
Mit der Addition ergibt das die Berechnung der Partialsummen.

Damit können wir schon die allgemeine explizite Euler-Methode zusammenbauen.
Der Einfachheit halber nehme ich zunächst 1 als Schrittweite.

Haben wir eine Gleichung erster Ordnung, also der Form
  y'(x) = f(x, y(x))
dann sieht das explizite Euler-Verfahren wie folgt aus.

> eulerForward1Step1 :: Num a => (a -> a -> a) -> a -> a -> [a]
> eulerForward1Step1 f x0 y0 =
>    let x   = iterate (1+) x0
>        y   = integrateStep1 y0  y'
>        y'  = zipWith f x y
>    in  y

Wollen wir das Ergebnis statt als Wertefolge als Potenzreihe haben,
müssen wir lediglich die integrateStep1-Funktion ersetzen durch:

> integratePowerSeries :: (Fractional a) => a -> [a] -> [a]
> integratePowerSeries c x = c : zipWith (/) x (map fromInteger [1 ..])

Die neue Koeffizientenfolge beginnt also mit c
und es folgen die Glieder von x dividiert durch 1, 2, 3, ... usw.
Addition und Skalierung von Folgen sind bei Potenzreihen genau so wie bei Wertefolgen.


Zurück zur Lösung der Differentialgleichung als Wertefolge:
Wir wollen eine Liste h der Schrittweiten vorgeben.

> integrate :: Num a => [a] -> a -> [a] -> [a]
> integrate h c y = integrateStep1 c (zipWith (*) h y)

Nun sieht das Euler-Verfahren erster Ordnung mit allgemeiner Verfahrensfunktion phi so aus:

> euler1 :: Num a => (a -> a -> a -> a) -> (a,[a]) -> a -> [a]
> euler1 phi (x0,h) y0 =
>    let x   = scanl (+) x0 h
>        y   = integrate h y0 y'
>        y'  = zipWith3 phi x h y
>    in  y

Beim expliziten Verfahren ist das Vektorfeld f
gleichzeitig die Verfahrensfunktion.
Lediglich die Schrittweite muss ignoriert werden (const . f).

> eulerForward1 :: Num a => (a -> a -> a) -> (a,[a]) -> a -> [a]
> eulerForward1 f = euler1 (const . f)

Zweite Ordnung:
  y''(x) = f(x, y(x), y'(x))

> euler2 :: Num a => (a -> a -> a -> a -> a) -> (a,[a]) -> a -> a -> [a]
> euler2 phi (x0,h) y0 y'0 =
>    let x   = scanl (+) x0 h
>        y   = integrate h y0  (tail y')
>        y'  = integrate h y'0 y''
>        y'' = zipWith4 phi x h y y'
>    in  y

> eulerForward2 :: Num a => (a -> a -> a -> a) -> (a,[a]) -> a -> a -> [a]
> eulerForward2 f = euler2 (const . f)

'tail' entfernt das erste Element einer Liste.
Ohne 'tail' bekämen wir eine unnötige Verzögerung.
Dann würden die ersten beiden Werte der Ergebnisfolge
ausschließlich von den Anfangswerten y0 und y'0 abhängen,
nicht aber von f.

Nun schreckt uns auch eine Gleichung beliebiger Ordnung nicht mehr.
  (D n y)(x) = f(x, (y(x),y'(x),..., (D (n-1) y)(x)))
Die Ordnung wird durch die Anzahl der Anfangswerte festgelegt.

> euler :: Num a => (a -> a -> [a] -> a) -> (a,[a]) -> [a] -> [a]
> euler phi (x0,h) y0s =
>    let integrateTail c y = integrate h c (tail y)
>        x   = scanl (+) x0 h
>        ys  = init (scanr integrateTail (0:yn) y0s)
>        yn  = zipWith3 phi x h (transpose ys)
>    in  head ys

> eulerForward :: Num a => (a -> [a] -> a) -> (a,[a]) -> [a] -> [a]
> eulerForward f = euler (const . f)

'scanr' wird hier nicht für das Akkumulieren von Zahlen verwendet,
sondern für die wiederholte Integration.
'integrateTail' integriert, überspringt aber den ersten Wert einer Folge
um wiederholte Verzögerungen zu vermeiden.


Jetzt können wir die Besselgleichungen
mit dem allgemeinen Euler-Verfahren lösen.

> besselField :: Fractional a => a -> a -> a -> a -> a
> besselField m x y y' = - (y'/x + (1-(m/x)^!2)*y)

> besselForward :: Fractional a => a -> (a,[a]) -> a -> a -> [a]
> besselForward m = eulerForward2 (besselField m)



Nun wollen wir uns noch an einem impliziten Euler-Verfahren versuchen.
Dafür brauchen wir zuerst das Newton-Verfahren.
Das Newton-Verfahren benötigt
  Funktion f, die Funktionswert und Ableitung an einer Stelle berechnet
  Startnäherung x0
Es gibt die (unendliche) Folge aller Näherungswerte zurück.

> newton :: (Fractional a) => (a -> (a,a)) -> a -> [a]
> newton f = iterate (newtonStep f)
> newtonStep :: (Fractional a) => (a -> (a,a)) -> a -> a
> newtonStep f x = let (y,y') = f x in x - y / y'

Als nächstes definieren wir ein paar numerische Grenzwertbildungen.
Mit diesen können wir aus der Liste der Zwischennäherungen vom Newton-Verfahren
die bevorzugte Näherung herausfischen.
Wir können also Näherungsverfahren und Abbruchkriterium trennen
und frei miteinander kombinieren.
Darüberhinaus können wir auch beliebig Konvergenzbeschleuniger dazwischenschalten.

Die einfachste numerische Grenzwertbildung ist sicherlich,
wenn wir nach einer festen Zahl von Schritten abbrechen.

> limFix :: Int -> [a] -> a
> limFix = flip (!!)

Wir können aber auch die Folge abbrechen,
wenn sich aufeinanderfolgende Glieder nur noch wenig unterscheiden.

> limDiff :: (Ord a, Num a) => a -> [a] -> a
> limDiff tol (x0:x1:xs) =
>    if abs (x1-x0) <= tol
>    then x1
>    else limDiff tol (x1:xs)
> limDiff _ _ = error "Endliche Folge, aber kein Glied gen\252gt Abbruchkriterium."

Oder wir brechen die Folge ab, sobald wir nahe genug an einer Nullstelle sind.
Dazu filtern wir aus der Folge alle die Werte heraus,
deren Bilder bzgl. f klein genug sind.
Von all diesen Werten nehmen wir den ersten (head).

> limZero :: (Ord a, Num a) => a -> (a -> a) -> [a] -> a
> limZero tol f = head . filter ((tol >=) . abs . f)

Das wollen wir nun für das implizite Euler-Verfahren verwenden.
Die Integration haben wir schon herausgelöst
und brauchen uns nur noch mit dem Teil des Verfahrens zu beschäftigen,
der spezifisch für die implizite Variante ist,
also mit der Verfahrensfunktion.

Eingabe:
   Vektorfeld f und Ableitung
   Startwert x0, Schrittweiten h, Anfangswert y0

Ausgabe:
   unendliche Liste der Näherungswerte

> backwardSlope1 :: (Fractional a) =>
>    (a -> a -> (a, a)) -> a -> a -> a -> a
> backwardSlope1 f xi hi yi =
>    limFix 20 (newton
>         (\ d ->
>               let (nd, ndy) = f (xi+hi) (yi+hi*d)
>               in  (nd - d, ndy * hi - 1))
>         (fst (f xi yi)))
>
> eulerBackward1 :: Fractional a =>
>    (a -> a -> (a,a)) -> (a,[a]) -> a -> [a]
> eulerBackward1 = euler1 . backwardSlope1

Hier das implizite Euler-Verfahren für Gleichungen zweiter Ordnung -
zumindest das, was ich dafür halte.

> backwardSlope2 :: (Fractional a) =>
>    (a -> a -> a -> (a, (a, a))) -> a -> a -> a -> a -> a
> backwardSlope2 f xi hi yi y'i =
>    limFix 20 (newton
>         (\ d2 ->
>               let d1 = y'i + hi*d2
>                   d0 = yi  + hi*d1
>                   (nd2, (nd2y, nd2y')) = f (xi+hi) d0 d1
>               in  (nd2 - d2, (nd2y * hi + nd2y') * hi - 1))
>         (fst (f xi yi y'i)))

> eulerBackward2 :: Fractional a =>
>    (a -> a -> a -> (a,(a,a))) -> (a,[a]) -> a -> a -> [a]
> eulerBackward2 = euler2 . backwardSlope2

Und jetzt für beliebige Ordnung.

> backwardSlope :: (Fractional a) =>
>    (a -> [a] -> (a, [a])) -> a -> a -> [a] -> a
> backwardSlope f xi hi yi =
>    limFix 20 (newton
>         (\ dn ->
>               let d = init (scanr (\yji dj -> yji + hi * dj) dn yi)
>                   (ndn, ndny) = f (xi+hi) d
>               in  (ndn - dn, foldl1 (\ndnyj dy -> hi * ndnyj + dy) (ndny ++ [-1])))
>         (fst (f xi yi)))
>
> eulerBackward :: Fractional a =>
>    (a -> [a] -> (a,[a])) -> (a,[a]) -> [a] -> [a]
> eulerBackward = euler . backwardSlope


Berechnet Vektorfeld für Besseldifferentialgleichung und seine Ableitungen.

> besselFieldD :: Fractional a => a -> a -> a -> a -> (a,(a,a))
> besselFieldD m x y y' =
>    let c0 = 1-(m/x)^!2
>        c1 = 1/x
>    in  (-(c0*y+c1*y'), (-c0, -c1))

> besselBackward :: Fractional a => a -> (a,[a]) -> a -> a -> [a]
> besselBackward m = eulerBackward2 (besselFieldD m)

> testBessel :: [[(Double, Double)]]
> testBessel =
>    let xh = (1, repeat 0.1)
>        x  = uncurry (scanl (+)) xh
>        y  = besselForward  1 xh 1 2
>        y' = besselBackward 1 xh 1 2
>    in  [zip x y, zip x y']

Allgemeine Testfunktionen

> test :: Bool
> test = and tests

> tests :: [Bool]
> tests =
>    map (all (uncurry (approx 1e-13)) . take 20) testSquareGen

> approx :: Double -> Double -> Double -> Bool
> approx eps x y = abs (x-y) <= abs (x+y) / 2 * eps

> testSolvers :: (Num a) =>
>    (a -> b) -> ((a, [a]) -> b -> [b]) -> ((a, [a]) -> b -> [b]) -> (a, [a]) -> [[(a, b)]]
> testSolvers yFunc yForwSol yBackSol (x0,h) =
>    let x     = scanl (+) x0 h
>        y     = map yFunc x
>        yForw = yForwSol (x0,h) (head y)
>        yBack = yBackSol (x0,h) (head y)
>    in  [zip x y, zip x yForw, zip x yBack]

> xhDefault :: Fractional a => (a, [a])
> xhDefault = 
>       (1, cycle [0.13,0.09,0.01])
> --      (1, repeat 0.1)

> testLinear :: [[(Double,Double)]]
> testLinear =
>    testSolvers
>       (\x -> x/2+1)
>       (eulerForward1  (\x y ->  (y-1)/x))
>       (eulerBackward1 (\x y -> ((y-1)/x, 1/x)))
>       xhDefault

> testSquare1 :: [[(Double,Double)]]
> testSquare1 =
>    testSolvers
>       (^!2)
>       (eulerForward1  (\x y ->  2*y/x))
>       (eulerBackward1 (\x y -> (2*y/x, 2/x)))
>       xhDefault

> testSquare2 :: [[(Double,Double)]]
> testSquare2 =
>    testSolvers
>       (\x -> x^!2)
>       (\xh y0 -> eulerForward2  (\x y y' ->  x*y'/y) xh y0 2)
>       (\xh y0 -> eulerBackward2 (\x y y' -> (x*y'/y, (-x*y'/y^!2, x/y))) xh y0 2)
>       xhDefault

> testSquare3 :: [[(Double,Double)]]
> testSquare3 =
>    testSolvers
>       (\x -> x^!3)
>       (\xh y0 -> eulerForward  (\_ [y, y', y''] ->  y'*y''/(3*y)) xh [y0, 3, 6])
>       (\xh y0 -> eulerBackward (\_ [y, y', y''] -> (y'*y''/(3*y), [y'*y''/(3*y^!2), y''/(3*y), y'/(3*y)])) xh [y0, 3, 6])
>       xhDefault

> testSquareGen :: [[(Double, Double)]]
> testSquareGen =
>    let xh = (1, cycle [0.1,0.2,0.05])
>    in  [zip (eulerForward1 (\x  y  -> 2*y/x) xh  1 )
>             (eulerForward  (\x [y] -> 2*y/x) xh [1]),
>         zip (eulerBackward1 (\x  y  -> (2*y/x,  2/x )) xh  1 )
>             (eulerBackward  (\x [y] -> (2*y/x, [2/x])) xh [1]),
>         zip (eulerForward2 (\x  y  y'  -> x*y'/y) xh  1 2 )
>             (eulerForward  (\x [y, y'] -> x*y'/y) xh [1,2]),
>         zip (eulerBackward2 (\x  y  y'  -> (x*y'/y, (-x*y'/y^!2, x/y))) xh  1  2 )
>             (eulerBackward  (\x [y, y'] -> (x*y'/y, [-x*y'/y^!2, x/y])) xh [1, 2])
>        ]