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:
<ul>
  <li>alpha: numero di Courant</li>
  <li>N+1: numero di nodi spaziali</li>
  <li>M+1: numero di nodi temporali</li>
  <li>Ui0(1:N+1,1): condizione iniziale</li>
  <li>U0j(1,1:M+1): condizione al bordo sinistro</li>
</ul>

Output:
<ul>
  <li>Uij(1:N+1,1:M+1): soluzione approssimata nei nodi</li>
</ul>

Inizializzazione del vettore delle approssimazioni
<pre><code>Uij = nan(N+1,M+1);</code></pre>

Condizione iniziale
<pre><code>Uij(:,1) = Ui0;</code></pre>

Condizione al bordo sinistro
<pre><code>Uij(1,:) = U0j;</code></pre>

Algoritmo del metodo upwind
<pre><code>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</code></pre>

<pre><code>%oppure
for j = 1:M
  Uij(2:N,j+1) = (1-alpha) * Uij(2:N,j) + alpha * Uij(1:N-1,j);
end
end</code></pre>

func upwind 2

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:
<ul>
  <li>xnodi(nnodi): vettore dei nodi</li>
  <li>fnodi(nnodi): funzione da approssimare calcolata nei nodi</li>
  <li>ndeg: grado del polinomio</li>
  <li>xeval: punto in cui valutare il polinomio</li>
</ul>

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 coefficientinnodi = length(xnodi);ak = nan(1,nord+1);bk = nan(1,nord);ak(1) = 2 / nnodi * sum( fnodi );for k=1:nordak(k+1) = 2 / nnodi * sum( fnodi .* cos(k*xnodi) );bk(k) = 2 / nnodi * sum( fnodi .* sin(k*xnodi) );end% Calcolo del polinomio trigonometricop_ndeg = ak(1) / 2;for k=1:nordp_ndeg = p_ndeg + ak(k+1) * cos(k*xeval) + bk(k) * sin(k*xeval);endendfunc approx trig 3METODO APPROSSIMAZIONEfunction [ 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

discretizzazione (nodi equispaziati)Ti = linspace(Tspan(1),Tspan(end),n_step+1);Ti = Ti';% Inizializzazione degli arrayYi = nan(n_step+1,1);% Condizione inizialeYi(1) = Y0;% Metodo di Heunfor i = 1:n_stepY_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;endendfunc heun1D 5FUNCTION HEUN Dfunction [ Ti, Yi ] = HeunD(ffun,I,Y0,n_step)% [ Ti, Yi ] = Heun(ffun,Tspan,Y0,n_step)% Soluzione di un sistema differenziale Y’=ffun(t,Y), Y(t0)=Y0% di D equazioni differenziali con il metodo di Heun%% INPUT% Tspan(1:2): estremi dell'intervallo di integrazione% n_step: numero di passi temporali% ffun(1,1:D): vettore contenente il termine noto del problema differenziale% Y0(1,1:D): condizioni iniziali.%% OUTPUT% Ti(1:n+1,1): nodi di discretizzazione% Yi(1:n+1,D): matrice contenente le approssimazioni% Passo di discretizzazioneh = (I(end)-I(1))/n_step;% Punti di discretizzazione (nodi equispaziati)Ti =

X = Thomas_diag(D,A,S,B)
-----------------------------------------
Input:
D: diagonale principale di A (vettore)
A: diagonale inferiore di A (vettore)
S: diagonale superiore di A (vettore)
B: termine noto

Output:
X: vettore soluzione

Verifica delle dimensioni dei vettori
in = length(D);
la = length(A);
ls = length(S);
lb = length(B);
if la~=ls || la~=n-1 || lb~=n
    error('Le dimensioni dei vettori non sono corrette')
end

Inizializzazione dei vettori Alfa, U, Y
Alfa = zeros(1,n-1);
U = zeros(1,n);
Y = zeros(1,n);

Fattorizzazione LU di A
U(1) = D(1);