Sistema equazioni differenziali - Integrazione numerica

Messaggioda Studente35 » 04/01/2019, 20:42

Buonasera a tutti, mi trovo ad affrontare un problema del quale non riesco a comprenderne pienamente la risoluzione.
Illustrerò il problema ed esporrò quanto ho svolto io.
Vorrei prima delinearlo matematicamente "a penna", per poi implementarlo in matlab.
Devo risolvere un sistema con il metodo di Heun, metodo che conosco abbastanza bene e che hoapplicato ad equazioni differenziali del primo ordine. La difficoltà, per me, sta nel modo in cui considerare il vettore m, costituito dai due vettori m1 ed m2, aventi 3 componenti ciascuno. Tali vettori dovrebbero rappresentare le velocità, nello spazio (ergo con 3 componenti) di due campi magnetici. Il valore di tutti i termini che compaiono sono riportati nelle immagini prese da matlab alla fine del post.
Il sistema è il seguente

Immagine

di seguito riporto alcune procedure svolte:

sviluppo di alcuni prodotti vettoriali
Immagine

che sostituiti:

Immagine

portando i termini con le derivate al primo membro:


Immagine

considerando solo il primo membro:


Immagine

che diventerà:


Immagine

la matrice che compare è quella che in matlab, ho chiamato "A"

dopo questo, ho svolto in matlab il secondo membro e ho portato la matrice A al secondo membro.

Ho iniziato ad applicare Heun ma mi sono bloccato al primo step (EULERO).

Qualunque dritta possiate darmi, è veramente ben accetta.
riporto di seguito la sintassi utilizzata in matlab.
le condizioni iniziali dovrebbero essere i vettori m1 ed m2 che riporto in matlab (righe 6 e 7).


Immagine


Immagine



Immagine

Vi ringrazio anticipatamente.
Studente35
Starting Member
Starting Member
 
Messaggio: 2 di 12
Iscritto il: 11/12/2018, 10:09

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda feddy » 05/01/2019, 15:21

Ciao,

premetto che non ho il tempo per mettermi a controllare tutti i dettagli, ma ho dato un'occhiata.
Mi pare di vedere che tu sia interessato alla risoluzione di un problema del tipo $A \frac{dm}{dt}= F(\vec{m})$, dove $F$ è ovviamente una funzione a valori vettoriali. $A$ viene chiamata solitamente "matrice di massa", e ci sono solver dedicati a questo tipo di problemi in MatLab/Octave (tipo ode23t).
In particolare, l'inversione della matrice e il passaggio da $A(\vec{m}) \frac{dm}{dt}=F(m)$ a $\frac{dm}{dt}=A(m)^{-1} F(m)$, tramite il comando backslash, non è in genere una buona scelta. In questo caso inoltre la tua matrice dipende esplicitamente da $m=(m_1, \ldots,m_6)$.

Se non hai molto tempo, una strategia potrebbe essere quella di leggere la documentazione di un solver per problemi con matrici di massa, e passare la funzione e la matrice del problema in esame, che è una scemata ovviamente perché fai fare tutto a MatLab, ma avresti la sicurezza dell'implementazione.
Se vuoi farlo da solo e bene, allora la faccenda si complica perché un'implementazione naive di quel tipo non è detto sia corretta

Per la documentazione, leggi qui, e questa utile pagina "choose an ODE solver": click

Bada che il metodo di Heun che tu proponi, è comunque un metodo esplicito, e quindi di facile implementazione1. A tal proposito, mi sembra che il tuo codice possa essere migliorato. In particolare, io fossi in te, definirei la funzione $F$ una volta per tutti. Ossia
Codice:
F=@(m) [m(1) - m(2);m(3)-m(5);m(6)-m(2);m(3);m(1)-m(4);]
dove ovviamente ho messo valori a caso in ogni entrata. La parte noiosa è quella di inserire in ognuna della sei componenti i vari prodotti vettoriali che dipenderanno da $m$, ma con comandi analoghi a quelli fatti da te. Insomma, alla fine ti ritroverai con una funzione vettoriale $F(\vec{m})$ ben precisa, un termine noto $m_0$ che definirai come un vettore a sei componenti che è noto e ti ritroverai a risolvere il problema

$A(\vec{m}) \frac{dm}{dt}=F(\vec{m})$


con dato iniziale $\vec{m_0}=[m_1,\ldots,m_6]$.

Bada bene che $A$ è una matrice "state dependent", quindi per l'eventuale inversione (a.k.a risoluzione di un sistema lineare ad ogni time-step) dovresti sapere a priori che questa è non singolare ad ogni time step, e questo può essere un problema. Ad ogni modo, seguendo il tuo ragionamento, ti riconduci al problema $\vec{m}'=A^{-1}(\vec{m}) F(\vec m)= \tilde{F}(\vec{m})$ e applichi Heun al problema $\vec{m}'=\tilde{F}(\vec{m})$, cioè

$y_{n+1}=y_n + \frac{h}{2}(\tilde{F}(y_n) + \tilde{F}(y_n+h \tilde{F}(y_n)))$


che con un banalissimo ciclo for si risolve.

Tirando lo somme, io farei così:

1. Prima, se hai effettivamente bisogna dei risultati, usa un solutore di MatLab. Sicuramente la tua soluzione sarà meno accurata, senza contare che MatLab utilizza un passo variabile,e mille altre cose che rendono il codice realemente efficiente.

2. Prova la risoluzione che ti ho proposto, che altro non è che la tua idea scritta in modo più schematizzato. A mio avviso, il tuo codice è un po' difficile da leggere, nel senso che quel ciclo per il metodo di Heun, facendo come l'ho proposto io, risulta effettivamente chiaro, e tutta la complessità della scrittura del codice viene mappata nello scrivere la funzione $F(\vec{m})$ una volta per tutte, ma non è la fine del mondo.

In bocca al lupo

Note

  1. fosse stato implicito, la faccenda sarebbe molto più difficile a causa della presenza della matrice di massa
Ultima modifica di feddy il 12/01/2019, 01:13, modificato 1 volta in totale.
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2350 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda Studente35 » 05/01/2019, 16:06

Grazie mille per la risposta e i consigli che mi hai dato. L'obiettivo era quello di non usare una funzione predisposta in matlab effettivamente. Quindi vedrò di seguire i tuoi consigli. In effetti era proprio questo che mi bloccava un poco, cioè il fatto che anche A sia dipendente da m e non mi convinceva la definizione di F, la vedevo caotica. Grazier ancora, ti aggionerò sul lavoro.
Studente35
Starting Member
Starting Member
 
Messaggio: 3 di 12
Iscritto il: 11/12/2018, 10:09

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda feddy » 07/01/2019, 10:11

Bene, buon lavoro! ;)
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2360 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda Studente35 » 11/01/2019, 21:30

feddy ha scritto:Bene, buon lavoro! ;)


Ciao, ho provato a lavorarci su. Analiticamente, per sommi capi, ho messo su carta i procedimenti da fare in maniera generica, sono corretti?

Immagine
ovviamente, il procedimento che ho fatto per t=0 e t=1 deve essere ripetuto per tutti i time-step attraverso il ciclo for.

Ho provato a implementare il tutto in Matlab, e ho ottenuto 6 curve, che mi aspettavo di ottenere, ma che sono costanti dopo t=1. Credo ci sia un problema nell'indicizzazione, utilizzo, per esempio, (:,ii) , essendo ii l'indice del cliclo for. Mi dice, però, che l'indice non può essere maggiore di 3. La mia idea era appunto quella di indicizzare "normalmente" chiedendo al ciclo for di spostarsi colonna dopo colonna sia per il vettore m, sia per il vettore che io ho chiamato Tm, ovvero A^-1*F (So a priori che la matrice sia invertibile ad ogni time-step), utilizzando indici del tipo di quello sopra e anche (:,ii+1). Il software però non gira.
le sei curve, costanti dopo t=1, le avevo ottenute utilizzando semplicemente, per esempio, con (:,1).

Non riesco comunque ad associare la schematizzazione della funzione F che hai fatto tu al mio problema.
Inoltre ho dubbi sull'utilizzo del modulo del vettore m. So che si deve impostare, ma in che modo?

Grazie ancora
Studente35
Starting Member
Starting Member
 
Messaggio: 4 di 12
Iscritto il: 11/12/2018, 10:09

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda feddy » 12/01/2019, 10:33

Ciao,

non dovresti postare immagini. E poi quello che hai scritto nell'immagine sono i passaggi per la ricorrenza del metodo di Heun, che non è nulla di difficile. Piuttosto posta il codice (possibilmente indentato bene).
Non mi è molto chiaro cosa stai dicendo quando dici che hai problemi nell'indicizzazione.

Se riesci a scrivere per bene la $F$, il gioco dovrebbe essere fatto, ma ho come il sospetto che ci sia qualche problema anche lì.
Ho dato un'altra occhiata. La funzione del problema, che tu avevi definito come $F=[-H1+D1;-H2+D2]$, ci sta anche come scrittura. Però non hai specificato le dipendenze da $m_1,m_2$ nei membri che compongono la funzione ! Una scrittura del tipo
Codice:
F=@(m1,m2) [-H1(m1,m2) +D1;-H2(m1,m2)+D2 ]

invece rende la ricorrenza di Heun più facile da scrivere.

Provo a farti vedere come costruire il primo pezzo della funzione (cioè le prime 3 componenti del vettore funzione).

$H_1$ è dato dal prodotto vettoriale tra $m_1$ (che è l'incognita, 3 componenti) e $h_1(m_1,m_2)$, dove, a sua volta, $h_1$ è funzione è $m_1$ e $m_2$.
Quindi basterebbe definire
Codice:
h1=@(m1,m2) [1.3*m1(1)-500*m2(1);-500*m2(2);-500*m2(3)];
Ora $h_1$ è una funzione di $m_1$ e $m_2$.

Quindi risulta che $H_1$ è
Codice:
H1=@(m1,m2) cross(m1,h1(m1,m2));
e come vedi è il primo pezzo della funzione.

Per quanto riguarda $D_1$, questo termine deriva da
Codice:
d1=@(m1) cross(m1,cross(m1,p));


Come vedi, anche qui l'ho scritto come una funzione di $m_1$. In particolare, $D_1$sarà
Codice:
D1=dj.*d1(m1);
(nota come abbia passato $m_1$ a $d_1$ !)

Ora il primo pezzo della $F$ è fatto (a meno della matrice di massa), e le prime 3 righe del vettore $F$ sono date da
Codice:
F=-H1(m1,m2)+dj.*d1(m1)
che è proprio quanto scritto all'inizio di questo esempio. Il secondo pezzo è analogo.
Ultima modifica di feddy il 12/01/2019, 15:49, modificato 1 volta in totale.
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2363 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda feddy » 12/01/2019, 13:30

Avrei due domande:
1. Come fai a sapere, a priori, che quella matrice è non-singolare?
2. Questo problema in che contesto l'hai trovato? E' un esercizio assegnato dal docente? Il meodo di Heun lo usi perché ti è stato detto di usarlo da chi di dovere?
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2364 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda Studente35 » 17/01/2019, 09:15

feddy ha scritto:Avrei due domande:
1. Come fai a sapere, a priori, che quella matrice è non-singolare?
2. Questo problema in che contesto l'hai trovato? E' un esercizio assegnato dal docente? Il meodo di Heun lo usi perché ti è stato detto di usarlo da chi di dovere?


Ciao, continuo a ringraziarti dell'aiuto.
Dunque, per rispondere alle tue domande: questo è un vecchio compito assegnato dal docente che possedeva delle linee guida sulla risoluzione come, appunto, utilizzare A\F, quindi invertendo $ A $ (per questo davo per "scontato" il fatto che fosse non singolare e quindi invertibile assegnando ai vari componenti quei determinati valori, magari mi sbaglio però), e utilizzando appunto il metodo di Heun.
Ho provato a seguire il tuo consiglio,considerando $ m1 $ ed $ m2 $ le variabili, ho dovuto dichiararle prima in questo modo:
m1 = zeros(3,1);
m2 = zeros(3,1);

Poi ho inserito la dipendenza dalle variabili nelle righe che lo necessitavano. L'ho aggiunto anche alla matrice $ A $, è corretto?
e infine, chiamando $ T $ il prodotto di $ A^-1*F $ :
T = @(m1,m2) A\F;

Fin qui tutto corretto?

grazie ancora
Studente35
Starting Member
Starting Member
 
Messaggio: 5 di 12
Iscritto il: 11/12/2018, 10:09

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda feddy » 19/01/2019, 15:20

Sì, mi pare vada bene.
Comunque, se posso dire, non mi pare una grande idea quella di "invertire" la matrice ad ogni passo. E sicuramente tutti i solver di MatLab che trattano questo tipo di equazioni non si lanciano nella risoluzione di un sistema lineare. Va anche detto però che tali solver utilizzano metodi "impliciti". Potresti provare a chiedere al tuo docente spiegazioni in merito.
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2368 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Sistema equazioni differenziali - Integrazione numerica

Messaggioda Studente35 » 20/01/2019, 14:18

feddy ha scritto:Sì, mi pare vada bene.
Comunque, se posso dire, non mi pare una grande idea quella di "invertire" la matrice ad ogni passo. E sicuramente tutti i solver di MatLab che trattano questo tipo di equazioni non si lanciano nella risoluzione di un sistema lineare. Va anche detto però che tali solver utilizzano metodi "impliciti". Potresti provare a chiedere al tuo docente spiegazioni in merito.

va bene. Purtroppo per ora è fuori sede e non ho avuto modo di chiedere delucidazioni in merito. Ma, se posso, in base a quanto mi hai consigliato, cioè di esplicitare il legame con m1 ed m2 ecc, come implementeresti il problema con Heun?

Grazie ancora.
Studente35
Starting Member
Starting Member
 
Messaggio: 6 di 12
Iscritto il: 11/12/2018, 10:09


Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite