Integration durch Spline-Interpolation

• Man berechnet die kubische Spline-Interpolierende s.
• Als Randbedingungen fordert man, daß s dreimal stetig differenzierbar in und sein soll.
• Die Integration von s erfolgt segmentweise mit Hilfe der Simpsonregel.

program Integration2(input, output);

    {Integration durch Splines}
    
    uses
        stdfile;
        
    const
        MaxArray = 255;
    
    type
        MyArray = array[0..MaxArray] of double;
        
    var
        t, y:    MyArray;
        x:       double;
        n, i:    integer;
        h0, h1:  double;
        a, b:    double;
        u, v, w: MyArray;
        nd, hd:  MyArray;
        c:       MyArray;
        q:       double;
        h, res:  double;
                
    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
            writeln('Integration durch Splines');
            writeln;
            writeln('Welche Datei enthaelt die Punktepaare? (In integr2.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                n := -1;
                repeat
                    n := n + 1;
                    if n > MaxArray then begin
                         writeln('Die Datei "', fName, '" hat zu viele ',
                             'Stuetzstellen (mehr als ', MaxArray, ').');
                         halt
                     end {if}
                 until not StdRead2(f, t[n], y[n]);
                 n := n - 1
            end else begin {if}
                writeln('Kann Datei "', fName, '" nicht oeffnen.');
                halt
            end; {if}
            close(f);
            writeln('Insgesamt ', n + 1, ' Stuetzstellen gelesen.');
            writeln
        end; {init}

    procedure Randbedingungen(var c: MyArray; u, v, w, y, t: MyArray);
        var
            h1, hn, h2, hm:     double;
            d1, dn, d2, dm:     double;
            b1, b2:             double;
            alpha, beta:        double;
            a11, a12, a21, a22: double;
            det:                double;
            i:                  integer;
        begin {Randbedingungen}
            h1 := t[1] - t[0];
            hn := t[n] - t[n - 1];
            d1 := (y[1] - y[0]) / h1;
            dn := (y[n] - y[n - 1]) / hn;
            b1 := 3 * d1 - u[1];
            b2 := 3 * dn - u[n - 1];
            h2 := t[2] - t[1];
            hm := t[n - 1] - t[n - 2];
            d2 := (y[2] - y[1]) / h1;
            dm := (y[n - 1] - y[n - 2]) / hm;
            a11 := (v[1] + 1) / sqr(h1) - (v[1] + v[2]) / sqr(h2);
            a12 := w[1] / sqr(h1) - (w[1] + w[2]) / sqr(h2);
            b1 := (u[1] + u[2] - 2 * d2) / sqr(h2) - (u[1] - 2 * d1) / sqr(h1);
            a21 := (v[n - 1] + v[n - 2]) / sqr(hm) - v[n - 1] / sqr(hn);
            a22 := (w[n - 1] + w[n - 2]) / sqr(hm) - (1 + w[n - 1]) / sqr(hn);
            det := a11 * a22 - a12 * a21;
            b2 := (2 * dm - u[n - 1] - u[n - 2]) / sqr(hm) + (u[n - 1] - 2 * dn) /sqr(hn);
            alpha := (b1 * a22 - a12 * b2) / det;
            beta := (a11 * b2 - b1 * a21) / det;
            c[0] := alpha;
            for i := 1 to n - 1 do begin
                c[i] := u[i] + alpha * v[i] + beta * w[i]
            end; {for}
            c[n] := beta
        end; {Randbedingungen}

    function Spline(x: double; t, y, c: MyArray): double;
        var
            k:              integer;
            h, delta:       double;
            a0, a1, a2, a3: double;
        begin {Spline}
            k := n + 1;
            repeat
                k := k - 1
            until (x > t[k]) or (k = 0);
            k := k + 1;
            h := t[k] - t[k - 1];
            delta := (y[k] - y[k - 1]) / h;
            a0 := y[k];
            a1 := c[k];
            a2 := (c[k] - delta) / h;
            a3 := (c[k] - 2 * delta + c[k - 1]) / sqr(h);
            Spline := ((a3 * (x - t[k - 1]) + a2) * (x - t[k]) + a1) * (x- t[k]) + a0
        end; {Spline}
        
    procedure result;
        begin {result}
            writeln('Das Integral liefert den Wert ', res:8:4)
        end; {result}

    begin {Splines}
        init;
        h0 := t[1] - t[0];
        a := (y[1] - y[0]) / sqr(h0);
        v[1] := -1 / h0;
        for i := 1 to n - 1 do begin
            h1 := t[i + 1] - t[i];
            b := (y[i + 1] - y[i]) / sqr(h1);
            u[i] := 3 * (a + b);
            hd[i] := 2 * (1 / h0 + 1 / h1);
            nd[i] := 1 / h1;
            a := b;
            h0 := h1
        end; {for}
        w[n - 1] := -1 / h1;
        {Vorwaertselimination}
        for i := 1 to n - 2 do begin
            q := -nd[i] / hd[i];
            hd[i + 1] := hd[i + 1] + nd[i] * q;
            u[i + 1] := u[i + 1] + q * u[i];
            v[i + 1] := q * v[i]
        end; {for}
        {Rueckwaertselimination}
        for i := n - 1 to 2 do begin
            q := hd[i];
            u[i] := u[i] / q;
            v[i] := v[i] / q;
            w[i] := w[i] / q;
            q := nd[i - 1];
            u[i - 1] := u[i - 1] - q * u[i];
            v[i - 1] := v[i - 1] - q * v[i];
            w[i - 1] := -q * w[i]
        end; {for}
        q := hd[1];
        u[1] := u[1] / q;
        v[1] := v[1] / q;
        w[1] := w[1] / q;
        Randbedingungen(c, u, v, w, y, t);
        res := 0;
        for i := 1 to n do begin
            h := t[i] - t[i - 1];
            res := res + h / 6 * (y[i - 1] + 4 * Spline(t[i - 1] + h/2, t, y, c) + y[i])
        end; {for}
        result     
    end. {Integration2}