Risoluzione numerica di un'equazione differenziale lineare

[pier.paolo]

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.

[feddy]

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

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.

[pier.paolo]

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]).
\]

Quando poi scrivi

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.

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.

[feddy]

Sì certo, hai approssimato le derivate con differenze finite, ma comunque ci sono ancora cose che non tornano. Ad esempio

\[ \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

[pier.paolo]

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.

[feddy]

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

[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.

[pier.paolo]

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.

[feddy]

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

[pier.paolo]

Se $x$ è campionato, puoi assumere tutta la regolarità che vuoi. Avevo già letto degli integratori di tipo Magnus, di che si tratta?