Appunti Metodi numerici - Parte 2

ODE IV_02

function [t,t xm_f, Vm_f] = funzione

end

ODE IV_06

function [t, y] = ode (dt, t₀, y₀)t = 0:dt:t_f; t = t'; NT = length(t);y = zeros (NT,1);y(1) = y₀;for it = 2:NT    y(it) = y(it-1) + dt * ( ... )end

end

METODO DI HEUN

Metodo più accurato del Metodo di EuleroFino ad ora abbiamo studiato problemi in questa forma.^dy/_dt = e(y, t) es. * ^dy/_dt = -alpha-y

• dxm/dt = ...
• dvm/dt = -^k/_m...

Con il Metodo di Eulero sviluppiamo il seguente algoritmo:...

Questa algoritmo si porta a fare degli errori: O(1)

L'interpretazione geometrica del metodo di Eulero è:

L'errore è dato dal fatto che ci immaginiamo su una retta che ha pendenza che è quella iniziale della funzione.

Per ridurre l'errore devo ridurre il Δt.

Altri metodi sono quelli di Runge-Kutta e Adam Bashforth.

Uno dei più semplici della categoria dei metodi di Runge-Kutta è il METODO DI HEUN.

L'idea dei metodi di Runge-Kutta è quello di fare un passo in avanti usando un punto "migliorato" e non quello iniziale (come invece facevamo con il Metodo di Eulero).

METODO DI HEUN:

Uso, cioè, una pendenza media (che ai primi non sapevamo), quindi si usa prima la pendenza di partenza e calcolo un y_int che è in colore vi Tantòra y_int^E (questo è l'errore per calcolare).

yh(it) = yh(it-1) + dt * 1/2 * (RHS_iulz + RHS_fil)

tempo Heun = toc

% % Disegno dei vettori

figure

plot(y,e,'c'); ...

hold on

plot(, y, quadrico,'b')

plot( + yh,'k')

title('Prima eq. differenziale')

legend('...','...', 'soluzione con metodo di Heun')

t(ic) -> lo parte lineare

(tempo-toc) -> funzione lineare e on da il calcolo del tempo

L'accuratezza del metodo di Heun è più che il doppio rispetto al metodo di Eulero.

save "dove_stepHeun_terodix_03"

% Integrazione nel tempo con metodo di Heun

t(ic)

yh(1) = yj;

for it = 2: Nt

RHS_iul = -a * yh(it-1);

correzione

ypredittore = yh(t':+1) + dt * RHS_iulz; % se ora si chiama predittore

for i=1:2

RHS_fiu = -o * ypredittore;

ypredittore = yh(it-1) + dt * 1/2 (RHS_iuln + RHS_fin)

end

yh(it) = ypredittore;

end

Il sistema del 1^o ordine lo posso scrivere come:

M d²x/dt² + Ddx/dt + Kx = {ß}

Il sistema del 1^o ordine, invece, è:

m dy/dt + Ky = F

Quindi, possiamo dire che:

m = ^[2x3]_[4x4]^[2+3]_[4+3](^{I    O}_{O    M})    e    K = (^{O    I}_{K    D})

ESEMPIO:

ODE-IV.09 - Edificio multipiano, Metodo di Eulero

31/10/2022

sistema di 3 ODE del 1^o ordine lineare

d²x₁/dt² = -k₄/m₁x₁ + k₃/m₁(x₂ - x₁)

d²x₂/dt² = k₂/m₂(x₁ - x₂) + k₃/m₂(x₃ - x₂)

d²x₃/dt² = k₃/m₃(x₂ - x₃)

X = (x1x2x3)

Devo scrivere il sistema del 1^o ordine. Introduco il vettore delle incognite:

M d²x/dt² + Ddx/dt + Kx = {∫∫}

Metodo dell'attraversamento dello zero

ALGORITMO 1°

x(1) ? x(2)

x(2) ? x(3)

...

x(N-1) ? x(N)

• for ie = 1 : N - 1
• Confronto tra x(ie) e x(ie+1)
• se cambio segno → memorizzo ie

end

Tra x(ie) e x(ie+1) c'è un cambio segno?

1. if (x(ie) < 0 ∧ x(ie+1) >= 0) → cambio!
2. if (x(ie) >= 0 ∧ x(ie+1) < 0) → cambio!

Operatore relazionale

1. > 0 se x(ie) e x(ie+1) hanno stesso segno
2. <= 0 se x(ie) e x(ie+1) hanno segno diverso, ossia tra di loro c'è un punto e' zero

if x(ie) * x(ie+1) cambio!

se cambio segno memorizzo ie in icambio

ncambio = 0

for ie = 1 : N - 1

• if x(ie) * x(ie+1) <= 0
• ncambio = ncambio + 1 → aumenta ncambi di 1
• icambio(ncambio) = ie

end

function [G]= prepara_matrice(v)

N= length(v)-1;

G= zeros(N);

G(1,1)= v(1)+v(2);

G(1,2)= -v(2);

G(N,N)= v(N)+v(N+1);

G(N,N-1)= -v(N);

for ic=2:N-1

G(ic, ic-1)= -v(ic);

G(ic, ic)= v(ic)+v(ic+1);

G(ic, ic+1)= -v(ic+1);

end

end

%% Preparazione matrice MM e kk

MM= [eye(N), zeros(N), zeros(N), M];

kk= [zeros(N), -eye(N), k, D];

%% preparazione variabili

t=0:dt:Tf; t=t'; Nt=length(t);

Y= zeros(2*N,Nt); F= zeros(2*,Nt);

%% Integrazione nel Tempo (metodo di Eulero)

Y(1:N,1)= x0;

for it=1:Nt-1

Y(:,it+1)= Y(:,it)+dt * inv(MM) * (-kk*Y(:,it)+F(:,it));

end

disp(' Finito')

%% Disegno dei risultati

figure, plot(t,Y(1:N,:));

end

t(s)

• 0
• 5
• 9
• 11
• 15
• t_f = 30

Φ(t)

• 0
• 5
• 0
• 5
• 0
• 0

t_p = [0, 5, 9, 11, 15, t_f]

Φ_p = [0, 5, 0, 5, 0, 0]

Funzione di Matlab utile:

• Interp1(t_p, Φ_p, t)
• perché Φ dipende solo da t

Metodo di Eulero

^dx_dt = u_cosα → x^m+1 = x^m + Δt ⋅ (ν^m ⋅ cos(σ^m))

^dy_dt = u_senα → y^m+1 = y^m + Δt ⋅ (ν^m ⋅ sen(σ^m))

^dσ_dt = u_Tαρ → σ^m+1 = σ^m + Δt ⋅ (ν^m ⋅ tg(Φ^m) / L)

In Matlab:

• sim(x)
• cos(x) => x in rad
• sinα(x)
• cosd(x) => x in gradi

MATLAB:

• Save "Veicolo"
• Clear all;
• V = 10; L = 4;
• xo = 0; yo = 0; t_el0 = 0;
• Δt = 0.1; t_f = 30;