H09 – natürliche Splineinterpolation

 

a )

Schreiben Sie ein m-file (’spline.m’), in dem zu gegebenen Stützstellen {x_0}, \ldots ,{x_n} und entsprechenden Werten {y_0}, \ldots ,{y_n} die Koeffizienten {C_i},{D_i},{M_i} des natürlichen kubischen Interpolationssplines S_4^2 berechnet werden. Verwenden Sie dazu folgenden Programmkopf:

function [C,D,M] = spline(x,y)
%x Stützstellen des Polynoms
%y Funktionswerte an den Stützstellen
%C,D,M Koeffizienten f¨ ur Sˆ2_4 Spline

b )

Schreiben Sie ein m-file (’evaluate.m’), in dem die Splinefunktion aus a) an vorgegebenen Stellen ausgewertet wird. Verwenden Sie dazu folgenden Programmkopf:

function z= evaluate(w,x,C,D,M)
%w Stellen an denen Spline ausgewertet werden soll (aufsteigend sortiert)
%x Stützstellen
%C,D,M Koeffizienten f¨ ur Sˆ2_4 Spline

c )

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

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

im Intervall \left[ {-1,1} \right] interpoliert wird. Verwenden Sie dazu äquidistante Stützstellen

{x_i} = -1+i \cdot 0,2\quad \quad i = 0, \ldots ,10

Zeichnen Sie anschließend die Werte der Interpolationssplines 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

u91.m

clear;
tic

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

% Stützstellen mit Abstand 0.2:
x = -1 : .2 : 1
y = 1 ./ (1 + 25 .* x .* x)

[C D M] = spline(x,y)
z = evaluate(xe, x, C, D, M);

plot(xe, fe , xe, z)
legend('exakt', 'Splines')

toc

spline.m

function [C, D, M] = spline(x,y)
% x: Stuetzstellen des Polynoms
% y: Funktionswerte an den Stuetzstellen
% C, D, M: Koeffizienten des S2_4 Splines

n = length(x) - 1;
h = zeros(n, 1);
for i = 1 : n
    h(i) = x(i + 1) - x(i);
end

% Aufstellen des Gleichungssystems
diag = ones(n - 1, 1);
for i = 1 : n - 1
    diag(i) = 2 * (h(i) + h(i + 1));
end

v2 = h(1 : n - 1);
A = spdiags([v2, diag, v2], -1: 1, n - 1, n - 1);

b = ones(n - 1, 1);
for i = 1 : n - 1
    b(i) = 6 * ((y(i + 2) - y(i + 1)) / h(i + 1) - (y(i + 1) - y(i)) / h(i));
end

M = A \ b;
M = [0; M; 0]; % natürliche Splines

C = zeros(n - 1, 1);
D = zeros(n - 1, 1);
for i = 1 : n
    C(i) = 1/2 * (y(i) + y(i + 1)) - h(i) * h(i) / 12 * (M(i) + M(i + 1));
    D(i) = 1/h(i) * (y(i + 1) - y(i)) - h(i) / 6 * (M(i + 1) - M(i));
end

evaluate.m

function z = evaluate(w, x, C, D, M)
% w: Stellen, an denen Spline ausgewertet werden soll
% x: Stuetzstellen
% C, D, M: Koeffizienten des S2_4 Splines

stellen = length(w);
z = zeros(stellen, 1);

n = length(x) - 1;
h = zeros(n, 1);
for i = 1 : n
    h(i) = x(i + 1) - x(i);
end

akt = 1;

for i = 1 : stellen
    while w(i) > x(akt + 1)
        akt = akt + 1;
    end
    z(i) = C(akt) + D(akt)*(w(i) -  (x(akt) + x(akt + 1)) / 2) ...
        + M(akt + 1) * (w(i) - x(akt))^3 / 6 / h(akt) ...
        - M(akt) * (w(i) - x(akt + 1))^3 / 6 / h(akt);
end

z

Ähnliche Artikel

Kommentar verfassen