ESERCIZI^{∗}

A. SOMMARIVA^{†}

1. Approssimazione di funzioni. Se si intende approssimare f(x) = x^{4} nell’intervallo [0, 1] con una retta g_{2}(x) = a_{1} + b_{1}x cosicchè sia minima la quantità

∫_{0}^{1}(f(x) − g_{2}(x))^{2} dx

si ottiene quale miglior approssimante

g_{2}^{∗}(x) = −\frac{1}{5} + \frac{4}{5}x.

Se invece si intende minimizzare con una retta g_{∞}(x) = a_{2} + b_{2}x la quantità

max_{x∈[0,1]} |x^{4} − g_{∞}(x)|

si ottiene

g_{∞}^{∗}(x) = \frac{−3 · 2^{1/3}}{16} + x.

E’ questa una contraddizione? Perchè ?

1. Usando Matlab/Octave e il comando quad (aiutarsi con l’help!) calcolare

∫_{0}^{1}(f(x) − g_{2}(x))^{2} dx,

∫_{0}^{1}(f(x) − g_{∞}^{∗}(x))^{2} dx.

2. Definiti 201 punti equispaziati x_{k} = k/200 ove k = 0, …, 200 nell’intervallo [0, 1] e aiutandosi con l’help del comando linspace, valutare

max_{x_{k}} |x_{k}^{4} − g_{2}(x_{k})|

max_{x_{k}} |x_{k}^{4} − g_{∞}(x_{k})|.

3. Disegnare il grafico delle funzioni f, g_{2}, g_{∞}. Inoltre plottare in scala semilogaritmica |f − g_{2}| e |f − g_{∞}|.

∗Ultima revisione: 4 marzo 2010 †Dipartimento di Matematica Pura ed Applicata, Universitá degli Studi di Padova, via Trieste 63, 35121 Padova,

Italia (alvise@euler.math.unipd.it).

Cosa attestano questi risultati numerici?

Suggerimenti:

1. Per l’utilizzo di quad in Octave può essere d’aiuto

http://www.math.uic.edu/h̃anson/Octave/OctaveIntegralEG.html

2. Per l’utilizzo di quad in Matlab può essere d’aiuto

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/quad.html

2. Interpolazione trigonometrica e Fast Fourier Transform.

1. Utilizzando il software fornito calcolare i coefficienti di Fourier per n = 8, relativi alla funzione f(x) = cos (2x). Giustificare tale risultato. Suggerimento: utilizzare la formula di Eulero e riscrivere f come combinazione lineare di due esponenziali in campo complesso.

2. Giustificare perchè esempio FFT calcola i coefficienti di Fourier. Dire come si poteva scrivere più semplicemente l’ultima riga di esempio FFT

   c_fft=fft(ft,n)*h/(2*pi);

3. Quadratura numerica. Si calcolino per N = 10, 20 con la formula composta dei trapezi, di Cavalieri-Simpson e un’appropriata formula gaussiana i seguenti integrali

∫_{−1}^{1} x^{20} dx = 2/21 ≈ 0.095238095238096801  (3.1)
∫_{−1}^{1} e^{x} dx = e − e^{−1} ≈ 2.3504023872876032  (3.2)
∫_{−1}^{1} e^{−x^{2}} dx = erf(1)·\sqrt{π} ≈ 1.4936482656248538  (3.3)
∫_{−1}^{1} 1/(1 + 16x^{2}) dx = 1/2 · atan(4) ≈ 0.66290883183401628  (3.4)
∫_{−1}^{1} e^{−1/x^{2}} dx ≈ 0.17814771178156086  (3.5)
∫_{−1}^{1} |x|^{3} dx = 1/2  (3.6)
∫_{0}^{1} \sqrt{x} dx = 2/3  (3.7)
∫_{−1}^{1} e^{x}\sqrt{1 − x} dx ≈ 1.7791436546919095  (3.8)

1. Quali delle due formule ha errori relativi inferiori?

2. Quali funzioni risultano più difficili da integrare numericamente?

3. La scelta N = 10, 20 nei singoli codici a quanti nodi corrisponde?

4. Suggerimento: ricordarsi che il peso w(x) ≡ 1 corrisponde ad un’opportuna scelta del peso di Jacobi.

4. Autovalori di una matrice. Data la matrice di Hilbert di ordine 5, ottenibile in Matlab col comando hilb(5) si calcolino col metodo delle potenze i suoi minimi e massimi autovalori in modulo. Da questi si calcoli il condizionamento della matrice in norma 2 e lo si confronti con cond(hilb(5),2). Eseguire lo stesso esercizio utilizzando il metodo QR.

5. Metodi iterativi per l’algebra lineare.

1. Modificare in demo algebra lineare
maxit=200; tol=10^(-8);
siz=5;
A = makefish(siz);      % MATRICE DI POISSON.

la dimensione della matrice di Poisson, variando il valore assunto da siz (senza che questo superi il valore 10). Se opportuno alzare il valore della variabile maxit, nel caso l’uscita del metodo di Jacobi sia per un eccesso di iterazioni e non per un’avvenuta convergenza (a tal proposito si osservi il valore assunto dalla variabile flag che è pari a uno quando il metodo esce per un numero eccessivo di iterazioni). Come si comportano i metodi? E’ il numero di iterazioni di Jacobi pari a circa il doppio di quelle eseguite da Gauss-Seidel (a parità di tolleranza)?

2. nel set di matrici offerte dal comando Matlab gallery si cerchi una matrice di piccole dimensioni (avente cioè al massimo grandezza 50 × 50 per motivi computazionali) che sia simmetrica e definita positiva. Per verificare ciò facciamo un esempio che spiega come usare la gallery di Matlab (si ricorda che non esiste questo comando in Octave).

>>A=gallery('lotkin',5)
A =
     1.0000     1.0000     1.0000     1.0000     1.0000
     0.5000     0.3333     0.2500     0.2000     0.1667
     0.3333     0.2500     0.2000     0.1667     0.1429
     0.2500     0.2000     0.1667     0.1429     0.1250
     0.2000     0.1667     0.1429     0.1250     0.1111
>> % COMMENTO: LA MATRICE NON E' SIMMETRICA.
>>
>>% VERIFICA. SE IL RISULTATO ASSUNTO DA norm(A'-A)  E' ZERO
>>% ALLORA LA MATRICE "A" E' SiMMETRICA.
>> norm(A'-A)
ans =
     1.3773
>>
>>% VERIFICA SE LA MATRICE SIA DEFINITA POSITIVA,
>>% CIOE' AUTOVALORI TUTTI POSITIVI (STRETTAMENTE).
>>  min(eig(A))
ans =

-0.2138
>> % COMMENTO: UN AUTOVALORE E' NEGATIVO.
>> % LA MATRICE NON E' DEFINITA POSITIVA.

Una volta trovata una tale matrice A simmetrica e definita positiva, si modifichi il file demo algebra lineare così da confrontare su un test il metodo di Jacobi, Gauss-Seidel, SOR ottimale e il gradiente coniugato. Converge il metodo di Jacobi?

Nota 5.1. Per vedere alcuni comandi di base aiutiamoci con delle matrici predefinite in Matlab/Octave. Digitiamo nella shell di Matlab/Octave >> help elmat. In Matlab 6.5 abbiamo

>> help elmat
Elementary matrices and matrix manipulation.
Elementary matrices.
zeros        - Zeros array.
ones         - Ones array.
eye          - Identity matrix.
repmat       - Replicate and tile array.
rand         - Uniformly distributed random numbers.
randn        - Normally distributed random numbers.
linspace     - Linearly spaced vector.
logspace     - Logarithmically spaced vector.
freqspace   - Frequency spacing for frequency response.
meshgrid     - X and Y arrays for 3-D plots.
:             - Regularly spaced vector and index into matrix.
...
Specialized matrices.
compan       - Companion matrix.
gallery      - Higham test matrices.
hadamard     - Hadamard matrix.
hankel       - Hankel matrix.
hilb         - Hilbert matrix.
invhilb      - Inverse Hilbert matrix.
magic        - Magic square.
pascal       - Pascal matrix.
rosser       - Classic symmetric eigenvalue test problem.
toeplitz     - Toeplitz matrix.
vander       - Vandermonde matrix.
wilkinson   - Wilkinson's eigenvalue test matrix.

Questo ci dice che Matlab ha predefinito un set di matrici di particolare interesse. Se possibile si suggerisce di provare i metodi che andremo ad introdurre con una matrice facente parte della gallery di Matlab. Ciò non appare possibile nelle recenti releases di Octave, come GNU Octave 2.1.73. Per gli utenti Octave, ci si aiuti con

http://users.powernet.co.uk/kienzle/octave/testmatrix-pak.tar.gz

Da Matlab 6.5

>> help gallery
GALLERY Higham test matrices.
[out1,out2,...] = GALLERY(matname, param1, param2, ...)
takes matname, a string that is the name of a matrix family, and
the family's input parameters. See the listing below for available
matrix families. Most of the functions take an input argument
that specifies the order of the matrix, and unless otherwise
stated, return a single output.
For additional information, type "help private/matname", where matname
is the name of the matrix family.
cauchy   Cauchy matrix.
chebspec Chebyshev spectral differentiation matrix.
chebvand Vandermonde-like matrix for the Chebyshev polynomials.
chow      Chow matrix -- a singular Toeplitz lower Hessenberg matrix.
circul   Circulant matrix.
...
poisson  Block tridiagonal matrix from Poisson's equation (sparse).
prolate  Prolate matrix -- symmetric, ill-conditioned Toeplitz matrix.
randcolu Random matrix with normalized cols and specified singular
values.
randcorr Random correlation matrix with specified eigenvalues.
randhess Random, orthogonal upper Hessenberg matrix.
rando     Random matrix with elements -1, 0 or 1.
randsvd  Random matrix with pre-assigned singular values and specified
          bandwidth.
redheff  Matrix of 0s and 1s of Redheffer.
riemann  Matrix associated with the Riemann hypothesis.
ris       Ris matrix -- a symmetric Hankel matrix.
smoke     Smoke matrix -- complex, with a "smoke ring" pseudospectrum.
toeppd   Symmetric positive definite Toeplitz matrix.
toeppen  Pentadiagonal Toeplitz matrix (sparse).
tridiag  Tridiagonal matrix (sparse).
triw      Upper triangular matrix discussed by Wilkinson and others.
wathen   Wathen matrix -- a finite element matrix (sparse, random
entries).
wilk      Various specific matrices devised/discussed by Wilkinson.
          (Two output arguments)
GALLERY(3) is a badly conditioned 3-by-3 matrix.
GALLERY(5) is an interesting eigenvalue problem.  Try to find
its EXACT eigenvalues and eigenvectors.
See also MAGIC, HILB, INVHILB, HADAMARD, WILKINSON, ROSSER, VANDER.

6. Equazioni differenziali ordinarie.

1. Approssimare per λ = −100 il valore assunto dalla soluzione \overline{y} del problema di Cauchy

\{ \begin{matrix}y′(x) = λy(x), x ≥ 0 \\ y(0) = 1 \end{matrix}  (6.1)

nel punto x = 0.2. A tal proposito si utilizzano i metodi di Eulero esplicito e Eulero implicito (nella versione con il metodo delle secanti) con passi h = 0.1, h = 0.05, h = 0.02, h = 0.01, h = 0.001. Al variare di h verificare quando |1 + hλ| < 1. Se eseguiti correttamente dovrebbero fornire quali risultati (cf. [?, p.396])

   1.0000e-01   8.1000e+01   8.2645e-03
   5.0000e-02   2.5600e+02   7.7160e-04
   2.0000e-02   1.0000e+00   1.6933e-05
   1.0000e-02   2.0612e-09   9.5161e-07
   1.0000e-03   1.3556e-09   3.2046e-09

in cui alla prima colonna si stabilisce il passo h, nella seconda si espongono i risultati ottenuti dal metodo di Eulero esplicito e nella terza quelli di Eulero implicito (nella versione con il metodo delle secanti).

2. Approssimare per λ = −500 il valore assunto dalla soluzione \overline{y} del problema di Cauchy

\{ \begin{matrix}y′(x) = λy(x), x ≥ 0 \\ y(0) = 1 \end{matrix}  (6.2)

nel punto x = 0.2. A tal proposito si utilizzano i metodi di Eulero esplicito e Eulero implicito (nella versione con il metodo delle secanti) con passi h = 0.1, h = 0.05, h = 0.02, h = 0.01, h = 0.001, h = 0.0005. Al variare di h verificare quando |1 + hλ| < 1. Il metodo di Eulero implicito per h piccoli può richiedere molto tempo di calcolo.