Risoluzione numerica di un'equazione differenziale lineare
Premetto che non ho un forte background in analisi numerica, per cui potrei commettere degli errori grossolani.
Motivato da alcuni problemi di elaborazione dei segnali, vorrei risolvere numericamente un ODE lineare di ordine $m$ del tipo:
\begin{equation*}
y(t)+ a_1(t) y'(t) + \dots + a_m(t) y^{(m)}(t) = b_0(t) x(t) + b_1(t) x'(t) + \dots + b_m(t) x^{(m)}(t),
\end{equation*}
dove $x$ è nota (input) e $y$ è incognita (output). Sia $f_s > 0$ la frequenza di campionamento e $y[n] = y(n/f_s)$ la soluzione approssimata da calcolare. Per brevità, possiamo porre:
\begin{equation*}
f(t) := b_0(t) x(t) + \dots + b_m(t) x^{(m)}(t).
\end{equation*}
Finora ho fatto così. Prima di tutto, ho ridotto l'equazione a un sistema differenziale lineare del primo ordine al solito modo:
\begin{equation*}
\begin{cases}
y_0(t) + a_1(t) y_1(t) + \dots + a_m(t) y_m(t) = f(t) \\
y'_0 (t) = y_1(t) \\
\dots \\
y'_{m-1}(t) = y_m(t). \\
\end{cases}
\end{equation*}
Poi ho usato il seguente metodo lineare 1-step:
\begin{equation} \label{eq:method}
y'[n] \simeq (\alpha+1) f_s \ y[n] - (\alpha+1) f_s \ y[n-1] - \alpha \ y'[n-1],
\end{equation}
dove $\alpha \in [0,1]$. Si ottiene così il sistema lineare
\begin{cases}
y_0[n] + a_1[n] y_1[n] + \dots + a_m[n] y_m[n] = f[n] \\
y_{i+1}[n] = (\alpha+1) f_s \ y_i[n] - (\alpha+1) f_s \ y_i[n-1] - \alpha \ y_{i+1}[n-1], \hspace{1 cm} i = 0, \dots, m-1 \\
\end{cases}
e la soluzione è data ricorsivamente da:
\begin{cases}
y_0[n] = \dfrac{f[n] + \sum_{i=1}^m \{ a_i[n] \sum_{j=1}^i (\alpha+1)^{i-j} f_s^{i-j}((\alpha+1) f_s \ y_i[n-1] + \alpha \ y_{i+1}[n-1])\}}{1 + \sum_{i=1}^m a_i[n] (\alpha+1)^i f_s^i} \\
y_{i+1}[n] = (\alpha+1) f_s \ y_i[n] - (\alpha+1) f_s \ y_i[n-1] - \alpha \ y_{i+1}[n-1], \hspace{1 cm} i = 0, \dots, m-1 \end{cases}
Ora vi espongo i miei dubbi.
1) Il procedimento è corretto?
2) La riduzione a un sistema di equazioni del primo ordine può inficiare la qualità della soluzione? Quali sono le alternative? Ad esempio, come si comportano i metodi multistep?
3) Il metodo lineare 1-step più generale è:
\begin{equation*}
y'[n] \simeq k_1 y[n] + k_2 y[n-1] + k_3 y'[n-1],
\end{equation*}
dove $k_1, k_2, k_3 \in \mathbb R$. L'espressione (\ref{eq:method}) si ottiene imponendo che il metodo sia almeno del primo ordine, mentre la condizione $\alpha \in [0,1]$ serve ad evitare problemi di stabilità. Per $\alpha = 0$ si ottiene Eulero all'indietro, per $\alpha = 1$ il metodo del trapezio (per $y'$) e per $\alpha \to \infty$ Eulero in avanti. Questo metodo è citato in alcuni articoli del campo applicativo su cui sto lavorando, ma non ne ho trovato traccia nella letteratura matematica. Conoscete qualche referenza? Come si comporta, a livello di efficienza e accuratezza, rispetto a metodi più usuali?
4) In generale, quali referenze mi consigliate per metodi numerici per ODE, in particolare lineari?
Grazie in anticipo.
Motivato da alcuni problemi di elaborazione dei segnali, vorrei risolvere numericamente un ODE lineare di ordine $m$ del tipo:
\begin{equation*}
y(t)+ a_1(t) y'(t) + \dots + a_m(t) y^{(m)}(t) = b_0(t) x(t) + b_1(t) x'(t) + \dots + b_m(t) x^{(m)}(t),
\end{equation*}
dove $x$ è nota (input) e $y$ è incognita (output). Sia $f_s > 0$ la frequenza di campionamento e $y[n] = y(n/f_s)$ la soluzione approssimata da calcolare. Per brevità, possiamo porre:
\begin{equation*}
f(t) := b_0(t) x(t) + \dots + b_m(t) x^{(m)}(t).
\end{equation*}
Finora ho fatto così. Prima di tutto, ho ridotto l'equazione a un sistema differenziale lineare del primo ordine al solito modo:
\begin{equation*}
\begin{cases}
y_0(t) + a_1(t) y_1(t) + \dots + a_m(t) y_m(t) = f(t) \\
y'_0 (t) = y_1(t) \\
\dots \\
y'_{m-1}(t) = y_m(t). \\
\end{cases}
\end{equation*}
Poi ho usato il seguente metodo lineare 1-step:
\begin{equation} \label{eq:method}
y'[n] \simeq (\alpha+1) f_s \ y[n] - (\alpha+1) f_s \ y[n-1] - \alpha \ y'[n-1],
\end{equation}
dove $\alpha \in [0,1]$. Si ottiene così il sistema lineare
\begin{cases}
y_0[n] + a_1[n] y_1[n] + \dots + a_m[n] y_m[n] = f[n] \\
y_{i+1}[n] = (\alpha+1) f_s \ y_i[n] - (\alpha+1) f_s \ y_i[n-1] - \alpha \ y_{i+1}[n-1], \hspace{1 cm} i = 0, \dots, m-1 \\
\end{cases}
e la soluzione è data ricorsivamente da:
\begin{cases}
y_0[n] = \dfrac{f[n] + \sum_{i=1}^m \{ a_i[n] \sum_{j=1}^i (\alpha+1)^{i-j} f_s^{i-j}((\alpha+1) f_s \ y_i[n-1] + \alpha \ y_{i+1}[n-1])\}}{1 + \sum_{i=1}^m a_i[n] (\alpha+1)^i f_s^i} \\
y_{i+1}[n] = (\alpha+1) f_s \ y_i[n] - (\alpha+1) f_s \ y_i[n-1] - \alpha \ y_{i+1}[n-1], \hspace{1 cm} i = 0, \dots, m-1 \end{cases}
Ora vi espongo i miei dubbi.
1) Il procedimento è corretto?
2) La riduzione a un sistema di equazioni del primo ordine può inficiare la qualità della soluzione? Quali sono le alternative? Ad esempio, come si comportano i metodi multistep?
3) Il metodo lineare 1-step più generale è:
\begin{equation*}
y'[n] \simeq k_1 y[n] + k_2 y[n-1] + k_3 y'[n-1],
\end{equation*}
dove $k_1, k_2, k_3 \in \mathbb R$. L'espressione (\ref{eq:method}) si ottiene imponendo che il metodo sia almeno del primo ordine, mentre la condizione $\alpha \in [0,1]$ serve ad evitare problemi di stabilità. Per $\alpha = 0$ si ottiene Eulero all'indietro, per $\alpha = 1$ il metodo del trapezio (per $y'$) e per $\alpha \to \infty$ Eulero in avanti. Questo metodo è citato in alcuni articoli del campo applicativo su cui sto lavorando, ma non ne ho trovato traccia nella letteratura matematica. Conoscete qualche referenza? Come si comporta, a livello di efficienza e accuratezza, rispetto a metodi più usuali?
4) In generale, quali referenze mi consigliate per metodi numerici per ODE, in particolare lineari?
Grazie in anticipo.
Risposte
2): Per i metodi principali per ODE è necessario avere un sistema del prim'ordine, e questo non inficia in alcun modo sull'accuratezza del risultato.
Inoltre, quando scrivi $y'[n]$, qui stai dicendo che la tua nuova incognita è la **funzione** $y[n]$, giusto? Intendo dire: non hai ancora discretizzato nulla? Perché se non è così allora quella derivata non ha senso.
Quando poi scrivi
onestamente non ho capito quale sia l'obiettivo. Dovresti essere più specifico e scrivere esattamente il sistema differenziale da risovlere.
Inoltre, quando scrivi $y'[n]$, qui stai dicendo che la tua nuova incognita è la **funzione** $y[n]$, giusto? Intendo dire: non hai ancora discretizzato nulla? Perché se non è così allora quella derivata non ha senso.
Quando poi scrivi
"pier.paolo":
Poi ho usato il seguente metodo lineare 1-step:
\begin{equation} \label{eq:method} y'[n] \simeq (\alpha+1) f_s \ y[n] - (\alpha+1) f_s \ y[n-1] - \alpha \ y'[n-1], \end{equation}
onestamente non ho capito quale sia l'obiettivo. Dovresti essere più specifico e scrivere esattamente il sistema differenziale da risovlere.
"feddy":
Inoltre, quando scrivi $ y'[n] $, qui stai dicendo che la tua nuova incognita è la **funzione** $ y[n] $, giusto? Intendo dire: non hai ancora discretizzato nulla? Perché se non è così allora quella derivata non ha senso.
Sono partito da un problema continuo e ora voglio discretizzare, per cui ho bisogno di un modo per approssimare le derivate. E' quello che si fa di solito: ad esempio, il metodo di Eulero in avanti consiste nell'approssimare
\[
y'[n] \simeq f_s (y[n+1]-y[n]).
\]
"feddy":
Quando poi scrivi [quote="pier.paolo"]
Poi ho usato il seguente metodo lineare 1-step:
\[ \begin{equation} \label{eq:method} y'[n] \simeq (\alpha+1) f_s \ y[n] - (\alpha+1) f_s \ y[n-1] - \alpha \ y'[n-1], \end{equation} \]
onestamente non ho capito quale sia l'obiettivo. Dovresti essere più specifico e scrivere esattamente il sistema differenziale da risovlere.[/quote]
Per ora non ho in mente un problema specifico da risolvere, anche se posso pensare a degli esempi e nel caso postarli. Intendevo che, utilizzando quell'espressione per la derivata e sostituendola a tutte le derivate del sistema differenziale, si ottiene un sistema lineare ad ogni passo. I valori trovati per $y[n]$ forniscono poi la soluzione approssimata.
Sì certo, hai approssimato le derivate con differenze finite, ma comunque ci sono ancora cose che non tornano. Ad esempio
questo non è un sistema differenziale del prim'ordine, a sinistra di ogni uguale devi avere una derivata prima.
Non ho letto nel dettaglio il tuo procedimento, ma una volta che uno ha un sistema di n ODE del prim'ordine $Y' = F(t,Y)$, dove $Y = [y_1(t),\ldots, y_n(t)]$, poi ogni metodo può (in teoria ) essere applicato. Ad esempio un Eulero esplicito sarebbe (partenedo da un dato iniziale $Y_0= [y_1(t_0), \ldots, y_n(t_0)]$)
$Y_{n+1} = Y_n + \Delta t F(t_n,Y_n)$, per $n \geq0$, dove $\Delta t$ credo che nel tuo caso sia $f_s$ e $F$ è la funzione vettoriale che descrive il membro di destra che conosci in vari istanti (visto che a quanto pare è qualcosa di campionato).
Ovviamente Eulero esplicito è da evitare per problemi di stabilità, mentre Eulero implicito questi problemi non ne ha e visto che l'equazione è lineare ad ogni passo dovrai risolvere un sistemi lineare, che però non è quello che fai te nel tuo schema
\[ \begin{equation*} \begin{cases} y_0(t) + a_1(t) y_1(t) + \dots + a_m(t) y_m(t) = f(t) \\ y'_0 (t) = y_1(t) \\ \dots \\ y'_{m-1}(t) = y_m(t). \\ \end{cases} \end{equation*} \]
questo non è un sistema differenziale del prim'ordine, a sinistra di ogni uguale devi avere una derivata prima.
Non ho letto nel dettaglio il tuo procedimento, ma una volta che uno ha un sistema di n ODE del prim'ordine $Y' = F(t,Y)$, dove $Y = [y_1(t),\ldots, y_n(t)]$, poi ogni metodo può (in teoria ) essere applicato. Ad esempio un Eulero esplicito sarebbe (partenedo da un dato iniziale $Y_0= [y_1(t_0), \ldots, y_n(t_0)]$)
$Y_{n+1} = Y_n + \Delta t F(t_n,Y_n)$, per $n \geq0$, dove $\Delta t$ credo che nel tuo caso sia $f_s$ e $F$ è la funzione vettoriale che descrive il membro di destra che conosci in vari istanti (visto che a quanto pare è qualcosa di campionato).
Ovviamente Eulero esplicito è da evitare per problemi di stabilità, mentre Eulero implicito questi problemi non ne ha e visto che l'equazione è lineare ad ogni passo dovrai risolvere un sistemi lineare, che però non è quello che fai te nel tuo schema
Ok, sono andato un po' veloce in effetti. Nel sistema
\begin{equation*}
\begin{cases} y_0(t) + a_1(t) y_1(t) + \dots + a_m(t) y_m(t) = f(t) \\
y'_0 (t) = y_1(t) \\
\dots \\
y'_{m-1}(t) = y_m(t). \\
\end{cases}
\end{equation*}
ho aggiunto l'ulteriore "incognita" $y_m$ giusto perché non mi piaceva lasciare $y'_{m-1}$ nella prima equazione. La forma normale si ottiene cancellando l'ultima equazione e ricavando $y_m = y'_{m-1}$ dalla prima. Bisogna supporre anche che $a_m(t)$ sia sempre non nulla affinché tutto torni.
\begin{equation*}
\begin{cases} y_0(t) + a_1(t) y_1(t) + \dots + a_m(t) y_m(t) = f(t) \\
y'_0 (t) = y_1(t) \\
\dots \\
y'_{m-1}(t) = y_m(t). \\
\end{cases}
\end{equation*}
ho aggiunto l'ulteriore "incognita" $y_m$ giusto perché non mi piaceva lasciare $y'_{m-1}$ nella prima equazione. La forma normale si ottiene cancellando l'ultima equazione e ricavando $y_m = y'_{m-1}$ dalla prima. Bisogna supporre anche che $a_m(t)$ sia sempre non nulla affinché tutto torni.
In ogni caso devi scriverlo prima nella forma $Y' = F(t,Y)$ se vuoi usare qualsiasi metodo classico per ODE. Se hai qualche referenza su quello che stai facendo magari riesco ad esserti più d'aiuto
Oltretutto, essendo $x$ campionato pure quello, non è per niente ovvio cosa sia $x^{(m)}$. Dovresti usare un'interpolante che ti garantisca un'ordine di derivazione almeno $m$, e perciò piuttosto elevato, per dare senso ad un espressione del genere.
"feddy":
Oltretutto, essendo $x$ campionato pure quello, non è per niente ovvio cosa sia $x^{(m)}$. Dovresti usare un'interpolante che ti garantisca un'ordine di derivazione almeno $m$, e perciò piuttosto elevato, per dare senso ad un espressione del genere.
Se sto discretizzando un sistema continuo, $x$ è una funzione nota e ne posso calcolare le derivate, sotto opportune ipotesi di regolarità. Se invece il mio sistema nasce discreto, per $x$ uso la stessa approssimazione che ho scritto per $y$, e per le derivate superiori itero semplicemente il procedimento. In effetti avevo anche il dubbio che questo potesse condurre a errori troppo grandi.
Beh ma se $x$ è discreto, doverne fare la derivata (ad esempio) quarta con differenze finite richiede che la funzione sia abbastanza regolare e soprattutto rende molto più difficile il tutto. E questa credo sia la situazione in cui ti trovi da quello che ho capito.
Secondo me, prima di qualsiasi cosa, è bene che tu abbia chiaro in mente come risulta il sistema differenziale discretizzato.
Poi, c'è una famiglia di metodi, integratori di tipo Magnus, che potrebbero fare al caso tuo. Questi usano l'esponenziale di matrice, ma per il momento è meglio prima sistemare quanto hai scritto sopra
Secondo me, prima di qualsiasi cosa, è bene che tu abbia chiaro in mente come risulta il sistema differenziale discretizzato.
Poi, c'è una famiglia di metodi, integratori di tipo Magnus, che potrebbero fare al caso tuo. Questi usano l'esponenziale di matrice, ma per il momento è meglio prima sistemare quanto hai scritto sopra
Se $x$ è campionato, puoi assumere tutta la regolarità che vuoi. Avevo già letto degli integratori di tipo Magnus, di che si tratta?
"pier.paolo":
Se x è campionato, puoi assumere tutta la regolarità che vuoi.
Onestamente questo non credo sia sempre vero, dipende da dove viene. Se lo assumi di sicuro hai una motivazione valida, però in generale non credo sia vero.
Sono integratori esponenziali per problemil lineari. Per usarli devi avere il problema scritto nella forma $Y'(t) = A(t) Y(t)$, dove $A(t)$ è la matrice che descrive il tuo sistema differenziale e può dipendere dal tempo. Ad ogni modo, per essere implementati in modo efficente richiedono delle accortezze particolari sulla gestione dell'esponenziale di matrice.
Al di là del metodo, però, devi scrivere il sistema come ho detto nei post precedenti, altrimenti non se ne fa niente
Ok, ho letto anche qui: https://en.wikipedia.org/wiki/Magnus_expansion. Tuttavia il mio problema non è omogeneo, c'è modo di adattare il metodo?
Sì, ma per iniziare devi scrivere per bene il sistema. Provare a usare un metodo semplice, per iniziare ed accertarsi di aver scritto le cose per bene. Poi magari puoi usare metodi più avanzati.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo