Script di Analisi numerica

Approssimazione di zeri di funzioni reali

Bisezione

function [x,xK,it] = bise2 (fname,a,b,tol)

fa = fname(a)

fb = fname(b)

if sign(fa)*sign(fb) >= 0 % fa e fb concordi

error ('intervallo non corretto');

else

maxit = ceil(log((b-a)/tol)/log(2));

it = 0;

while it< maxit && abs(b-a) >= tol && eps * max((abs(a),abs(b)))

it = it + 1;

xK(it) = a + (b-a)*0.5;

fprintf ('it= %d x= %0.15f \n', it, xK(it));

fxK = fname(xK(it));

if fxK ==0

break

elseif sign(fxK)* sign(fa)> 0

a = xK(it);

fa = fxK ;

elseif sign(fxK)* sign(fb)> 0

b = xK(it);

fb = fxK ;

end

end

x = xK(it);

end

REGULA FALSI

(falsa posizione)

Function [x, xK, it] = regula_falsi (Fname, a, b, tol, nmax) fa = Fname(a); fb = Fname(b); if sign(fa)*sign(fb) >= 0 error ('intervallo non corretto'); else it = 0; err_rel = 1; fXK = Fname(a); while it < nmax && err_rel > tol && abs(fXK) >= tol it = it + 1; xK(it) = a - fa*(b-a)/(fb-fa); fprintf('it=%2d x=%8.15f \n', it, xK(it)); fXK = Fname(xK(it)); if fXK == 0 break elseif sign(fXK)*sign(fa) > 0 a = xK(it); fa = fXK; elseif sign(fXK)*sign(fb) > 0 b = xK(it); fb = fXK; end if it > 1 err_rel = abs(xK(it) - xK(it-1)/abs(xK(it))); else err_rel = 1; end end if it == nmax error ('Raggiunto numero massimo di iterazioni'); end end

break

end

end

if it == nmax

fprintf('raggiunto numero massimo di iterazioni \n');

end

NEWTON (modificato)

Function [x1, xK, it] = newtonm (Fname, Fpname, m, x0, tolx, tolf, nmax)

Fx0 = Fname (x0);

dFx0 = Fpname (x0);

if abs(dFx0) > eps

d = Fx0 / dFx0 ;

x1 = x0 - m*d ;

Fx1 = Fname (x1);

it = 1;

xK(it) = x1 ;

fprintf ('it: %d \t x: %9.8 15F \n', it, x1) ;

else

fprintf ('derivato nulla in x0 - EXIT \n')

return

end

while it < nmax && abs(Fx1) >= tolf && abs(d) >= tolx * abs(x1)

x0 = x1 ;

Fx0 = Fname (x0);

dFx0 = Fpname (x0) ;

if abs (dFx0) > eps

d = Fx0 / dFx0 ;

x1 = x0 - m*d ;

Fx1 = Fname (x1) ;

it = it + 1 ;

PIVOT TOTALE

Function [P, Q, L, U] = pivot_totale(A)

If det(A) == 0

warning('A è singolare')

elseif det(A) < 10^-6

warning('A è vicina a non essere singolare')

End

n = size(A, 2)

P = eye(n)

Q = P

For K = 1 : n - 1

[M, vec_idx_row] = max(abs(A(K:n, K:n)))

[~, idx_col] = max(M)

idx_row = vec_idx_row(idx_col) + K - 1

idx_col = idx_col + K - 1

A([K, idx_row], :) = A([idx_row, K], :) ;

P([K, idx_row], :) = P([idx_row, K], :) ;

A(:, [K, idx_col]) = A(:, [idx_col, K]);

Q(:, [K, idx_col]) = Q(:, [idx_col, K]);

A(K+1:n, K) = A(K+1:n, K) / A(K, K);

A(K+1:n, K+1:n) = A(K+1:n, K+1:n) - A(K+1:n, K) * A(K, K+1:n);

End

L = eye(n) + tril(A, -1);

U = triu(A);

End

HORNER NEWTON

function yy = HornerN(a, x, xx) N = length(a); yy = a(N) * ones(size(xx)); for K = N-1 : -1 : 1 yy = (xx - x(K)) .* yy + a(K); end end

INTERPOLANTE NEWTON

function yy = interp_N(x, y, xx) N = length(x); a = y; for K = 2 : N a(K:N) = (a(K:N) - a(K-1:N-1)) ./ (x(K:N) - x(1:N-K+1))i; end yy = HornerN(a, x, xx); end

PLOT

plot(xx, F(xx), 'r-', x, y, 'ro', xx, yy, 'b-') legend('vero', 'punti', 'polinomio interpolante')