Parte 9¶

Sperimentazione numerica relativa ll'interpolazione di dati e funzioni.

Per stampare il grafico corrente si può impiegare: print -dpng foo.png

Interpolante di Lagrange¶

Calcola il polinomio di interpolazione in forma di Lagrange

yy = interpLagrange(x, y, xx)

Input:

x: vettore dei nodi di interpolazione, con elementi distinti
y: vettore dei valori assunti dalla funzione nei nodi di interpolazione,
xx: vettore dei punti in cui si vuole valutare il polinomio interpolante.

Output:

yy: vettore contenente i valori assunti dal polinomio interpolante,

In [1]:

function [yy] = interp_lagrange(x, y, xx)

n = length(x);

for j = 1:n
    % calcolo coefficienti j-esimo polinomio fondamentale di Lagrange
    x_zeri = [x(1:j-1) x(j+1:end)];
    p = poly(x_zeri);
    p = p / polyval(p, x(j));
    % calcolo valori assunti dal polinomio
    L(j,:) = polyval(p, xx);
end

yy = y * L;

end

Interpolante di Newton¶

Calcola il polinomio di interpolazione in forma di Newton. Per la valutazione del polinomio si impiega l'algoritmo di Horner.

yy = interpNewton(x, y, xx)

Input:

x: vettore dei nodi di interpolazione con elementi distinti
y: vettore dei valori assunti dalla funzione nei nodi di interpolazione,
xx: vettore dei punti in cui si vuole valutare il polinomio interpolante.

Output:

yy: vettore contenente i valori assunti dal polinomio interpolante,

In [2]:

function [yy] = interp_newton(x, y, xx)

% FASE 1 creazione della tabella delle differenze divise
% Calcolo dei coefficienti del polinomio di Newton

n = length(x);
%for k = 2 : n % si calcolano n - 1 colonne (una per coeff.)
 %   y(k:n) = (y(k:n) - y(k-1:n-1))./(x(k:n) - x(1:n-k+1));
%end
for k = 1 : n-1
    y(k+1:n) = (y(k+1:n) - y(k)) ./ (x(k+1:n) - x(k));
end

coeff = y; % diagonale della tabella delle differenze divise.

% si ottiene il vettore dei coefficienti per cui vale p(x(i)) = y(i)
% dove p(x) = c(1) + c(2)(x - x(1)) + ... + c(n)(x - x(1))...(x - x(n - 1))

% FASE 2 applicazione algoritmo di Horner
% Calcolo dei valori del polinomio di Newton nel vettore y.

n = length(coeff);
yy = coeff(n) * ones(size(xx)); 

for k = n - 1: -1 : 1
    yy =  (xx - x(k)).*yy + coeff(k);
end


% Si ottiene yy(i) = p(xx(i))
% dove p(x) = c(1) + c(2)(x - x(1)) + ... + c(n)(x - x(1))...(x - x(n - 1))

end

Esercizio 1¶

In [3]:

x = [0, 1/4, 1/2, 3/4, 1];
%x = [-1 -0.7 0.5 2];
y = [];

xx = linspace(min(x), max(x));
cols = ['b-', 'c-', 'r-', 'g-', 'r-'];

n = length(x);

for j = 1:n
    % calcolo coefficienti j-esimo polinomio fondamentale di Lagrange
    x_zeri = [x(1:j-1) x(j+1:end)];
    p = poly(x_zeri);
    p = p / polyval(p, x(j));
    % calcolo valori assunti dal j-esimo polinomio
    L(j,:) = polyval(p, xx);
    % rappresentazione grafica
    figure(j)
    hold on
    plot(xx, L,cols(j))
    plot(x, zeros(1,n), 'b-')
    plot(x(j), 1, 'c*')
end

Esercizio 2¶

In [4]:

for k = 0:4
    x(k + 1) = (k * pi) / 2;
    y(k + 1) = sin(x(k + 1));
end

xx = linspace(min(x), max(x));
yy = interp_lagrange(x, y, xx);
figure(1)
plot(xx, yy, "r-", x, y, '*', xx, sin(xx), 'b-')
legend("interpolante di Lagrange","punti di interpolazione", "y = sin(x)");
box off

In [5]:

for k = 0:4
    x(k + 1) = (k * pi) / 2;
    y(k + 1) = cos(x(k + 1));
end

xx = linspace(min(x), max(x));
yy = interp_lagrange(x, y, xx);
figure(2)
plot(xx, yy, "r-", x, y, '*', xx, cos(xx), 'b-')
legend("interpolante di Lagrange","punti di interpolazione", "y = cos(x)");
box off

Esercizio 3¶

La temperatura T in prossimità del suolo subisce una variazione dipendente dalla latitudine L secondo la seguente tabella:

$$ \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline L & -55 & -45 & -35 & -25 & -15 & -5 & 5 & 15 & 25 & 35 & 45 & 55 & 65 \\ \hline T & 3.7 & 3.7 & 3.52 & 3.27 & 3.2 & 3.15 & 3.15 & 3.25 & 3.47 & 3.52 & 3.65 & 3.67 & 3.52 \\ \hline \end{array} $$

Si vuole costruire un modello che descriva la legge $T = T(L)$ anche per latitudini non misurate. A tal fine si scriva uno script che fornisca la variazione di temperatura alle latitudini $L = \pm 42$ utilizzando il polinomio interpolante.

Visualizzare in un grafico i dati assegnati, il polinomio interpolante e le stime di T ottenute per $L = \pm 42.$

In [26]:

x = [-55 -45 -35 -25 -15 -5 5 15 25 35 45 55 65];
y = [3.7 3.7 3.52 3.27 3.2 3.15 3.15 3.25 3.47 3.52 3.65 3.67 3.52];

xx = linspace(min(x), max(x), 200);
figure(3)
hold on
[yy] = interp_newton(x, y, xx);
l = length(yy);
plot(xx, yy, 'b');
plot(x, y, 'r*'); % dati sperimentali
plot(42, yy(1), 'go'); % stima L = 42°
plot(-42, yy(l), 'ro'); % stima L = -42°
legend("interpolante di Newton", "dati", "stima L = 42°", "stima L = -42°", ...
    "location", "southoutside")

Esercizio 4¶

Scrivere uno script che calcoli il polinomio interpolante un insieme di punti $P_i = (x_i , y_i), i = 1, \dotsc, n + 1$, nella forma di Newton con $x_i$ scelti dall’utente come:

punti equidistanti in un intervallo $[a, b]$,
punti definiti dai nodi di Chebyshev nell’intervallo $[a, b]$, ossia

$$ x_{i}=\frac{(a+b)}{2}+\frac{(b-a)}{2} \cos \left(\frac{(2 i-1) \pi}{2(n+1)}\right), \quad i=1, \ldots, n+1 $$

e y i = f (x i ) ottenuti dalla valutazione nei punti x i di una funzione test f : [a, b] → R. Testare lo script sulle funzioni

$f(x) = \sin(x) - 2sin(2x)$,
$f(x) = \sinh(x)$,
$f(x) = |x|$,
$x \in [-\pi, \pi]$,
$x \in [-2, 2]$,
$x \in [-1, 1]$,
$f(x) = 1/(1 + x_2)$,
$x \in [-5, 5]$ (funzione di Runge).

Calcolare l’errore di interpolazione $r(x) = f (x) - p(x)$, tra la funzione test $f(x)$ e il polinomio di interpolazione $p(x)$. Visualizzare il grafico di $f(x)$ e $p(x)$, ed il grafico di $|r(x)|$. Cosa si osserva? Cosa accade all’aumentare del grado $n$ di $p(x)$?

In questo esercizio si sfrutta l'interpolazione non per dedurre dati sperimentali ma per approssimare funzioni. Il grado generalmente è pari al numero di punti forniti meno uno ma in questo caso l'utente può scegliere anche il grado (per studiare la variazione della norma della funzione resto).

In [7]:

text = ["1. Scegli il grado: troppo basso o troppo alto" ...
    " e l'approssimazione non sarà buona"];
disp(text);
n = input(" grado: ");

text = ["2. Indica l'intervallo [a,b] entro cui estrarre" ... 
    " i nodi di interpolazione:"];
disp(text);
a = input(" a: ");
b = input(" b: ");

text = ["3. Definisci la strategia di calcolo dei nodi di" ... 
    " interpolazione \n [1] Equispaziati \n [2] Chebyshev"];

disp(text);
scelta = input(" strategia: ");
switch scelta
    case 1 % equispaziati
        x = linspace(a, b, n + 1);
    case 2 % Chebyshev
        for i = 1 : n + 1
            x(i) = (a + b)/2 + (b - a)/2 * cos(((2*i - 1)) ...
                * pi / ((2 * (n + 1))));
        end
    otherwise
        close all
        return
end

text = ["4. Scegli la funzione da interpolare \n" ...
    " [1] sin(x) - 2sin(2x) \n [2] sinh(x) \n" ... 
    " [3] |x| \n [4] 1 / (1 + x^2)"];
disp(text);
scelta = input(" funzione: ");

xx = linspace(a, b, 301);
switch scelta
    case 1
        f = sin(xx) - 2 * sin(2 * xx); y = sin(x) - 2 * sin(2 * x);
    case 2
        f = sinh(xx); y = sinh(x); 
    case 3
        f = abs(xx); y = abs(x); 
    case 4
        f = 1 ./ (1 + xx.^2); y = 1 ./ (1 + x.^2);
    otherwise
      close all
      return
end


yy = interp_newton(x, y, xx);

figure(1)
plot(xx, yy, "b-", x, y, "ro", xx, f, "r--");
legend("Interpolante di Newton", "Dati", "Funzione test", ...
    "location", "southoutside");
box off

1. Scegli il grado: troppo basso o troppo alto e l'approssimazione non sarà buona

2. Indica l'intervallo [a,b] entro cui estrarre i nodi di interpolazione:

3. Definisci la strategia di calcolo dei nodi di interpolazione 
 [1] Equispaziati 
 [2] Chebyshev

4. Scegli la funzione da interpolare 
 [1] sin(x) - 2sin(2x) 
 [2] sinh(x) 
 [3] |x| 
 [4] 1 / (1 + x^2)

In [8]:

f_resto = abs(f - yy);
figure(2);
plot(xx, f_resto, "b-");
legend("abs err. interpolazione")
box off

In [9]:

norm_inf = max(f_resto)

norm_inf =  0.000044177

Esercizio 5¶

In [10]:

clc

a = -1; 
b = 1; 
spacing = 200;
xx = linspace(a, b, spacing);

LLe = zeros(1,4);
LLc = zeros(1,4);

i=0;

for n = 5:5:20  
    i = i + 1;
    
    xe = linspace(a, b, n + 1);
    xc = cos(((2*[0:n] + 1)) * pi / ((2 * (n + 1))));
    
    Le = zeros(1, spacing);
    Lc = zeros(1, spacing);
    
    for j = 1 : n + 1
        x_zeri = [xe(1:j-1) xe(j+1:end)];
        p = poly(x_zeri);
        pe = p / polyval(p, xe(j));
    
        Le = Le + abs(polyval(pe, xx));
        
        x_zeri = [xc(1:j-1) xc(j+1:end)];
        p = poly(x_zeri);
        pc = p / polyval(p, xc(j));
        
        Lc = Lc + abs(polyval(pc, xx));
    end
    
    LLe(i) = max(Le);
    LLc(i) = max(Lc);
end

LLe
LLc

LLe =

       3.10494      29.89431     508.71131   10759.64903

LLc =

   2.1044   2.4894   2.7278   2.9008

Si nota come la costante di Lebesque assuma valori notevolmente più grandi quando si impiegano nodi equispaziati mentre quando si utilizzano i nodi di Chebyshev la costante risulta molto vicina all'unità. Si deduce che con il set di dati fornito l'indice di condizionamento è particolarmente elevato nel caso della scelta di nodi equispaziati. Qualsiasi algoritmo si scelga di impiegare, sarà difficile ottenere delle soluzioni accurate.

Esercizio 6¶

In [19]:

clc
clearvars
a = -1; b = 1;
x = linspace(a, b, 22);

xx = linspace(min(x), max(x), 301);
f = @(z) sin(2 * pi * z);

yy_es = f(xx);

y = f(x);
yy = interp_lagrange(x, y, xx);

y_t = y + 0.0002 * randn(1, 22);
yy_t = interp_lagrange(x, y_t, xx);

In [20]:

figure
plot( xx, yy, "r-", ...
      x, y, "bo", ...
      xx, yy_es, "g--")
title("Interpolante di Newton con nodi equispaziati")
legend("intepolante di Newton", "dati sperimentali", "funzione test")
box off

In [21]:

figure
plot( xx, yy_t, "r-", ...
      x, y, "co", ...
      x, y_t, "bo", ...
      xx, yy, "g--");
title("Interp. di Newton con nodi equispaziati e perturbazione sui dati")
legend("intep. perturbato di Newton", "dati", ...
    "dati perturbati", "interp. Newton", "location", "southoutside")
box off

Con il set di dati perturbato si verifica il fenomeno di Runge.

In [14]:

format long
err_rel_dati = norm(y_t - y, inf) / norm(y, inf)

err_rel_dati =    3.932814314593389e-04

In [15]:

err_rel_ris = norm(yy_t - yy, inf) / norm(yy, inf)

err_rel_ris =    2.920313467184383e-01

Esercizio 7¶

Si calcoli il polinomio interpolatore nella forma di Newton di grado $n = 59$ della funzione $f(x) = log(x)$ nei nodi $x_i = 10 + i/59$, $i = 0, 1, \dotsc, 59$. Si calcoli poi il polinomio interpolatore di grado $n = 59$ per i nodi $x_i = 10 + \frac{i}{59}, i = 59, \dotsc, 1, 0$, ossia ordinati in ordine inverso.

Per entrambi i polinomi si visualizzi in un grafico in scala semilogaritmica il comportamento dell’errore di interpolazione in valore assoluto, commesso in $1000$ punti equispaziati sull’intervallo $[10, 11]$.

In [62]:

xx = linspace(10, 11, 1000);
f = log(xx);

x = linspace(10,11,60);
y = log(x);
yy = interp_newton(x, y, xx);

figure
plot( xx, yy,...
      xx, log(xx), ...
      x,y, "*")
title("Confronto interpolatori di Newton");
box off

In [54]:

x_t = fliplr(x); % flip array left to right
y_t = log(x_t);
yy_t = interp_newton(x_t, y_t, xx);

figure
plot( xx, yy_t,"r", ...
      xx, log(xx), ...
      x,y, "o")
title("Confronto interpolatori di Newton con nodi invertiti");
box off

In [55]:

err_ass = abs(f - yy);
err_ass_t = abs(f - yy_t);

figure
semilogy( xx, err_ass, "b.-", ...
          xx, err_ass_t, "r.-")
warning ("off", "Octave:negative-data-log-axis");
box off