Esercizi Informatica

FUNCTION SECANTE VARIABILE

function [xnew,it,vettscarti]=funsecvar(f,x0,x1,toll,itmax)

% function [xnew,it,vettscarti]=funsecvar(f,xo,x1,toll,itmax)

% Input f = funzione data come function handle per risolvere il problema

% f(x)=0

% x0,x1 = approssimazioni iniziale / punti iniziali della successione,

corrisponde all'iterazione 0

% toll = tolleranza prefissata per l'approssimazione

% itmax = n° massimo di iterazioni

% Output xnew = soluzione approssimata/valore approssimato per il punto fisso oppure

% valore distante da esso, approssimazione della soluzione all'iterazione

corrente (ultimo valore che si ottiene dall'algoritmo)

% it = n° di iterazioni effettuate

% vettscarti = vettore degli scarti, fornisce il valore dello scarto ad ogni

iterazione

it=1; % contatore delle iterazioni, parte da 1 perché consideriamo i primi due punti dati

in input

xold=x1;

xoldold=x0;

vettscarti=zeros(itmax,1); % preallocazione di memoria per il vettore degli scarti

scarto=abs(x1-x0); % valore fittizio per poter entrare nel ciclo

% Lo scarto è una buona approssimazione dell'errore

% L'importante è che la quantità iniziale dello scarto sia maggiore della tolleranza

% Pongo lo scarto = ad un valore fittizio per entrare nel ciclo while e applicare

l'algoritmo

% Se lo scarto diventa maggiore della tolleranza --> algoritmo non converge

vettscarti(it)=scarto; % primo valore dello scarto

while scarto>toll && it<itmax

% il ciclo viene eseguito fintantoché lo scarto è maggiore della tolleranza & il n°

di iterazioni è minore del n° massimo di iterazioni

it=it+1; % incremento it di 1, vado avanti di un'iterazione nel metodo

cn=(f(xold)-f(xoldold))/(xold-xoldold);

xnew=xold-f(xold)/cn;

scarto=abs(xnew-xold); % aggiorno lo scarto all'iterazione it

vettscarti(it)=scarto; % aggiorno la componente it del vettore scarti ad ogni

iterazione

xoldold=xold;

xold=xnew; % aggiorno le variabili per poter applicare il passo successivo

end

vettscarti=vettscarti(1:it); % taglio il vettscarti fino alle prime it componenti (quelle

che ci servono) in modo da togliere tutti gli 0 che non mi servono, così avrà la

lunghezza delle iterazioni effettivamente eseguite

% Gli scarti servono per fare il grafico di convergenza

end % chiudo la function SCRIPT SECANTE VARIABILE

Primo caso test con f(x)=0 dove f(x) = e^(-x)-sin(pi x) partendo da 0 e 0.1 in [a,b]=[0 0.5]

Stimare M facendo il rapporto tra gli scarti con p=1.618 come esponente del denominatore

Secondo caso test con f(x) = (x-4)^2 con f(x)=0 trova la radice x=4 che e' doppia cioe' va contata 2 volte in [2 6]

partendo da x0=3 e x1=3.1

Come N-R diventa lineare perché la radice è multipla anche la secante variabile diventa lineare

casotest=2;

switch casotest

case 1

f=@(x) exp(-x)-sin(pi*x);

df=@(x) -exp(-x)-pi*cos(pi*x);

ddf=@(x) exp(-x)+pi^2*sin(pi*x);

a=0; b=0.5; % estremi dell’intervallo

x0=0; x1=0.1;

case 2

f=@(x) (x-4).^2;

df=@(x) 2.*(x-4);

ddf=@(x) 2;

a=2; b=6; % estremi dell’intervallo

x0=3; x1=3.1;

end

fplot(f,[a b]); % grafico della funzione nell’intervallo [a,]

hold on

yline(0,'--'); % grafico delle ascisse

itmax=100; % n° massimo di iterazioni

toll=1.e-16; % tolleranza prefissata

[xnew,it,vettscarti]=funsecvar(f,x0,x1,toll,itmax);

fprintf(1,'iterazioni effettuate = %i \n',it);

fprintf(1,'approssimazione raggiunta = %15.9g \n',xnew);

plot(xnew,f(xnew),'ro');

figure

semilogy(1:it,vettscarti,'o',1:it,vettscarti); % grafico semilogaritmico di convergenza

del metodo con: ascisse i valori delle iterazioni, ordinate i logaritmi degli scarti

% Stimo il fattore di convergenza usando il vettore degli scarti

% M = fattore di convergenza/costante asintotica dell'errore/rapporto tra gli scarti

p=0.618; % fattore di convergenza superlineare

M=(abs(0.5*(ddf(xnew)./df(xnew)))).^p;

fprintf(1,'Approssimazione del valore di convergenza se p=0.618 \n');

fprintf(1,'%g \n',M); INTERPOLAZIONE

SCRIPT POLYVANDER

% Quattro coppie di punti (0 1) (1 2) (1.5 2.5) (2 4) --> n = 4

% Grado nel polinomio m = n-1 = 3

x=[0 1 1.5 2]; % vettore della ascisse dei punti da interpolare

y=[1 2 2.5 4]; % vettore delle ordinate dei punti da interpolare

n=length(x); % n = n° delle coppie di punti da interpolare

m=n-1; % m = grado al più del polinomio da interpolare

A=vander(x); % costruisco la matrice di Vandermonde: le colonne sono potenze del vettore

x

% Dobbiamo risolvere il sistema A*p=y dove p è il vettore incognito che ci darà i

coefficienti del polinomio di interpolazione dal valore corrispondente ad a_n fino ad

arrivare al valore di a_0

y=y(:); % rende colonna un vettore riga e lascia colonna un vettore colonna

p=A\y; % risolvo il sistema lineare

a=min(x); % a = il più piccolo valore del vettore x

b=max(x); % b = il più grande valore del vettore x

xval=linspace(a,b); % creare un vettore di 100 punti equidistanti nell'intervallo di

interpolazione [a,b]

yval=polyval(p,xval); % valuto il polinomio di interpolazione nel vettore xval, cioè

valuto il polinomio nelle 100 componenti del vettore xval

plot(x,y,'ko') % grafico dei punti di interpolazione

hold on

plot(xval,yval,'r--'); % grafico del polinomio

SCRIPT INTERPOLAZIONE POLINOMIALE

Date alcune coppie di dati vediamo cosa succede con polyfit. Interpolazione di una funzione

x=[-1 2 4 6 9]; % vettore della ascisse dei punti da interpolare

f=@(x) exp(-2*x);

y=f(x); % vettore delle ordinate dei punti da interpolare

n=length(x); % n = n° delle coppie di punti da interpolare

m=n-1; % m = grado al più del polinomio da interpolare

p=polyfit(x,y,m); % vettore con i coefficienti del polinomio di interpolazione,

restituisce i coefficienti in ordine di grado decrescente

a=min(x); % a = il più piccolo valore del vettore x

b=max(x); % b = il più grande valore del vettore x

xval=linspace(a,b); % creare un vettore di 100 punti equidistanti nell'intervallo di

interpolazione [a,b]

yval=polyval(p,xval); % valuto il polinomio di interpolazione nel vettore xval, cioè

valuto il polinomio nelle 100 componenti del vettore xval

yval1=f(xval); % valuto il polinomio in xval con la funzione ?

plot(x,y,'ko',xval,yval,'b',xval,yval,'r--');

legend('Nodi di interpolazione','Polinomio di interpolazione','Funzione interpolata');

errore=abs(yval1-yval);

figure

plot(xval,errore); SCRIPT UNISCI I PUNTI

Script che permette di scegliere i punti e di creare la curva che li unisce usando un'interpolazione di tipo lineare

m=input('scrivi quanti punti vuoi prendere '); % m = numero di punti

figure % apro la figura

axis([0 1 0 1]); % apro una figura con 0 <= x <= 1 e 0 <= y <= 1 cioe' il quadrato

[0 1] x [0 1] fissiamo gli assi in modo da prendere i punti su

questa finestra

hold on

x=zeros(m,1); y=zeros(m,1);

for i=1:m

[x(i),y(i)]=ginput(1); % con questo comando scegliamo il punto

plot(x(i),y(1),'r.') % metto il punto scelto come punto sulla figura

text(x(i),y(i),num2str(i)); % vai sul punto che ho scelto e scrivi il numero del

numero: scriviamo che numero di nodo è

(permette di scrivere un testo sulla figura)

end

plot(x,y); % unisce i punti FUNCTION BASE LAGRANGE

function [yval]=fun_baseLagrange(i,x,xval)

% function [yval]=fun_baseLagrange(i,x,xval)

% Input i = indice della funzione base di Lagrange (L1, L2, L3, ..., Lm)indice del

% polinomio base di Lagrange

% x = vettore delle ascisse da interpolare

% xval = vettore in cui si valuta il polinomio base di interpolazione

% (vettore in cui valuto la funzione base)

% Output yval = vettore con il valore del polinomio base i-esimo di Lagrange valutato

% in xval, valuta Li in xval --> Li(xval)

% In questa function gli indici vanno da 1 a n+1 (non da 0 a n) perchè ho n+1 coppie di

punti. Quindi i polinomi della formula di Lagrange da applicare partono dall'indice 1

fino ad arrivare all'indice n e viene tutto traslato di 1 negli indici

% Formula di Lagrange L0(x)= (x-x1)(x-x2)/(x0-x1)(x0-x2) da applicare nel ciclo for

xval=xval(:); % rendere xval vettore colonna

m=length(x); % m = lunghezza del vettore

yval=ones(size(xval)); % vettore di tutti 1, ha = dimensione a xval

for j=1:i-1

yval=(yval.*(xval-x(j)))./(x(i)-x(j)); % faccio la prima parte della produttoria

% se xval e yval sono vettori bisogna vettorizzare le operazioni

end

% Separo i due cicli per saltare l'indice i

for j=i+1:m

yval=(yval.*(xval-x(j)))./(x(i)-x(j)); % faccio la seconda parte della produttoria e

salto l'indice i

% se xval e yval sono vettori bisogna vettorizzare le operazioni

end

end % chiudo la function FUNCTION POLY LAGRANGE

function [yval]=fun_polyLagrange(x,y,xval)

% function [yval]=fun_polyLagrange(x,y,xval)

% Input x = vettore delle ascisse da interpolare

% y = vettore delle ordinate da interpolare

% xval = vettore in cui si valuta il polinomio base di interpolazione

% Output yval = vettore con il valore del polinomio base i-esimo di lagrange valutato

in xval, valuta Li in xval --> Li(xval) pn(xval), il polinomio proprio

% In questa function gli indici vanno da 1 a n+1 (non da 0 a n) perchè ho n+1 coppie di

punti. Quindi i polinomi della formula di Lagrange da applicare partono dall'indice 1

fino ad arrivare all'indice n e viene tutto traslato di 1 negli indici

m=length(x); % m = lunghezza del vettore

n=m-1; % n = grado del polinomio di interpolazione

yval=zeros(size(xval)); % vettore di tutti 0, ha = dimensione di xval

for i=1:n+1

yval=yval+fun_baseLagrange(i,x,xval).*y(i); % sommo tutti i contributi delle funzioni

base di Lagrange moltiplicati per il valore della corrispondente ordinata di

interpolazione e ottengo il polinomio di interpolazione valutato nei punti xval

end

end % chiudo la function SCRIPT POLY LAGRANGE

Si implementi la function poly lagrange che valuta il polinomio di interpolazione di Lagrange.

Salvare i dati in un file da chiamare datiex4.dat disposti su due righe

A=load('dati_ex4.dat');

x=A(1,:);

y=A(2,:);

m=length(x); % m = lunghezza del vettore

n=m-1; % n = grado del polinomio

p=polyfit(x,y,n); % vettore con i coefficienti del polinomio di interpolazione,

restituisce i coefficienti in ordine di grado decrescente

xval=linspace(min(x),max(x));

yval=polyval(p,xval); % polinomio di interpolazione

yvalI=interp1(x,y,xval); % interpolazione a tratti

yvalS=spline(x,y,yval); % spline

plot(x,y,'ro',xval,yval,'k--',xval,yvalI,'b',xval,yvalS,'g')

legend('dati','polinomio',