Teoria completa Numerical Methods in Engineering Sciences, prof.Sangalli

Numerical methods

Reminders on matrices A ∈ R^n×n:

• A symmetrica se A = A^T, autovalori di una simmetrica sono reali.
• Matric A definita positiva se x^TAx > 0 ∀x ∈ Rⁿ con x ≠ ρ
• A non singolare → ATA simmetrica e positiva definita

Natural norm: ........ ||A|| .......... ||Av||

Verifica → || Ay || ≤ || A || || y ||

|| A + B || ≤ || A || + || B ||

Solving linear systems: the problem → given b ∈ Rⁿ, c ∈ R^r x ∈ R^r solution of

Ax = b → sol → x = A^-1 b

Solution unica se e solo se A e' non singolare (o invertible)

Simplest systems are diagonal systems → Dx = b con D =

• d₁₁ 0
• 0 d₂₂

x_i = b_i/d_jj i= 1, 2, ..., n

n operations

Lower triangular matrices → L = { l1 l2... 0 }

System can be solved “forward”

• x₁ =
• x₂ = b₂ - l₂₁x₁
• x_n = b_n - ∑ (from j=1 to n-1) l_ij x_j

A special lower triangular matrix (L has only 1 on the diagonal) →

• x₁ = b₁
• x₂ = b₂ - ∑ l_ijx_j

Upper triangular matrices → U =

• u₁₁ u₁₂ u₂₁
• 0 u₂₂ u₂₁

Execution time: t ≈ # operations * ... flops

Numerical methods

There are two classes of methods for solving linear systems:

• Direct: They give the exact solution in a finite number of operations
• Iterative: Starting from an initial guess x⁽⁰⁾, they construct a sequence {x^(k)} such that x = lim_{k → ∞} x^(k)

The most important direct method is GEM (Gaussian Elimination Method)

The original system is transformed into an equivalent system having the same solution

Ux = b con U matrice triangolare superiore

Il costo totale della costruzione c_c ≈ 1/3 n³ + 2n² n ≈ 2/3 n³

A^(t+1) x = A^(t) b^{(t) = b}

Eliminazione la prima colonna:

Continuo per n step, ottenendo A^(n,k)x_b⁽ⁿ⁾ con A triangolare superiore

GEM algorithm

Input A ∈ ℝ^n×m and b ∈ ℝⁿ

For k = 1, 2, ..., n for i = k+1, ..., n

End end

Define U = A⁽ⁿ⁾, b⁽² = b. Then solve Ux = b⁽²⁾ with back substitution

La condizione det(A) ≠ 0 no garantisce che la procedura di eliminazione sia un successo. Per evitare questo usare il metodo con pivot.

Primo step: Prima di eliminare la prima colonna, trovare il coefficiente della colonna con il pi. Alto valore assoluto e scambiare la riga dove si trova questo coefficiente con la prima riga.

Secondo step: Prima di eliminare le secondy colonne, trovare il coeff della colonna con pi alto valore assoluto e scambiarlo con la seconda riga. Il pattern della seconda riga.

Step j: Prima di eliminare le j-esima, il coefficiente che entra in quella riga j con pi alto valore assoluto e scambiarlo con riga j.

The pivoting procedure corresponds then to solve the system PAx = Pb

Permutation matrix

Usando il metodo di Gaussiano, risolvi:

• -5x₁ + 3x₂ + 4x₃ = 1
• -2x₁ - x₃ = 7
• 9x₃ = 39

Fatto Sol. x₃ = 13/3, x₂ = -13/3, x₁ = 2/15

Usando il metodo di Gaussiano, asolvi:

• x₁ - 3x₂ + 4x₃ = 12
• 2x₂ + x₃ = 8
• 5x₃ = 10

Fatto Sol. x₃ = 2; x₂ = -1; x₁ = 1

Parte II: Metodo Iterativi

Consideriamo un sistema lineare Ax = b. I metodi iterativi partono da una supposizione iniziale x⁽⁰⁾ e costruiscono una sequenza approssimata di soluzioni tale che

x = lim_k→∞ x^(k)

Splitting Methods

Se guida la matrice A viene divisa: A = M - N

x⁽ⁿ⁾ dopo risolvere Mx^(k+1) = Nx^(k) + b k = 1,2,...

Con i metodi iterativi approssimativi prima di trovare l'esatta soluzione, si vogliamo un basso costo computazionale. Invece, vogliamo che per ogni supposizione iniziale x⁽⁰⁾ si converga velocemente, q (x^k+1)

1. Attendere spetto unirezion di N divano vive di anessi uccesi ive metodi iterativi
2. Jacobi O

Mf si digit[sp] aaplicabile se a_ii ≠ 0 ∀i. A classica matrice, lo consciene rissolver un sistrama iterizione

Con i = 1, ..., n

end

SC = 5.3

PC = 9.07

GAUSS.m

function [x, res, j] = Gauss (A, b, x0, imax, tol)

n = size (b, 1);

x = zeros(n, 1);

for k = 1:n-1,

x(k+1:n) = b(k+1:n) - U((k+1):n, x(k+1:i n)) * x((k+1):n)) / U(k, k);

end

end

Consider the system Ax=b with A=gallery('poisson',n) and b=(1,...,1)^T

For n=10 solve the system using the Jacobi Gauss and CG (as implemented by the prog functional) methods compare the history of convergence of the three methods

JACOBISOLVER.m

n = 10

A = gallery('poisson',n);

b = ones(n 1);

x0 = zeros(n^2, 1);

tol = 1e-8;

maxit = 130

sc=((abs( eigs(A, 1) );

x = special focuss if to converge is guaranteed

c=condest (A); %condition number of A

x = A\b;

[~, res_J, its_J] = Jacobi(A, b, tol, maxit, x0);

[~, res_g, its_g] = Gauss (A, b, tol, maxit, x0);

[X, CG, r, res_CG, its_CG, res_CG] = pcg(A, b, tol, maxit);

fies(2)

semilog (res_J, hold on

semilog(res g); hold on

semilog(res cg) on

xlabel('iterations'); ylabel('Relative residuals');

legend('Jacobi', 'Gauss', 'Conjugate Gradient')

GRID on;

-for a n 10 20 for 20x solve the sys using CG (check how res; decreases by axiting flongs) (check also how the C ondsnun number c makes!!!GC1 dung using function eigs, and try to estimate it goth the Hkp of tol its

n = (10:10:100);

for r = 1:length(n),

A = gallery('poisson',ni ni);

b = ones(n(i), 1);

x0 = ones(n(i), 1); holp on;

Aiget) = eigs(A, b, tol, maxlit);

item = @(conv(res_JH)

end

this is finally equivalent to the existence of a vector v_i so that (N - I_d)v = 0.

this happens for y = 1/N; indeed, recalling that M = > * > , we have

d = 1 + 1 + . . . + 1 = > then > H - I_d(d)|1 = 0

1 - 1

Algoritmo PageRank

function [y] = pagerank (H, v, maxit)

N = size(H,1);

for k = 1:maxit

   y = H*y; y = y/||y||₁

end

This step (ct. Hoopes) y has already ||y||₁ = 1

Already normalized

dangling ones = d = H * ones (N,1)

d = d + N * dangling

H = y + x (x * dangling) * ones (1,N)

Algoritmo da usare in caso di esame PageRank

function [y] = pagerankbis (H, gamma, maxit, tol, v)

N = size(H,1);

d = sum(H,1);

dangling = (d == 0);

for i = 1:N

   H(i,:) = zeros(1,N)

   H(i,i) = ones(1,N);

end

dhat = d + N * dangling;

v = 1/N * ones (N,1);

y = rand(N,1);

r = y/sum(y);

erc = tol + 1;

while k <= maxit && erc > tol

   y_new= [ (gamma) * ((H*y)./dhat).+ + ((y./dhat) * diag(1) * ones(1,N)) * (1-gamma) ]

    * ones(N,1) * v^T ] + v ]

   erc = norm(y_new - y,1) / n];

   y = y_new;

   r = r+1;

end

end