Appunti Calcolo Numerico via HW2

ylabel('Ordinata y(t)')

%%

% Svolti i due passaggi di intermezzo, si calcola, come richiesto, la

% lunghezza del perimetro con la function sopra definita. Inoltre, si

% utilizzano i due vettori _xSpline_ e _ySpline_ come argomenti della

% plot. In questo modo è possibile ottenere il polinomio a tratti di terzo

% grado che interpola sui nodi, i quali sono anch'essi mostrati (tramite

% pallini verdi).

lunghezzaSplineCubica=perimetro(xSpline,ySpline);

plot(xSpline,ySpline,'r',x,y,'og')

xlabel('Ascisse')

ylabel('Ordinate')

%% 1 - Parte seconda: Interpolazione con cubica di Hermite

% Per interpolare con la cubica di Hermite, il ragionamento da seguire è il

% medesimo applicato all'interpolazione tramite spline. L'unica differenza

% consiste nella function da utilizzare che questa volta corrisponde a

% _pchip()_. Quindi, si procede prima con l'interpolazione delle ascisse

% rispetto al solito parametro _t_ e se ne mostra il grafico:

xHermite=pchip(t,x,tt);

plot(t,x,'oc',tt,xHermite,'c')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro t')

ylabel('Ascissa x(t)')

%%

% Valutata la cubica di Hermite sulle ascisse, si applica lo stesso

% ragionamento per le ordinate:

yHermite=pchip(t,y,tt);

plot(t,y,'oc',tt,yHermite,'c')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro t')

ylabel('Ordinate y(t)')

%%

% Anche in questo caso, si calcola la lunghezza della cubica di Hermite che

% ricostruisce il grafico tramite la function _perimetro()_. Infine, viene

% visualizzato il grafico passando come parametri di input della plot i due

% valori _xHermite_ e _yHermite_.

lunghezzaCubicaHermite=perimetro(xHermite, yHermite);

plot(xHermite,yHermite,'c',x,y,'og')

xlabel('Ascisse')

ylabel('Ordinate')

%% 2 - Parte prima: Approssimazione dei minimi quadrati

% Sarà ora necessario procedere con l'approssimazione della curva nel senso

% dei minimi quadrati. In questa prima parte, per approssimare si utilizzerà

% il polinomio di grado 6. Anche in questo caso, trattandosi di una curva

% chiusa, sarà necessario procedere approssimando in funzione del parametro

% _t_ prima le ascisse e poi le ordinate.

% Per evitare intrecci nel grafico finale, è stata definita una nuova griglia

% fitta _ttMQ_. I coefficienti sono stati calcolati tramite _polyfit()_, a

% cui sono stati passati come parametri:

%

% t, il solito parametro

%

% x, il vettore delle ascisse

%

% 6, il grado del polinomio che approssimerà la soluzione.

%

% È stata poi utilizzata la _polyval()_ che si è occupata di valutare il

% polinomio di grado 6, con i coefficienti precedentemente calcolati, sulla

% griglia di valutazione _ttMQ_. Infine, ne è stato visualizzato il grafico:

ttMQ=linspace(t(1),t(end-1),500);

coefficientiX=polyfit(t,x,6);

polinomio6MQx=polyval(coefficientiX,ttMQ);

plot(t,x,'ob',ttMQ,polinomio6MQx,'b')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro t')

ylabel('Ascissa x(t)')

%%

% Analogamente, si procederà allo stesso modo con il vettore delle

% ordinate, infatti:

coefficientiY=polyfit(t,y,6);

polinomio6MQy=polyval(coefficientiY,ttMQ);

plot(t,y,'ob',ttMQ,polinomio6MQy,'b')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro t')

ylabel('Ordinata y(t)')

%%

% I valori contenuti in _polinomio6MQx_ e in _polinomio6MQy_ permettono di

% rappresentare una curva che però resta aperta. Per questo motivo, sarà

% necessario unire i punti con il procedimento descritto all'inizio di

% questa relazione.

polinomio6MQx=[polinomio6MQx polinomio6MQx(1)];

polinomio6MQy=[polinomio6MQy polinomio6MQy(1)];

%%

% Adesso, è possibile calcolare, come richiesto, il perimetro del polinomio

% di sesto grado che approssima la curva, per poi visualizzarne tramite

% plot il grafico:

lunghezzaPolinomio6MQ=perimetro(polinomio6MQx,polinomio6MQy);

plot(polinomio6MQx,polinomio6MQy,'b',x,y,'og')

xlabel('Ascisse')

ylabel('Ordinate')

%% 2 - Parte seconda: Approssimazione con modello trigonometrico

% È ora possibile procedere con l'approssimazione tramite il modello

% trigonometrico

% a1+a2*cos(x)+a3*cos(2x)+a4*cos(3x)+a5*sin(x)+a6*sin(2x)+a7*sin(3x).

% Si ricorda che, anche in questo caso, sarà necessario procedere prima sul

% vettore delle ascisse e poi su quello delle ordinate. Inoltre, come

% accadeva per il polinomio di sesto grado, anche qui si utilizzerà la griglia

% _ttMQ_ per evitare intrecci di sorta. Per approssimare correttamente con

% il modello trigonometrico di ordine 3, sarà necessario ricostruire la

% matrice (unica) B dei coefficienti. Sarà poi possibile calcolare le a(i)

% e valutare il polinomio sulla griglia fitta _ttMQ_.

% Sebbene in Matlab sia possibile risolvere un sistema rettangolare con il

% semplice comando _B\x_, si è scelto di scrivere in forma "estesa" le

% operazioni necessarie alla risoluzione del problema (nel senso dei minimi

% quadrati). Calcolato _xTrig_, si può visualizzare il grafico di _t_ e

% _x(t)_.

B=[ones(size(t)) cos(t) cos(2*t) cos(3*t) sin(t) sin(2*t) sin(3*t)];

a=(B'*B)\(B'*x);

xTrig= a(1)+a(2)*cos(ttMQ)+a(3)*cos(2*ttMQ)+a(4)*cos(3*ttMQ)+a(5)*sin(ttMQ)...

+a(6)*sin(2*ttMQ)+a(7)*sin(3*ttMQ);

plot(t,x,'ok',ttMQ,xTrig,'k')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro t')

ylabel('Ascissa x(t)')

%%

% Per procedere con le ordinate, quindi, sarà necessario ripartire dal

% calcolo dei coefficienti _a_ (la matrice B dei coefficienti del sistema

% rettangolare è la medesima). Analogamente, si ricalcolerà il polinomio

% trigonometrico sulla griglia fitta _ttMQ_ con il nuovo vettore _a_ e si

% visualizzerà questo secondo grafico (_t_ e _y(t)_).

a=(B'*B)\(B'*y);

yTrig= a(1)+a(2)*cos(ttMQ)+a(3)*cos(2*ttMQ)+a(4)*cos(3*ttMQ)+a(5)*sin(ttMQ)...

+a(6)*sin(2*ttMQ)+a(7)*sin(3*ttMQ);

plot(t,y,'ok',ttMQ,yTrig','k')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro t')

ylabel('Ordinate y(t)')

%%

% Prima di procedere al calcolo del perimetro e alla visualizzazione finale

% del grafico, si ricollegano i punti in modo da avere una curva chiusa.

xTrig=[xTrig xTrig(1)];

yTrig=[yTrig yTrig(1)];

%%

% È ora possibile calcolare il perimetro con la function apposita e

% utilizzare il comando plot con input _xTrig_ e _yTrig_ per visualizzare

% il grafico del polinomio trigonometrico approssimante.

lunghezzaPolTrig=perimetro(xTrig,yTrig);

plot(xTrig,yTrig,'k',x,y,'og')

xlabel('Ascisse')

ylabel('Ordinate')

%% 3 - Interpolazione polinomiale sui nodi di Chebychev

% Infine, occorre interpolare i punti con il polinomio di grado 39 sui nodi

% di Chebychev. Prima di procedere, si ricorda per l'ennesima volta che

% sarà necessario applicare il ragionamento che a breve sarà descritto sia

% sulle ascisse sia sulle ordinate.

% Si dovranno generare 40 nodi di Chebychev nell'intervallo [-1,1] con

% l'apposita formula. Questi punti costituiranno un nuovo vettore per la

% valutazione di ascisse e ordinate di Chebychev per l'interpolazione

% parametrica (chiamato _tCheby_). Tali coordinate di Chebychev saranno

% ottenute con la function _pchip()_. I parametri di input, saranno:

%

% t, il parametro per interpolare

%

% x, il vettore delle ascisse

%

% tCheby, tale parametro è utilizzato come griglia di valutazione di 40

% punti per ottenere esattamente 40 ascisse di Chebychev.

%

% Ottenute queste nuove ascisse, quindi, si potranno calcolare i

% coefficienti del polinomio di grado 39 con la _polyfit()_, che costituirà

% il primo dei due input per _polyval()_. _polyfit()_ è utilizzata

% prendendo in input:

%

% tCheby, il parametro

%

% xCheby, il vettore delle ascisse di Chebychev

%

% 39, il grado del polinomio.

%

% _polyval()_, invece, valutato sulla griglia fitta _tt_, restituisce un

% vettore contenente i valori del polinomio di grado 39 che interpola le

% ascisse al variare di _tCheby_. Se ne stampa poi il grafico:

tCheby=-cos((0:39)*pi/39)';

xCheby=pchip(t,x,tCheby);

polChX=polyval(polyfit(tCheby,xCheby,39),tt);

plot(tCheby,xCheby,'oy',tt,polChX,'y')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro tCheby')

ylabel('Ascissa x(tCheby)')

%%

% In modo totalmente analogo, si ripete il procedimento con il vettore

% delle ordinate, come è possibile vedere di seguito:

yCheby=pchip(t,y,tCheby);

polChY=polyval(polyfit(tCheby,yCheby,39),tt);

plot(tCheby,yCheby,'oy',tt,polChY,'y')

axis([-1.1 1.1 -1.1 1.1])

xlabel('Parametro tCheby')

ylabel('Ordinate y(tCheby)')

%%

% Anche nel caso del polinomio di grado 39 che interpola sui nodi di

% Chebychev, la curva risultante nella maggioranza dei casi (come questo)

% non sarà chiusa. Pertanto, si ricongiungono i punti con il solito metodo:

polChX=[polChX polChX(1)];

polChY=[polChY polChY(1)];

%%

% Si calcola poi il perimetro del polinomio interpolante di Chebychev di

% grado 39 con la function apposita e se ne mostra il grafico passando come

% argomenti della plot _polChX_ e _polChY_. Sul grafico, inoltre, sono

% mostrati sia i nodi di Chebychev (quelli segnati da un triangolo rosso),

% sia i nodi da interpolare di partenza (segnati da un pallino verde).

lunghezzaPolCheby=perimetro(polChX,polChY);

plot(polChX,polChY, 'y',x,y,'og', xCheby, yCheby, '^r')

xlabel('Ascissa')

ylabel('Ordinata')

%% Conclusioni - Il grafico

% In questa ultima sezione, come già anticipato, si mostreranno i grafici di

% tutte le ricostruzioni eseguite nei punti precedenti, in modo da

% osservarne le differenze e trarne le opportune considerazioni.

plot(x,y,'og',x,y,'g')

hold on

plot(xSpline,ySpline,'r')

plot(xHermite,yHermite,'c')

plot(polinomio6MQx,polinomio6MQy,'b')

plot(xTrig,yTrig,'k')

plot(polChX,polChY, 'y',xCheby, yCheby, '^r')

title('Grafico completo delle curve')

legend('nodi','Spezzata','Spline','Hermite','Pol 6°', 'Pol Trig',...

'Chebychev', 'Nodi Chebychev', 'Location','BestOutside')

xlabel('Ascisse')

ylabel('Ordinate')

%%

% Come &eg