Metodi Numerici con Elementi di Programmazione · 2019. 12. 14. · Scrivere la funzione Matlab...

Post on 16-Oct-2020

5 views 0 download

Transcript of Metodi Numerici con Elementi di Programmazione · 2019. 12. 14. · Scrivere la funzione Matlab...

Introduzione al MatLab

(9)

1

Metodi Numerici con Elementi di Programmazione

A.A. 2019-2020

2

Docente: Vittoria Bruni Email: vittoria.bruni@sbai.uniroma1.it Ufficio: Via A. Scarpa, Pal. B, I piano, Stanza n. 16 Tel. 06 49766648 Ricevimento: consultare la sezione Avvisi nella pagina web dedicata al corso

3

Il materiale didattico è disponibile sul sito http://ingaero.uniroma1.it/ nella pagina dedicata al corso Metodi Numerici con Elementi di Programmazione Per consultazione: Getting Started with MatLab – The mathworks www.mathworks.com

Esercizi

4

5

Esercizio 1

6

function [P] = successione(a,M) %[P] = successione(a,M) calcola i primi M elementi della successione così definita % P(1) = 2*a % P(2) = 4*a^2-1 % P(k) = 2*a*P(k-1)-P(k-2), k>2 % % e la confronta con la successione così definita % Q(k) = 1/(2*sqrt(a^2-1) )* (abs(a+sqrt(a^2-1))^(k+1) - abs(a-sqrt(a^2-1))^(k+1)); % % % INPUT % a = numero reale tale che |a|>1. La funzione offre tre possibilità per % correggere il valore della variabile data in input, altrimenti % interrompe l'esecuzione % M = numero intero positivo % OUTPUT % P = vettore contenente gli elementi della successione Pn

7

% controllo del numero delle variabili di input e di output if (nargin<=1) | (nargin>2) error('il numero degli input non e'' corretto!!!') end if (nargout>1) error('il numero degli output non e'' corretto!!!') end % controllo del primo input tentativi = 0; while abs(a)<=1 & tentativi <3 a= input('il modulo di a deve essere maggiore di 1! a= '); tentativi = tentativi + 1; end

Nota: All’uscita del ciclo while si possono avere le due situazioni seguenti:

• Il valore di a è in modulo maggiore di 1 e tentativi<=3 • Il numero di tentativi ha raggiunto il valore 3 e il valore

di a è in modulo minore o uguale a 1

8

if abs(a)>1 % controllo del secondo input if (round(M)~=M) | (M<=0) error('il secondo input deve essere un numero intero positivo!!!') end % calcolo del vettore P P(1) = 2*a; P(2) = 4*a^2-1; for k=3:M P(k)=2*a*P(k-1)-P(k-2); end % calcolo del vettore Q Q(1:M) = 1/(2*sqrt(a^2-1)) * (abs(a+sqrt(a^2-1)).^(2:M+1) - abs(a-sqrt(a^2-1)).^(2:M+1));

% grafico dei risultati figure, plot(P,'k-') hold on, plot(Q,'r--') xlabel('k'), ylabel('successioni'), title('confronto successioni') legend('vettore P','vettore Q') else error('gli input inseriti non sono corretti!') end

9

Dal Command Window

>> [P] = successione(2,-10); ??? Error using ==> successione at 38 il secondo input deve essere un numero intero positivo!!! >> [P] = successione(2,1.5); ??? Error using ==> successione at 38 il secondo input deve essere un numero intero positivo!!! >> [P] = successione(2,-1.6); ??? Error using ==> successione at 38 il secondo input deve essere un numero intero positivo!!! >> [P] = successione(6); ??? Error using ==> successione at 22 il numero degli input non e' corretto!!!

10

Dal Command Window

>> [P] = successione(-.5,100); il modulo di a deve essere maggiore di 1! a= .5 il modulo di a deve essere maggiore di 1! a= .2 il modulo di a deve essere maggiore di 1! a= .4 ??? Error using ==> successione at 53 gli input inseriti non sono corretti! >> [P] = successione(2,10); >> [P] = successione(2,100); >> [P] = successione(6,100); >> [P] = successione(-6,100);

A queste chiamate la funzione non da messaggi di errore e grafica i risultati

11

Si consideri la funzione definita per Scrivere la funzione Matlab esercizio2.m che riceva in input un numero reale y e un numero naturale N10 e restituisca in output: • la variabile S contenente il valore di F(y) ottenuto usando i primi N

termini della serie; • il vettore C dei coefficienti fn usati. Se N non viene dato in input, porre

N=1000. La funzione deve graficare le componenti del vettore C usando una linea continua rossa per i primi 10 elementi e un marcatore di punto a scelta di colore blu per i successivi elementi. Il grafico deve essere completo di titolo, etichette per gli assi e legenda.

n=1

Esercizio 2

12

function [S,C] = esercizio2(y,varargin) % function [S,C] = esercizio2(y,varargin) % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. if nargin==1 N=10000; elseif nargin==2 N = varargin{1}; if N<10 | (round(N) ~=N) error(‘il valore di N non e’’ ammissibile’); end

else

error('controllare il numero delle variabili date in input!!!')

end

C(1:2) = 1;, S = 0;

for k=1:N

S = S + C(k+1)*y^k;

C(k+2) = C(k+1)+C(k); end

C([1 length(C)]) = [];

figure, plot(1:10,C(1:10),'r'), hold on, plot(11:length(C), C(11:length(C)),'b*')

Scrivere la funzione Matlab esercizio3.m che: • legga in input una matrice A; • effettui, se possibile, un opportuno scambio di righe per cui il

numero di elementi non nulli su ogni riga di A cresca con l’indice di riga;

• restituisca in output la matrice A, eventualmente modificata in accordo al punto precedente;

• determini il numero di elementi nulli di A e restituisca in output il coefficiente C dato dal rapporto tra il numero di elementi nulli di A e la dimensione di A definito;

• Se C =0 stampi il messaggio ‘La matrice non è sparsa’ Richiamare la funzione scritta sulle seguenti matrici e salvare tutti gli input e gli output nel file risultaties3.mat.

Esercizio 3

14

function [A,C] = esercizio3(A) % function [A,C] = esercizio3(A) % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. [m,n] = size(A); C = length(find(A(:)==0))/(m*n); if C==0 disp('la matrice non e'' sparsa'); else for i=1:m nnz(i) = length(find(A(i,:))); end [S,P] = sort(nnz); A=A(P,:); end

Esercizio 4 Scrivere la funzione Matlab esercizio4.m che: • legga in input le funzioni f, df, df2, con df e df2 rispettivamente

derivata prima e seconda di f, due numeri reali a e b (con a < b) e due interi positivi M e nmax, con nmax non superiore a 1000. Se M non viene dato in input, assegni ad M il valore 5; se nmax non viene dato in input, assegni ad nmax il valore massimo consentito;

• determini un intervallo di separazione I = [a1, b1] [a, b] della radice di f con ampiezza non superiore a 0.25;

• approssimi la radice di f in I con il metodo delle tangenti, scegliendo come approssimazione iniziale x0 l'estremo di Fourier e arrestando il procedimento quando la approssimazione prodotta xn ha M decimali esatti o non appena il numero di iterazioni eseguite diventa pari a nmax. Se l'estremo di Fourier non esiste, ponga x0 uguale all'estremo positivo più piccolo di I;

• restituisca in output xn, la variabile fxn contenente il valore di f in xn e l'approssimazione iniziale x0.

Utilizzare la funzione scritta per approssimare la radice della equazione non lineare x3-cos(x)+ 1/3 = 0 nell'intervallo [0.25, 1] e salvare sia gli input che gli output nel file risultaties4.mat.

function [xn,fxn,x0] = esercizio4(f,df,df2,a,b,varargin) % function [xn,fxn,x0] = esercizio4(f,df,df2,varargin) % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. % ………………………………………………………………………………………………………. switch nargin case 5 M = 5; nmax = 1000; case 6 M = varargin{1}; nmax = 1000; case 7 M = varargin{1}; nmax = varargin{2}; otherwise error('controllare il numero delle variabili di input!!!') end

if (M~=round(abs(M))) error('M deve essere intero positivo!!!') end if (nmax~=round(abs(nmax))) | nmax>1000 error('nmax deve essere intero positivo e non superiore a 1000!!!') end if a>=b error('a deve essere strettamente minore di b!!!') end x = a:.2:b; y=f(x); pos = find(abs(diff(sign(y)))); if length(pos)~=1 error('la radice non e'' unica') else a1 = x(pos); b1 = x(pos+1); end

if f(a1)*df2(a1)>0 x0 = a1; elseif f(b1)*df2(b1)>0 x0 = b1; else if sum(sign([a1 b1]))<=0 error('non e'' possibile fissare la condizione iniziale') elseif a1>0 x0 = a1; else x0 = b1; end end eps = 0.5*10^(-M); err = eps+1;, iter = 0; while (err>eps) & (iter<nmax) xn = x0-f(x0)/df(x0); iter = iter + 1; err = abs(xn-x0); x0=xn; end fxn = f(xn);

Switch

Switch-case-Otherwise:

switch nome_variabile (o espressione)

case valore1

1° blocco di istruzioni

case valore2

2° blocco di istruzioni

……….

otherwise

ultimo blocco di istruzioni

end

Se il valore di nome_variabile è valore1 viene eseguito il 1° blocco di

istruzioni;

Se è valore2, il 2° blocco di istruzioni, e così via,

altrimenti esegue l’ultimo blocco di istruzioni

Switch

Esempio: a = rem(b,3); % resto della divisione per 3

switch a

case 0

c = b/3;

case 1

c = (b-a)/3;

otherwise

c = b-a/2;

end

function [c] = esercizio(y,z,m) % ............................................... % ............................................... % ............................................... C = zeros(m);, v = zeros(m,1);, s = zeros(2*m-1,1); for i = 1:length(y) temp = z(i); for j = 1:m

v(j) = v(j) + temp;, temp = temp*y(i); end temp = 1; for j = 1:2*m-1 s(j) = s(j) + temp;, temp = temp*y(i); end end for i = 1:m for j = 1:m C(i,j) = s(i+j-1); end end c = C\v;

Esercizio 5 Scrivere lo help della funzione descrivendo

anche tutte le variabili di input e di output.

Commentare le istruzioni.

Esercizio 6 Scrivere lo help della seguente funzione; commentare e completare le istruzioni. function [V,M] = fun_esame(t,g) %............................................... %............................................... %............................................... %............................................... %............................................... J = length(g); D = (t(end)-......................); if ............. V = (g(1)+2*sum(g(2:2:J-1))+2*sum(g(2:J-1))+g(end))*D/3; M = 0; else V = D*sum(g) - 0.5*D*(g(1)+ .....................); M = 1; end

function [In,Nn,E] = my_fun(x1,x2,fun,varargin) %................................................................ %................................................................ %................................................................ if nargin<3 error('input non sufficienti') elseif ............................... tl=0.5*10^-5; elseif nargin == 4 tl = ............................. else error('troppe variabili di input') end D = (x2-x1)/2;, I0 =D/3*(fun(x1)+4*fun(x1+D)+fun(x2));, E = tl+1; while E >= tl

D = D/2; Nn = ............................ xn = linspace(x1,x2,Nn);, fn = fun(xn); In = D/3*(fn(1)+4*sum(fn(2:2:Nn-1))+2*sum(fn(3:2:Nn-2))+fn(Nn)); E = abs(In-I0)/15; I0 = ...........................

end

Esercizio 7 Completare e

commentare le istruzioni della

funzione e scriverne lo help

Scrivere la funzione Matlab esercizio8.m che • riceva in input tre numeri a, b, eps e la funzione f; se eps non

viene data in input, assegni il valore 0.510-5

• approssimi l'integrale con la formula di Cavalieri-

Simpson; • approssimi I(f) con la formula delle parabole generalizzata

dimezzando ogni volta il passo h e arrestandosi quando il modulo del resto Rh(f) diventa minore di eps;

• restituisca in output l' approssimazione finale dell' integrale, il numero totale di nodi utilizzati per calcolarla e la stima di Rh(f) ottenuta tramite il criterio di Runge.

Utilizzare la funzione realizzata per approssimare con 6 decimali esatti e salvare gli output nel file risultaties8.mat.

Esercizio 8

Esercizio 9

function [F,varargout] = fun_bruni(t,M,varargin) % [F] = fun_bruni(t,M); % descrizione…… % % INPUT % inserire descrizione delle variabili di input…. % % OUTPUT % inserire descrizione delle variabili di output…. switch nargin case 2 x=0; m=5; case 3 x=varargin{1}; m=5; if x<-1 | x>1 error('x deve appartenere all''intervallo [-1,1]') end

case 4 x=varargin{1}; m=varargin{2}; if (x<-1 | x>1) | length(x)>1 error('x deve essere uno scalare e appartenere all''intervallo [-1,1]') end if m<=0 | m~=round(m) error('m deve essere un intero positivo') end otherwise error('Controllare il numero delle variabili di input') end if M<=0 | M~=round(M) error('M deve essere un intero positivo') end

eps = 0.5*10^(-m); g = (1-t*x)/(1-2*x*t+t^2); F=[]; T(1) = 1; f(1) = T(1); if abs(f(1)-g) < eps F = [F f(1)]; end T(2) = x; f(2) = sum(T.*t.^(0:1)); if (abs(f(2)-g) < eps) & (length(F)<M) F = [F f(2)]; end n=3; while (length(F)<M) & n<100 T(n) = 2*x*T(n-1)-T(n-2); f(n) = sum(T.*t.^(0:n-1)); if (abs(f(n)-g) < eps) F = [F f(n)] end n=n+1; end

ind = find(T>=0); ind2 = find(T<0); figure title('grafico esercizio') subplot(1,2,1), plot(F,'r--') ylabel('vettore F') subplot(1,2,2), plot(ind-1, T(ind),'k*',ind2-1,T(ind2),'ro') xlabel('n') ylabel('Tn') if nargout ==2 V=ind-1; varargout{1}=V; end

Esercizio 10 • Scrivere la funzione lagrange.m che, dato un vettore di nodi

xnodi, un numero intero i e uno scalare xeval, valuti lo i-esimo polinomio della base di Lagrange nel punto xeval.

• Scrivere la funzione interp_lagrange.m che, dati due vettori, xnodi e fnodi aventi la stessa lunghezza, e uno scalare xeval, valuti il polinomio interpolatore nella forma di Lagrange costruito sui dati xnodi e fnodi nel punto xeval. La funzione deve richiamare opportunamente la funzione lagrange.m

Le funzioni devono contenere opportuni controlli sulle variabili di input. Usare la funzione interp_lagrange.m per costruire il polinomio che interpola i dati riportati nella seguente tabella e valutarlo nei punti xeval1 = 0.2, xeval2 = 0.6 e xeval3 = 1

function [yeval] = lagrange(xnodi, i, xeval) % function [yeval] = lagrange(xnodi,i,xeval) % La function valuta l'i-esimo polinomio della base di Lagrange costruito sul vettore di nodi % xnodi nel punto xeval % INPUT: % xnodi = vettore dei nodi % i = numero intero % xeval = scalare % % OUTPUT: % yeval = valore dell'i-esimo polinomio della base di Lagrange in xeval %numero di nodi n = length(xnodi); %controllo sull'input i if i~=floor(i) error('Il secondo input deve essere un numero intero.') end %costruzione e valutazione in xeval dell'i-esimo polinomio della base di Lagrange yeval = 1; for j=1:n if (i~=j) yeval=yeval*(xeval-xnodi(j))/(xnodi(i)-xnodi(j)); end end

function [y] = interp_lagrange(xnodi, fnodi, xeval ) %function [y] = interp_lagrange(xnodi, fnodi, xeval ) %La function valuta, nel punto xeval, il polinomio interpolatore nella forma di Lagrange %costruito sui dati xnodi e fnodi % INPUT: % xnodi = vettore dei nodi % fnodi = valori del polinomio corrispondenti ai nodi % xeval = scalare % % OUTPUT: % y = valore del polinomio interpolatore in xeval %numero di nodi n = length(xnodi); if length(fnodi)~=n error('I vettori xnodi ed fnodi devono avere la stessa lunghezza!') end % costruzione e valutazione in xeval del polinomio interpolatore nella forma di Lagrange, % richiamando la function lagrange y=0; for i=1:n y=y+fnodi(i)*lagrange(xnodi,i,xeval); end

%script lagrange_example.m xnodi = [0.0 0.4 0.8 1.2]; fnodi = [0.00000 0.42839 0.74210 0.91031]; xeval1 = 0.2; xeval2 = 0.6; xeval3 = 1; yeval(1)=interp_lagrange(xnodi, fnodi, xeval1 ); yeval(2)=interp_lagrange(xnodi, fnodi, xeval2 ); yeval(3)=interp_lagrange(xnodi, fnodi, xeval3 ); figure, plot(xnodi,fnodi,'ob') hold on plot([xeval1 xeval2 xeval3],yeval,'*r') xlabel('x') ylabel('y') legend('valori nei nodi', 'valori stimati') hold off

Dal command window >> lagrange_example

35

Scrivere la function Matlab lebesgue.m che: 1. riceva in input i vettori xnodi e xpunti; 2. costruisca il vettore fleb contenente i valori della funzione di Lebesgue

calcolata in xpunti; 3. se xpunti contiene più di un elemento, grafichi la funzione di Lebesgue

calcolata in xpunti etichettando opportunamente gli assi; 4. restituisca in output il massimo valore di fleb e il punto p_Mfleb in cui è

realizzato. Utilizzare la function per il seguente vettore di nodi [1.0; 1.5; 2.0; 2.5; 3.0] e salvare gli output nel file risultati_esercizio.mat

Esercizio 11

36

function [Mfleb,p_Mfleb] = lebesgue(xnodi,xpunti) % [xpunti,Mfleb,p_Mfleb] = lebesgue(xnodi,xpunti) valuta e grafica la funzione di % Lebesgue costruita sui nodi contenuti in xnodi in corrispondenza dei punti in % xpunti, ne calcola il massimo e il punto in cui quest'ultimo è realizzato. % % INPUT % - xnodi = vettore dei nodi % - xpunti = vettore dei punti in cui valutare la funzione di Lebesgue % % OUTPUT % - Mfleb = valore massimo assunto dalla funzione di Lebesgue sul vettore di % punti in xpunti % - p_Mfleb = punto di xpunti in cui la funzione di Lebesgue assume il valore % massimo % controllo degli input if length(xnodi)<=1

error('xnodi deve essere un vettore di almeno due elementi!'); end NOTA: Modificare inserendo il seguente controllo su xpunti:

gli elementi di xpunti devono essere contenuti nell’intervallo [xnodi(1), xnodi(end)]

37

% Calcolo i polinomi di base di Lagrange e li valuto in ogni elemento di xpunti n_xnodi = length(xnodi); n_xpunti = length(xpunti); for i=1:n_xpunti % per ogni elemento di xpunti for n=1:n_xnodi % per ogni nodo in xnodi plag(i,n) = prod(xpunti(i)-xnodi(1:n-1))*prod(xpunti(i)-xnodi(n+1:n_xnodi))/… (prod(xnodi(n)-xnodi(1:n-1))*prod(xnodi(n)-xnodi(n+1:n_xnodi))); end end

NOTA: La variabile plag è una matrice con n_xpunti righe e n_xnodi colonne Ogni elemento in posizione (i,n) contiene il valore che il polinomio della base di Lagrange relativo al nodo n-simo assume in corrispondenza dello i-simo elemento di xpunti

38

% calcolo il valore della funzione di Lebesgue in ogni elementi di xpunti fleb = sum(abs(plag')); % calcolo il massimo di fleb [Mfleb,p_Mfleb] = max(fleb); % grafico della funzione di Lebesgue if n_xpunti>1 figure, plot(xpunti,fleb) xlabel('xpunti') ylabel('funzione di Lebesgue') end

39

Dal Command Window >> xnodi=[1.0, 1.5, 2.0, 2.5, 3.0]; >> [Mfleb,p_Mfleb] = lebesgue(xnodi,2.2); >> save risultati_esercizio.mat Mfleb p_Mfleb

>> [xpunti,Mfleb,p_Mfleb] = lebesgue(xnodi(1),2.2); ??? Error using ==> lebesgue at 24 xnodi deve essere un vettore di almeno due elementi!

Non produce grafici

Alcune funzioni utili

43

Esempio

44

Esercizio 12

%script esercizio_molla.m x=[2 3 4 5 6 8 10]; y=[7 8.3 9.4 11.3 12.3 14.4 15.9]; a = polyfit(x,y,1); disp(a) yn = polyval(a,x); plot(x,y,'-ob'), hold on, plot(x,yn,'r'); xlabel('x') ylabel('y') legend('valori nei nodi','P1(x)','Location','northwest') hold off

Dal command window >> esercizio_molla 1.1383 5.0491

Costante elastica della molla

Esercizio 13 (compito d’esame gennaio2014)

Oss: • La funzione F(x,y) dipende anche dall’indice naturale n • I coefficienti ak dipendono anche dall’indice naturale n • Per il vettore A, poiché non specificato, si considera l’ultima sequenza calcolata

Si consideri la seguente funzioni di due variabili Scrivere la funzione Matlab fun_cognome.m che riceva in input un intero strettamente positivo N non superiore a 100, un numero reale x e una funzione f (anonymous function) dipendente da x. Se f non viene data in input, ponga f(x) = 1 + x. La funzione valuti F in corrispondenza della coppia di valori (x; f(x)) per ogni numero naturale n : 0 n N; n = 3h; h N. La funzione restituisca in output la matrice P così definita: i) la prima riga contiene i valori di n che hanno prodotto valori di F positivi; ii) la seconda riga contiene i corrispondenti valori positivi di F; iii) la terza riga contiene i valori di n che hanno prodotto valori di F negativi; iv) la quarta riga contiene i corrispondenti valori negativi di F. La funzione grafichi, e restituisca in output se richiesto, il vettore A contenente al piu (N+1)/3 elementi della successione ak usando un pallino rosso per gli elementi di indice pari e un asterisco blu per quelli dispari. Il grafico deve essere completo di titolo, etichette per gli assi e legenda.

function [P, varargout] = fun_bruni(N, x, varargin) % [P, varargout] = fun_bruni(N, x, varargin) % La funzione fun_bruni valuta la funzione F(x,y) che è la somma di n+1 termini il cui (k+1)- % simo addendo è a(k) = (n-k+1)/k a(k-1) se k>=1, e a(0) = 1 se k=0. % La funzione calcola F(x,y) per tutti i valori di 0<=n<=N e multipli di 3 e con y=f(x). % % INPUT % N = intero positivo <=100 % x = ascissa del punto in cui valutare la funzione % f = funzione che definisce l'ordinata del punto in cui valutare la funzione come y=f(x). % Se f non viene data in input, f=1+x % % OUTPUT % inserire descrizione degli output % controllo sul numero degli input if nargin == 2 f = @(x)(1+x); elseif nargin == 3 f = varargin{1}; else error('numero di input non corretto') end

% controllo input if N <= 0 || N > 100 || N ~= round(N) error('N deve essere un intero strettamente positivo non superiore a 100') end % successione y = f(x); for n=0:3:N a(1) = 1; F(n/3+1)=a(1)*y^n; for k=2:n+1 a(k) = (n-k+2)/(k-1)*a(k-1); F(n/3+1) = F(n/3+1) + a(k)*x^(k-1)*y^(n-k+1); end end pos = find(F>=0); neg = find(F<0); P = zeros(4,max(length(pos), length(neg))); P(1,1:length(pos)) = (pos-1)*3; P(2,1:length(pos)) = F(pos); P(3,1:length(neg)) = (neg-1)*3; P(4,1:length(neg)) = F(neg);

A = a(1:round((N+1)/3)); figure, hold on plot(1:2:length(A), A(1:2:length(A)), '*b') plot(2:2:length(A), A(2:2:length(A)), '*r') if nargout==2 varargout{1} = A; end

Esercizio 14 (compito d’esame gennaio2014)

function [K, S, D, V, varargout] = bruni(A, b, varargin) %INSERIRE HELP % controllo sugli input if nargin < 2 || nargin > 3 error('il numero degli input non e'' corretto') elseif nargin == 2 C = {'J', 'GS', 'SOR', 'MEG', 'Th'}; elseif nargin == 3 C = varargin{1}; end % controllo sulle dim di A [m,n] = size(A); if m ~= n error('A deve essere quadrata') else if m < 3 error('la dim di A deve essere non inferiore a 3x3') end end

K = false; S = false; D = false; if cond(A) < 7 K = true; end % S vale 1 se il numero degli elementi nulli di A è superiore alla metà dei % suoi elementi n_zeri = length(find(A==0)); if n_zeri > (n^2)/2 S = true; end % D vale 1 se esiste una permutazione di righe per cui la matrice del % sistema diventa diagonale [r, c] = find(A~=0); if length(r) <= n if (length(unique(c)) == length(unique(r))) D = true; if nargout == 5 [temp p] = sort(r); varargout{1} = p; end end end

fprintf('%s \n', 'quali metodi possono essere usati per risolvere Ax=b?') riga1 = zeros(1,length(C)); riga2 = zeros(1, length(C)); eps = 10^(-8); MD = diag(diag(A));, L = tril(A,-1);, U = triu(A,1); for i=1:length(C) switch C{i} case {'J'} Minv = inv(MD);, CJ = -Minv*(L + U); rho = max(abs(eig(CJ))); if rho < 1 fprintf('%s \n', 'il metodo di Jacobi converge'); riga1(i) = 1; QJ = Minv*b; iter = 0;, X0 = ones(n,1);, err = 100; while err > eps X = CJ*X0 + QJ; err = norm(X-X0, inf); X0 = X; iter = iter + 1; end riga2(i) = iter; end

case {'GS'} Minv = inv(tril(A)); CGS = -Minv*(triu(A,1)); rho = max(abs(eig(CGS))); if rho < 1 fprintf('%s \n', 'il metodo di Gauss Seidel converge'); riga1(i) = 1; QGS = Minv*b; iter = 0; X0 = ones(n,1); err = 100; while err > eps X = CGS*X0 + QGS; err = norm(X-X0, inf); X0 = X; iter = iter + 1; end riga2(i) = iter; end

case {'SOR'} T = (A == A'); if all(T(:)) && all(eig(A)) fprintf('%s \n', 'il metodo di SOR converge'); riga1(i) = 1; omega = 0.5; Minv = inv(L + MD/omega); C_sor = -Minv*(U-(1-omega)*MD/omega); Q_sor = Minv*b; iter = 0; X0 = ones(n,1); err = 100; while err > eps X = C_sor*X0 + Q_sor; err = norm(X-X0, inf); X0 = X; iter = iter + 1; end riga2(i) = iter; end

case {'MEG'} if det(A) ~= 0 fprintf('%s \n', 'e'' possibile usare il metodo di eliminazione di Gauss'); riga1(i) = 1; riga2(i) = n^2/2 + n^3/3; end case {'Th'} B = triu(A,-1) - triu(A,2); T = (B == A); if det(A) ~= 0 && all(T(:)) fprintf('%s \n', 'e'' possibile usare l'' algoritmo di Thomas'); riga1(i) = 1; riga2(i) = 5*n-4; end otherwise error('C contiene metodi sconosciuti') end end %% costruisco la matrice V V(1,:) = riga1; V(2,:) = riga2;

Esercizio 15 (compito prova preliminare dicembre 2015)

function [S] = fun_bruni(x,varargin) %SOLUZIONE ES 1 PROVA PRELIMINARE 18 DICEMBRE 2015 %... %..._ switch nargin case 1 N = 100; M = 5; case 2 N = varargin{1}; M = 5; case 3 N = varargin{1}; M = varargin{2}; otherwise error('controllare input') end if (M~=round(M)) | M<0 error('M deve essere intero positivo') end if (N~=round(N)) | N<0 | N>100 error('N deve essere intero positivo non superiore a 100') end

eps = 0.5*10^(-M); F = exp(x^2); S=[]; ind = []; H(1)= 1; G(1)= H(1); if abs(G(1)-F) <= eps S = [S [G(1); 1]]; end m=1; while ((size(S,2)<N) && (m<=100)) if m==1 H(2) = 2*x; else H(m+1) = 2*x*H(m)-2*m*H(m-1); end G(m+1) = G(m) + x^m/factorial(m)*H(m+1); if abs(G(m+1)-F) <= eps S = [S [G(m+1); m]]; else ind = [ind m]; end m = m+1; end

Numero massimo di iterazioni consentite --- stabilito in fase di programmazione

figure plot(ind,G(ind+1),'ro') hold on, if ~isempty(S) plot(S(2,:),S(1,:),'go') end xlabel('m') ylabel('Gm') title('Successione') legend('elementi di G non in S','elementi di G in S')

Esercizio 16 (compito prova preliminare dicembre 2015)

function [A,P,varargout] = bruni(A,varargin) %SOLUZIONE ES 2 PROVA PRELIMINARE 18 DICEMBRE 2015 %... %..._ if nargin==0 error('pochi input') elseif nargin == 1 if (length(size(A))==2) & size(A,1)>1 & size(A,2)>1 B = zeros(size(A,2),1); else error('controllare la dimensione di A') end elseif nargin == 2 if (length(size(A))~=2) | size(A,1)==1 | size(A,2)==1 error('controllare la dimensione di A') end B = varargin{1}; B = B(:); if length(B)~=size(A,1) error('controllare la dimensione di B') end else error('troppi input') end

[nr,nc] = size(A); if nr ~= nc error('A deve essere quadrata') end if det(A) == 0 error('la matrice e'' singolare') end if all(tril(A)==A) P=1:nr; else for i=1:nr ind = find(A(i,:)~=0); P(i) = ind(end); end end

if length(unique(P))~=nr P = []; A = []; x = []; disp('non esiste permutazione che rende A triangolare inferiore') else A(P,:) = A; X(1) = B(1)/A(1,1); for i=2:nr X(i) = (B(i)-sum(A(i,1:i-1).*X(1:i-1)))/(A(i,i)); end X=X(:); end if nargout==3 varargout{1}=X; elseif nargout<=1 error('pochi output') elseif nargout >3 error('troppi output') end

Esercizio 17

Scrivere la funzione Matlab diff_div.m che riceva in input i vettori xnodi e fnodi e restituisca in output la matrice mat_diff che contiene tavola delle differenze divise riferita ad una funzione i cui valori in corrispondenza dei nodi in xnodi sono gli elementi del vettore fnodi. Le colonne della matrice mat_diff si riferiscono alle differenze divise di uno stesso ordine. La prima colonna di mat_diff deve contenere le differenze divise di ordine 0 Richiamare opportunamente la funzione diff_div.m nella funzione pol_difdiv.m. Quest’ultima riceve in input i vettori xnodi e fnodi e il vettore di punti xeval in cui si vuole valutare il polinomio di Newton alle differenze divise costruito usando i dati in xnodi e fnodi; restituisce in output il vettore pn dei valori del polinomio interpolatore calcolato in xeval e la matrice tavdif contenente la tavola delle differenza divise calcolata usando la funzione diff_div.m

function [mat_diff] = diff_div(xnodi,fnodi) % mat_diff = diff_div(xnodi,fnodi) % Calcolo della tavola alle differenze divise a partire dai valori dei nodi % e dalla funzione nei nodi (memorizzati nella prima colonna della % matrice mat_diffdiv). % % Input: % xnodi: vettore delle coordinate dei nodi % fnodi: valori della funzione nei nodi % Output: % mat_diff: tavola delle differenze divise (per colonne) nella prima % colonna ci sono i valori fnodi nnodi=length(xnodi); % ordine massimo delle differenze divise kord = nnodi-1;

% inizializzazione della tavola alle differenze divise la prima colonna % contiene i valori della funzione nei nodi mat_diff = zeros(nnodi,kord+1); mat_diff(:,1) = fnodi'; % calcolo della tabella for j = 1:kord, % indice per le colonne della tavola (ordine delle diff. divise) for i = 0:nnodi-1-j, % indice per le righe della tavola (nodi coinvolti) mat_diff(i+1,j+1) = (mat_diff(i+1,j)-mat_diff(i+2,j))/(xnodi(i+1)-xnodi(i+j+1)); end

end

function[pn,tavdif] = pol_difdiv(xnodi,fnodi,xeval) % Polinomio interpolatore alle differenze divise: % [pn,tavdif] = pol_difdiv(xnodi,fnodi,xeval) % Input: % xnodi: vettore delle coordinate dei nodi % fnodi: valori nei nodi della funzione da interpolare % xeval: vettore delle coordinate dei punti in cui calcolare %il polinomio interpolatore % Output: % pn: vettore dei valori del polinomio interpolatore calcolato in xeval % tavdif: tavola delle differenza divise calcolata usando la funzione % diff_div.m % controllo degli input if length(xnodi)~=length(fnodi) error('il vettore dei nodi ha dimensione diversa dal vettore dei valori della funzione!!!')

end

% Calcolo della tavola alle differenze divise tavdif = diff_div(xnodi,fnodi); % Calcolo del polinomio interpolatore di grado nnodi-1 nnodi = length(fnodi); pn = tavdif(1,1); % primo elemento della tavola (f[x0]) pnod = 1; %inizializzazione dei polinomi nodali for i=1:nnodi-1 pnod = pnod.*(xeval-xnodi(i)); pn = pn + pnod*tavdif(1,i+1); end

Dal Command Window: >> xnodi = [1.90 2.15 2.20 2.25]; >> fnodi = [-0.10536 0.13976 0.18232 0.22314]; >> [pn,tavdif] = pol_difdiv(xnodi(1:3),fnodi(1:3),2); >> pn pn = -8.479999999996919e-004 Calcolando il polinomio interpolatore usando tutti i nodi della tabella e valutandolo nel punto x = 2 si ha p3(2) =-0.00014 >> xnodi = [1.90 2.15 2.20 2.25]; >> fnodi = [-0.10536 0.13976 0.18232 0.22314]; >> [pn,tavdif] = pol_difdiv(xnodi,fnodi,2); >> pn pn = -1.371428571418221e-004

Esercizio 18 Modificare la funzione pol_difdiv.m in modo che restituisca come ulteriore variabile di output una stima dell'errore di troncamento qualora venga richiesto. In tal caso, la funzione deve prevedere come ulteriore variabile di input un nodo aggiuntivo e il valore della funzione in questo punto.