Risolvere equazioni differenziali con MatLab
Sintassi generale per un file-M per la soluzione di un problema di cauchy:
function F=nomefunzione(t,y)
dove la funzione F, F=F(t,y) vale per il seguente problema di Cauchy $ => \begin{cases}
y'=F(t,y)\\
y(t_0)=y_0
\end{cases}$
se abbiamo il seguente problema di Cauchy:
$ \begin{cases}
y'=F(t,y)\\
y(t_0)=y_0
\end{cases}$
la risoluzione del problema è:
[t,y]=ode45('nomefunzione',tstan,y0)
dove t0 è l'estremo sinistro dell'intervallo.
Esempio 1:
Consideriamo il caso banale:
$ \begin{cases}
y'=t\\
y(0)=1
\end{cases}$
La soluzione cercata è $y= \frac {t^2}{2} +1$
Per risolvere il problema con matlab:
- Definisco una funzione di due variabili in un M-File come qui sotto:
- Definisco un intervallo di integrazione. In questo esempio cerchiamo una soluzione nell'intervallo [0,10]:
- L'istruzione usata , ode 23 è basata sul metodo di Runge-Kutta del 2 e 3 ordine:
- L’istruzione usata, ode45, è bastata sul metodo di Runge-Kutta del 4 e 5 ordine:
- Sono mostrati a schermo i valori di t e di y (soluzione numerica dell'equazione). Il programma non fornisce l'espressione della soluzione del problema di Cauchy ma un insieme di valori numerici da essa assunti.
- Plottiamo la soluzione trovata:
$ plot(t,y) $
- Poi troviamo il polinomio interpolante di secondo grado:
$ polyfit(t,y,2)
Esempio 2:
Studiamo il seguente problema di Cauchy:
$ \begin{cases} y'=y-t \\
y(0)=1 \end{cases0}$
Si tratta di un'equazione differenziale del primo ordine.
L'integrale generale è $ y=c e^x +x + 1 $
L'integrale particolare è y=x+1.
- Creo un file M che lo salvo con il nome: funzionesempiouno.m
- Ricerca della soluzione: (si preferisce l'ode45 perchè ha una maggiore precisione anche se fa un calcolo più complesso e lungo anche se con i processori moderno non si nota la differenza)
Esempio 3:
Questa volta prendiamo un'equazione differenziale del secondo ordine:
$ \begin{cases} y''+y=t^3 -2 \\
y(0)=0 \\
y'(0)=0
\end{cases}$
Per risolverla occorre trasformarla in un sistema di due equazioni del primo ordine ne seguente modo:
\begin{cases} y_1'=y_2 \\
y_2'=-y_1+t^3-2 \\
y_1(0)=0\\
y_2(0)=0
\end{cases}
Dove abbaimo posto $y_1=y e y_2=y'$
L'integrale generale è:
$ y(t)=A cos(t) + sin(t) + t_3 -6t-2$.
Risoluzione:
- e risolviamo come abbiamo fatto sopra...