H08 – Polynominterpolation, verschiedene Stützstellen

 

Gesucht ist ein Polynom p, das eine Funktion f in den paarweise verschiedenen Punktepaaren \left( {{x_i},f\left( {{x_i}} \right)} \right) interpoliert. Zur Bestimmung der Koe?zienten sollen die dividierten Differenzen verwendet werden. Die anschließende Auswertung soll mit dem Hornerschema erfolgen.

a )

Schreiben Sie ein m-File (’divdiff.m’), in dem die Koe?zienten des Interpolationspolynoms
mittels dividierter Differenzen bestimmt werden. Dabei sollen die Stützstellen x = \left( {{x_i}} \right)
und Funktionswerte y = \left( {f\left( {{x_i}} \right)} \right) vektoriell übergeben werden. Verwenden Sie dazu folgenden Programmkopf:

function a=divdiff(x,y)
%x Stuetzstellen
%y Funktionswerte an den Stuetzstellen
%a Polynomkoeffizienten (dividierte Differenzen)

b )

Schreiben Sie ein m-File (’horner.m’), in dem das Interpolationspolynom mit dem Hornerschema an der Stelle \xi ausgewertet wird. Verwenden Sie dabei die Terme \left( {x-{x_i}} \right) als Faktoren. Es soll mit folgenden Zeilen beginnen:

function f=horner(x,a,xi)
%x Stuetzstellen des Polynoms
%a Polynomkoeffizienten (dividierte Differenzen)
% xi Auswertungsstelle
%f Ergebnis der Auswertung

c )

Schreiben Sie ein m-File (’u81.m’), in dem die Funktion

f\left( x \right) = \frac{1}{{1+25{x^2}}}

im Intervall \left[ {-1,1} \right] interpoliert wird. Bestimmen Sie die Koe?zienten für das Polynom mit

1. den Stützstellen {x_i} = -1+i \cdot 0,2\quad \quad i = 0, \ldots ,10
2. den Stützstellen

-\frac{6}{5}+\frac{1}{5},\:-\frac{6}{5}+\frac{1}{4},\:-\frac{6}{5}+\frac{1}{3},\:-\frac{6}{5}+\frac{1}{2},-\frac{6}{5}+1,\:\:0,\:\:\frac{6}{5}-1,\frac{6}{5}-\frac{1}{2},\frac{6}{5}-\frac{1}{3},\frac{6}{5}-\frac{1}{4},\:\frac{6}{5}-\frac{1}{5}

3. den Nullstellen des Tschebyscheff-Polynoms T11 als Stützstellen

Zeichnen Sie anschließend die Werte der Interpolationspolynome zusammen mit der exakten Funktion in ein Diagramm, wobei diese an den Stellen

{\xi _i} = -1+i \cdot 0,01,\quad \quad i = 0, \ldots ,200

ausgewertet werden sollen.

Lösung

u81.m

clear;
tic

% Exakte Funktion:
xe = -1 : .01 : 1;
fe = 1 ./ (1 + 25 .* xe .* xe);

% Stützstellen mit Abstand 0.2:
x1 = -1 : 0.2 : 1;
y1 = 1 ./ (1 + 25 .* x1 .* x1);
a1 = divdiff(x1, y1);
f1 = ones(length(xe), 1);

% Stützstellen mit variablem Abstand:
x2 = [-1, -6/5 + 1/4, -6/5 + 1/3, -6/5 + 1/2, -6/5 + 1, ...
    6/5 - 1, 6/5 - 1/2, 6/5 - 1/3, 6/5 - 1/4, 1];
y2 = 1 ./ (1 + 25 .* x2 .* x2);
a2 = divdiff(x2, y2);
f2 = ones(length(xe), 1);

% Tschebyscheff-Nullstellen
x3 = 11 : -1 : 1;
x3 = cos((2 .* x3 - 1) / 22 * pi);
y3 = 1 ./ (1 + 25 .* x3 .* x3);
a3 = divdiff(x3, y3);
f3 = ones(length(xe), 1);

for i = 1 : length(xe)
    f1(i) = horner(x1, a1, xe(i));
    f2(i) = horner(x2, a2, xe(i));
    f3(i) = horner(x3, a3, xe(i));
end

plot(xe, fe , xe, f1, xe, f2, xe, f3)
legend('exakt', 'gleiche Abstaende', 'variable Abstaende', 'Tschebyscheff')

toc % Rechenzeit: 0.093344 Sekunden

divdiff.m

% Bestimmung der Koeffizienten mit Hilfe dividierter Differenzen

function a = divdiff(x, y)
% x: Stuetzstellen
% y: Funktionswerte an den Stuetzstellen
% a: Polynomkoeffizienten (dividierte Differenzen)

n = length(x);
f = zeros(n, n);

for i = 1 : n
    f(i,i) = y(i);
end

for k = 2 : n % Zeile
    for i = 1 : n - k + 1
        % disp(['f[' num2str(i) ', ' num2str(i+k-1) '] = (f[' num2str(i+1) ...
        %    ', ' num2str(i+k-1) '] - f[' num2str(i) ', ' num2str(i+k-2) ...
        %    ']) / (x(' num2str(i+k-1) ') - x(' num2str(i) '))'])
        f(i, i+k-1) = (f(i+1, i+k-1) - f(i, i+k-2)) / (x(i+k-1) - x(i));
    end
end

a = f(1, :) ';

horner.m

% Auswertung eines Polynoms mit dem Hornerschema

function f = horner(x, a, xi)
% x: Stuetzstellen des Polynoms
% a: Polynomkoeffizienten (dividierte Differenzen)
% xi: Auswertungsstelle
% f: Ergebnis der Auswertung

f = a(length(a));

for j = length(a) : -1 : 2
    f = f * (xi - x(j - 1)) + a(j - 1);
end