Matlab - esame calcolo numerico con elementi di programmazione

Function Trapezi

function [ T_new, Nnodi, err ] = Trapezi(a, b, toll, f)

%%%% Approssimazione di un integrale definito con la formula dei trapezi

% e il criterio di Runge: [ T_N, Nnodi, err ] = Trapezi(a, b, rho, f_int)

Input:

• a, b: estremi dell'intervallo di integrazione
• rho: accuratezza
• f_int(x): funzione integranda

Output:

• T_N: approssimazione finale dell'integrale
• Nnodi: numero di nodi finale
• err: errore finale

%%% Approssimazione con la formula del trapezio

T_1 = (b-a) / 2 * ( f(a) + f(b) );

% Inizializzazioni

N = 1; % Numero di nodi iniziale: N+1

h = (b-a)/N; % Passo di integrazione

err = 10*toll;

T_old = T_1;

% Ciclo iterativo

while abs(err) > toll

N = N*2;

h = h/2;

xnodi = linspace(a,b,N+1);

T_new = h / 2 * ( f(xnodi(1)) + 2*sum(f(xnodi(2:N))) + f(xnodi(N+1)) );

err = (T_new - T_old) / 3;

T_old = T_new;

end

Nnodi = N+1;

end

Function UpWind

function [ Uij ] = UpWind ( alpha, N, M, Ui0, U0j )

%% Uij = UpWind(alpha, N, M, Ui0, U0j)

% Soluzione dell'equazione del trasporto con il metodo upwind

Input:

• alpha: numero di Courant
• N+1: numero di nodi spaziali
• M+1: numero di nodi temporali
• Ui0(1:N+1,1): condizione iniziale
• U0j(1,1:M+1): condizione al bordo sinistro

Output:

• Uij(1:N+1,1:M+1): soluzione approssimata nei nodi

%% Inizializzazione del vettore delle approssimazioni

Uij = nan(N+1,M+1);

% Condizione iniziale

Uij(:,1) = Ui0;

% Condizione al bordo sinistro

Uij(1,:) = U0j;

% Algoritmo del metodo upwind

for j = 1:M

for i = 2:N+1

Uij(i,j+1) = (1-alpha) * Uij(i,j) + alpha * Uij(i-1,j);

end

end

% oppure

% for j = 1:M

% Uij(2:N,j+1) = (1-alpha) * Uij(2:N,j) + alpha * Uij(1:N-1,j);

% end

end

Func Approx Trigonometrica

function [ p_ndeg, ak, bk ] = PolApproxTrig(xnodi, fnodi, nord, xeval)

% Costruisce il polinomio trigonometrico di grado ndeg per approssimare

% dati periodici: [p_ndeg,ak,bk] = PolApproxTrig(xnodi,fnodi,ndeg,xeval)

Input:

• xnodi(nnodi): vettore dei nodi
• fnodi(nnodi): funzione da approssimare calcolata nei nodi
• ndeg: grado del polinomio trigonometrico
• xeval(nx): vettore dei punti in cui calcolare il polinomio trigonometrico

Output:

• p_ndeg(nx): polinomio trigonometrico calcolato in xeval
• ak(ndeg+1),bk(ndeg): coefficienti del polinomio

%% Calcolo dei coefficienti

nnodi = length(xnodi);

ak = nan(1,nord+1);

bk = nan(1,nord);

ak(1) = 2 / nnodi * sum( fnodi );

for k=1:nord

ak(k+1) = 2 / nnodi * sum( fnodi .* cos(k*xnodi) );

bk(k) = 2 / nnodi * sum( fnodi .* sin(k*xnodi) );

end

% Calcolo del polinomio trigonometrico

p_ndeg = ak(1) / 2;

for k=1:nord

p_ndeg = p_ndeg + ak(k+1) * cos(k*xeval) + bk(k) * sin(k*xeval);

end

end

Metodo Approssimazione

function [ x_new, dist_k, k_iter ] = Metodo_approx( x_old, epsilon, N_max, phi )

% Metodo_approx: approssima il punto unito di una funzione con il metodo delle

% approssimazioni successive

%% [ x_fin, dist_fin, k_fin ] = Metodo_approx( x0, epsilon, N_max, phi )

Input:

• x0: approssimazione iniziale
• epsilon: accuratezza richiesta
• N_max: numero iterazioni
• phi: funzione di iterazione

Output:

• x_fin: approssimazione del punto unito
• dist_fin: accuratezza raggiunta
• k_fin: iterazioni esguite

% Inizializzazione del ciclo iterativo

dist_k = 10*epsilon;

k_iter = 1;

% Ciclo iterativo

while dist_k > epsilon && k_iter <= N_max

x_new = phi(x_old);

dist_k = abs(x_new-x_old);

k_iter = k_iter + 1;

x_old = x_new;

end

k_iter = k_iter - 1;

end

Function Heun 1D

function [ Ti, Yi ] = Heun1D(ffun,Tspan,Y0,n_step)

% [ Ti, Yi ] = Heun1D(ffun,Tspan,Y0,n_step)

% Soluzione di un'equazione differenziale del primo ordine

% ai valori iniziali con il metodo di Heun

Input:

• ffun(t,y): termine noto del problema differenziale (function)
• Tspan(1:2): estremi dell'intervallo di integrazione
• Y0(1): condizioni iniziali (vettore colonna)
• n_step: numero di passi temporali

Output:

• Ti(1:n_step+1): nodi di discretizzazione (equispaziati)
• Yi(1:n_step+1): matrice delle approssimazioni

% Passo di discretizzazione

h = (Tspan(end)-Tspan(1))/n_step;

% Punti di discretizzazione (nodi equispaziati)

Ti = linspace(Tspan(1),Tspan(end),n_step+1);

Ti = Ti';

% Inizializzazione degli array

Yi = nan(n_step+1,1);

% Condizione iniziale

Yi(1) = Y0;

% Metodo di Heun

for i = 1:n_step

Y_old = Yi(i);

K1 = ffun(Ti(i),Y_old);

K2 = ffun(Ti(i+1),Y_old+h*K1);

Yi(i+1) = Y_old + h * ( K1 + K2 ) / 2;

end

end