Appunti di Metodi numerici - Parte 1

PROF. GIORGIO BELLOTTI 19/09/2022

2 ESONERI: NOVEMBRE E DICEMBRE

Un'equazione differenziale semplice (soluzione analitica)

^dy/_dt = -a.y

y(t) e' reale

t e' reale

a e' reale positivo

la grandezza y dipende dal tempo

variabile dipendente

variabile indipendente

parametro del problema

Quando y→0 la variazione di y nel tempo tende ad essere molto piccola.

Per risolvere l'eq. differenziale, e quindi calcolare la variazione dy, dobbiamo specificare una condizione iniziale:

t=0 y(0)=y₀ → condizione iniziale IC

Per questa eq. c'è una soluzione analitica perchè è un'eq. differenziale ordinaria lineare (ODE) ed è:

y_a(t) = y_f e^-a.t → soluzione analitica

Non sempre abbiamo una soluzione analitica.

Tecniche di soluzione numerica

• La soluzione è approssimata (produce un certo errore)
• Le soluzioni numeriche sono funzioni discrete mentre la soluzione analitica è una funzione continua nel tempo.

Cerca la successione di valori: y₁, y₂, y₃, y₄, y₅ e non cerca per la funzione y(t)

Mi devo riscrivere l'equazione differenziale in funzione dei valori y₄, y₃, y₂, y₁, y₀.

RHS scritto al livello temporale 3:

RHS scritto al livello temporale 3:2y → -c. y₃

LHS scritto al livello temporale 3:

_dy = lim (y(t₃+Δt)-y(t₃)) _dt Δt→∞ Δt

= lim (y(t₄)-y(t₃))Δt→∞ Δt

= lim (y₄-y₃ ≃ y₁-y₃)Δt→∞ Δt Δt

Allora _dy = -c. y₃ si può scrivere in maniera approssimata come: _dt

y₄-y₃ = -c. y₃ → abbiamo approssimato e discretizzato Δt

In generale: y_i+1-y_i = -c. y_i → FORMULA ALGEBRICA Δt

La soluzione è: y_i+1=y_i+Δt (-c. y_i) → derivativo di avanzato del tempo→ METODO DI EULERO

Sottraggo i membri di sinistra con quelli di destra tra le 2 equazioni (1 e 2) e ottengo:

Esplicito _dy/_dx

_dy/_dx|xi = (_g(xi+Δx) – _g(xi–Δx)) / 2Δx

– (2/3)_Δx²/2 x³_i Δx³ / 3!

+ R_xi

Formula differenza approssimata per il calcolo delle derivate derivate

Per calcolare la derivata in x_i: considero un punto davanti x_i+Δx e un punto dietro x_i–Δx. → SCHEMA CENTRATO (DEL 2° ORDINE)

L’errore in questo schema di calcolo approssimato della derivata ha un ordine di Δx² quindi in questo caso l’errore è più piccolo rispetto al caso precedente (l’errore si annunisce in maniera quadratica).

Per calcolare le derivate seconde devo sommare le 2 equazioni, ma ci ritorneremo più avanti! E.g. Per calcolare salvo le 2 equazioni ottengo:

Esplicito _d²y/_dx²

_d²y/_dx²|xi:

(y(x_i+Δx)+y(x_i–Δx)–2g(x_i)) + (x⁴) / 2!

Δx⁴/x_i⁴ + R_xi

Non ci sono grandi differenze nella struttura dei 2 problemi, quindi possiamo applicare il Metodo di Eulero anche per il sistema di eq. differenziali che dobbiamo appena risolvere; l'unico problema è che quando avviamo altri nel tempo dobbiamo procedere in avanti nel tempo risolvendo entrambe le equazioni mentre per risolvere le 2^a eq. ho bisogno della 1^a eq.

Scrivo in forma discreta il sistema di eq. differenziali sfruttando quello che abbiamo studiato precedentemente il calcolo delle derivate in forma approssimata, quindi non cerchiamo la soluzione del problema nella sua dipendenza continua dal tempo ma cerchiamo piuttosto il valore che le funzioni assumono a specifici istanti di calcolo (cioè discretizzo il dominio del tempo in t₁, t₂, t₃, ...).

• X_m(t) → X_m(1), X_m(2), ...
• U_m(t) → U_m(1), U_m(2), ...

Scrivo le 2 eq. in forma discreta sfruttando il calcolo approssimato delle derivate:

^dmX_m^|_dt^|_ti = (X_m(i+1) - X_m(i)) / Δt + ERRORE o(Δt)

(U_m(i+1) - U_m(i)) / Δt = (-k/m X_m(i) - μ/m U_m(i) + S(l)/m) → anula e risolvi algebraica (numerica)

Esplicito U_m(i+1).

U_m(i+1) = U_m(i) + Δt · (-k/m X_m(i) - μ/m U_m(i) + S(l)/m)

Faccio la stessa cosa per la 1^a eq:

dX_m / dt = U_m → (X_m(i+1) - X_m(i)) / Δt = U_m(i)

MATLAB

script = file di testo con estensione .m

Poco scrivere uno script da 0: New script -> save -> "Primo script"

Cose prime cosa pulisco la memoria: clear -> clc -> close

dt = 0.1;

a = 2.1;

tf = 4.0;

Le dico Run mando in esecuzione lo script oppure posso scrivere sul Command Window il nome della script "Primo script"

Costruisco un vettore t che va da 0, ogni dt, fino a 10.

tf = 10;

t = 0 : dt : tf; -> vettore riga

Preferisco lavorare con un vettore colonna:

t = t' -> vettore colonna

CICLO FOR

for ii = [3,5]

disp (ii)

end

ii = variabile ausiliaria contatore per il ciclo che assume i valori 3 e 5

Quindi ciclo sarà ripeto 2 volte e c'è l'ombre si eseguono i valori 3 e 5

Ad ogni iterazione del ciclo sarà eseguita istruzione che prevede di far vedere e disploy (nella Command window) il valore di ii

Mando in esecuzione e beppo 3 e 5 sulla command window

for ii = t

disp(ii)

end .. -> ottengo i valori del vettore t

Se metto uno smorzamento piccolo l’ampiezza delle oscillazioni rimane quasi costante e crea un cerchio nel piano delle forze dove si conserva l’energia.

Se metto μ = 0 ad esempio che il sistema conservi energia e che nel piano delle forze ho una circonferenza, invece essendo che lo stato del sistema si adatta, ho una spirale che diverge all’infinito e le oscillazioni aumentano nel tempo.

Che metto μ = 0.1 e dt = 0.05 e inserisco una forza esterna che varia nel tempo, ad esempio una forzante armonica.

f(t) = Δ・sen(ωt), dove ω = 2π/T

ω ≈ frequenza o frequenza angolare

2π/T ≈ periodo delle oscillazioni

Δ=1;

f=Δ・sin(ω_exc・t);

Y = zeros(2, Nt)

Y(:, 1) → caricare le coord. iniziali

_{estraggo la -esima colonna della matrice Y}

Y₂ → Y(:, 2)

Y_i → Y(:, i_t)

X_t+1 → Y(:, i_t+1)

m → H = [1 0, 0, m]

oppure: comando "diag"

k → k = [0, -1, k, μ]

F → F

_{matrice con 2 righe e Nt colonne}

f → F = zeros(2, Nt)

_trasposta

F(1, iⱼ) = f^j

_inizializzo

Ciclo di avanzamento nel tempo

Y(1, 1) = X_mi; Y(2, 1) = V_mi

for i