Questa è la funzione:
function f = fun(x)
f = 2*x.*exp(-x) + cos(2*x);
Il codice rimanente dell'esercizio è questo:
n=8;
xx=linspace(0,2*pi,50);
yv=fun(xx);
x=linspace(0,2*pi,n);
y=fun(x);
pN= newton(x,y,xx);
[a,b,pS] = splineper3(x,y,xx);
x=linspace(0,2*pi-2*pi/n,n-1);
y=fun(x);
pF= interpft(y,50);
figure(1)
plot(xx,yv,'k',xx,pN,'b',xx,pS,'g',xx,pF,'r');
title('Confronto grafici');
legend('Funzione esatta','Newton','Spline perodica','Interpolazione trigonometrica');
Nel codice vengono richiamate le seguenti funzioni descritte precedentemente nel progetto:
function p = newton(x,y,xx)
n = length(x-1);
m = length(xx);
p = zeros(1,m);
diff = diag(diffdivise(x,y));
% calcolo del polinomio di Newton
for l = 1 : m
x0 = xx(l);
p0 = diff(n);
% Metodo Ruffini-Horner
for j = n:-1:2
p0 = p0*(x0-x(j-1)) + diff(j-1);
end
p(l) = p0;
end
function D = diffdivise(x,y)
n = length(x-1);
D = zeros(n);
D(:,1) = y';
for j=2:n
for i=j:n
D(i,j) = (D(i,j-1)-D(i-1,j-1)) / (x(i)-x(i-j+1));
end
end
function [xxs,yys,yval] = splineper3(x,y,xval)
n = length(x);
clf;
plot(x,y,'or');
%calcolo del nuovo nodo
x(n+1) = x(n)+x(2)-x(1);
y(n+1) = y(2);
% calcolo distanze dei nodi
h = NaN;
for i=1:n
h(i) = x(i+1)-x(i);
end
% calcolo delle costanti gamma(i) e delta(i)
gamma = NaN;
delta = NaN;
gamma(1) = NaN;
delta(1) = NaN;
for i=2:n
gamma(i) = h(i-1)/(h(i-1)+h(i));
delta(i) = h(i)/(h(i-1)+h(i));
end
% calcolo della matrice A
A = zeros(n-1,n-1);
for i=1:n-1
A(i,i) = 2;
end
A(1,n-1) = gamma(2);
A(n-1,1) = delta(n);
for i=1:n-2
A(i,i+1) = delta(i+1);
end
for i=2:n-1
A(i,i-1) = gamma(i+1);
end
% calcolo delle differenze divise f[x(i), x(i+1), x(i+2)]
f = NaN;
for i=1:n-1
nodix = [x(i) x(i+1) x(i+2)];
nodiy = [y(i) y(i+1) y(i+2)];
c = diffdivise(nodix,nodiy);
f(i) = 6*c(3,3);
end
% calcolo dei momenti
mu = NaN;
muaux = A\f';
mu(1) = 0;
for i=1:n-1
mu(i+1) = muaux(i);
end
mu(1) = muaux(end);
% calcolo dei coefficienti alfa(i) e beta(i)
alfa = NaN;
beta = NaN;
for i=1:n-1
beta(i) = y(i) - mu(i)*h(i)^2/6;
alfa(i) = (y(i+1)-y(i))/h(i) - h(i)*(mu(i+1)-mu(i))/6;
end
% creazione e disegno delle cubiche che compongono la spline periodica
xxs = zeros(n-1,100);
yys = zeros(n-1,100);
hold on
title('Spline cubica periodica');
for j=1:n-1
xxs(j,:) = linspace(x(j),x(j+1),100);
for i=1:100
yys(j,i) =(((x(j+1)-xxs(j,i))^3)*mu(j))/(6*h(j)) + (((xxs(j,i)-x(j))^3)*mu(j+1))/(6*h(j)) + alfa(j)*(xxs(j,i)-x(j)) + beta(j);
end
plot(xxs(j,:),yys(j,:),'-');
end
% calcolo dei valori della spline sui punti del vettore xval
yval = NaN;
if nargin<3
return
else
for i=1:length(xval)
if xval(i)>x(n) || xval(i)<x(1) %se xval è fuori dall'intervallo
yval(i) = NaN;
else
int = 1;
for j=2:n-1
if xval(i)>x(j) && xval(i)<=x(j+1)
int = j;
end
end
yval(i) =(((x(int+1)-xval(i))^3)*mu(int))/(6*h(int)) + (((xval(i)-x(int))^3)*mu(int+1))/(6*h(int)) + alfa(int)*(xval(i)-x(int)) + beta(int);
plot(xval(i),yval(i),'*r');
end
end
end
hold off
Con tutto questo codice dovrebbe partire il programma e quindi uscire anche il grafico.
Le altre funzioni però dovrebbero stare bene perchè i grafici escono corretti.
L'errore sta sicuramente nell'utilizzo dell'iterpft perchè nel grafico quella è l'unica curva a non passare per tutti i punti.