Terzo esercizio sui metodi ONE-STEP e introduzione metodi alle differenze finite

Esercizio 3

Il delta9-tetraidrocannabinolo (THC), principale cannabinoide della cannabis, può essere somministrato a scopi terapeutici ma può causare dipendenza se assunto in dosi inappropriate. La farmacocinetica del THC medico e dei suoi metaboliti, dopo la somministrazione orale, può essere descritta da un modello compartimentale rappresentato in Figura 1.

Anche la farmacocinetica utilizza metodi numerici per l'approssimazione di soluzioni di grandi sistemi di equazioni differenziali. Quello che segue, da implementare su Matlab, ne è un esempio:

Il modello è descritto dal seguente sistema di equazioni differenziali:

dC_a(t)/dt = k_a·C_b(t) - (k_a + k_e)·C_a(t)
dC_b(t)/dt = -k_a·C_b(t)
dC₁(t)/dt = k_a·C_b(t) - (k₂₁ + k₂₃ + k₃₁ + k₁₇ + k₁₀ + k_e)·C₁(t) + k₁₂·C₂(t) + k₃₁·C₃(t) + k₅·C₅(t) + k₆₁·C₆(t)
dC₂(t)/dt = k₂₁·C₁(t) - (k₁₂ + k₂₃ + k₄₅ + k₃)·C₂(t) + k₃₂·C₃(t) + k₄₆·C₄(t)
dC₃(t)/dt = k₂₃·C₁(t) - (k₃₁ + k₃₂ + k₄₃)·C₃(t) + k₃₁·C₁(t) + k₅₃·C₅(t)
dC₄(t)/dt = k₃₁·C₃(t) - (k₄₃ + k₃₄ + k₃)·C₄(t) + k₃₄·C₆(t) + k₇₆·C₇(t)
dC₅(t)/dt = k₁₅·C₁(t) - (k₃ + k₂₅ + k₄₆ + k₅₃)·C₅(t) + k₃₆·C₃(t)
dC₆(t)/dt = k_weg>·(C_x(t) + k₃₀·C_x(t) + k_4c)·C₆(t) - C_c(t)·(k_5e + k₃₅·C_x(t))
dC₇(t)/dt = k_cdc·C_y(t) - k_rc·C_a(t)

Risoluzione numerica con Matlab

Per la risoluzione numerica del problema differenziale mediante i metodi One-Step e la function di Matlab Ode45, ho definito un intervallo di discretizzazione con l’ampiezza dell’intervallo di discretizzazione β pari a 720 minuti e numero di step n_step pari a 800 così da valutare le approssimazioni con un passo di discretizzazione h = 0.9. Come da richiesta, ho graficato le concentrazioni di THC (Delta9-tetraidrocannabinol) e dei metaboliti THC-OH (11-hydroxy-delta9-tetraidrocannabinol) e THC-COOH (11-nor-9-carboxy-delta9-tetraidrocannabinol) nel plasma in funzione del tempo:

Concentrazione di THC-COOH nel plasma

Anche dall'osservazione dei grafici precedenti si nota come il metodo di Eulero è quello che si discosta un po’ di più dalla soluzione più accurata (Ode45) anche se, a causa del piccolo passo di discretizzazione, non si riesce ad avere una netta distinzione delle curve. Per questo motivo, ho calcolato gli errori globali dei metodi rispetto alla soluzione ottenuta con la Ode45 nel caso del calcolo della concentrazione di THC nel plasma:

Il metodo di Eulero ha l’errore più alto in funzione del tempo rispetto al metodo di Heun e Runge-Kutta del 4° ordine, sottolineando che si tratta del metodo meno accurato.

Concentrazione di THC-OH

Concentrazione di THC-COOH

Esercizi sulle differenze finite

I metodi alle differenze finite consistono nell’approssimare ciascuna derivata che compare nelle equazioni differenziali con una opportuna formula alle differenze finite. Consiste nello sviluppare la funzione generica y(x) in serie di Taylor in un intorno di un i-esimo punto e combinare linearmente varie espressioni relative a diversi sviluppi ottenuti considerando vari intorni del punto. Questo processo consente di trasformare un problema continuo in un problema discreto per facilitare l’analisi.

Per applicare questi metodi, è necessario introdurre una discretizzazione dell’intervallo [a, b] suddividendolo in N+1 sottoinsiemi uguali, definendo così una griglia di nodi equispaziati:

x_i = a + ih
Con i = 0,1,...N+1 e h = (b-a)/(N+1) è il passo di discretizzazione.

Nel nostro caso abbiamo utilizzato il metodo delle differenze finite centrate per i nodi interni e delle differenze finite all’indietro e in avanti per le condizioni al bordo. Quindi, analizzando un problema descritto dalla seguente equazione differenziale:

y″(x) = p(x)y′(x) + q(x)y(x) − r(x)

verificato che:

• p(x), r(x) e q(x) ∈ C([a, b]);
• q(x) > 0;

possiamo affermare che la soluzione del problema lineare esiste ed è unica ed è quindi possibile approssimare le derivate nei nodi con uno sviluppo di serie di Taylor di ordine 2 per la derivata prima e di ordine 3 per la derivata seconda.

Nell’intorno successivo di Xi ho che:

y(x_i+1) = y(x_i + h) = y(x_i) + hy′(x_i) + (1/2!) h²y″(x_i) + (1/3!) h³y⁽³⁾(η)_i, η_i ∈ (x_i, x_i+1).

Mentre nell’intorno precedente ho che:

y(x_i−1) = y(x_i − h) = y(x_i) − hy′(x_i) + (1/2!) h²y″(x_i) − (1/3!) h³y⁽³⁾(η)_i, η_i ∈ (x_i−1, x_i).