Jacobiano funzione che discretizza un problema

[feddy]

Buonasera a tutti,

sto risolvendo un problema ai limiti tramite differenze finite, che è il seguente:

$ { ( u''(x)=1/8(32 + 2x^3 -u(x)u'(x)) ),( u(1)=17 ),( u(3)=43/3 ):}// x \in (1,3) $

Il vero problema non è la risoluzione con questo metodo, bensì scrivere la matrice jacobiana del problema, infatti dopo aver scritto il problema come $F(u)=vec0$, devo trovare $JF(u)$.

Quindi, dette $D2,D1$ le matrici che discretizzano la derivata seconda e la derivata prima rispettivamente,ho:

$F(u)=D2*u -4 -x^3/4 +D1*u^2 - vecb$, con $b$ un vettore che non dipende da $u$

Perciò, direi che $JF(u)=D2 + 2*diag(D1*u) * D1$.

E' corretto?

[dissonance]

Buh, devi chiedere agli analisti numerici. Io, che non ci capisco niente, leggo la formula

$F(u)=D2*u -4 -x^3/4 +D1*u^2 - vecb$

e subito protesto perché non mi trovo con le dimensioni: se \(u\in \mathbb R^n\) allora \(D2 u\in \mathbb R^n\) ma allora cosa significa \(D2u-4\)? Che cosa sarebbe \(x^3\)? Eccetera.

Sospetto che la linearizzazione (=il passaggio alla matrice Jacobiana) vada fatta sul problema continuo, PRIMA di discretizzare. Quindi, senza sapere né leggere né tantomeno scrivere, io credo che per prima cosa tu debba riscrivere il problema *continuo* nella forma
\[
\frac{d}{dx}\mathbf{v} = F(\mathbf{v}(x), x). \]
E questo si fa con la posizione
\[
\mathbf{v}= \begin{bmatrix} u \\ u_x \end{bmatrix}, \quad F(\mathbf v, x)=\begin{bmatrix} v_2 \\ \frac18(32 + 2x^3-v_1v_2)\end{bmatrix}.\]
Spero ti sia utile.

[feddy]

Innanzitutto grazie della risposta

Buh, devi chiedere agli analisti numerici.

E' quello che ho fatto oggi in università . Mi è stato detto che, si potrebbe usare l'usuale regola della catena, cosa che però formalmente è scorretta. Il modo rigoroso, detta $F(u)$ la funzione di cui vogliamo calcolare lo Jacobiano, applicandolo in $v$ si ha $JF(u)*v=lim_(\epsilon -> 0) (F(u+\epsilonv) - F(u))/(\epsilon) $. Per semplicità di notazione, $v,u$ sono da intendersi come vettori in $\mathbb{R^n}$.

subito protesto perché non mi trovo con le dimensioni: se $u \in RR^n$ allora $D2u \in RR^n$ ma allora cosa significa $D2u−4$? Che cosa sarebbe $x^3$? Eccetera

Con $D2*u -4$: il prodotto matrice vettore mi restituisce un vettore in $RR^n$ a cui sottraggo il valore $4$ in ogni componente. La sintassi MatLab lo accetta già così, ma avrei dovuto specificare.

Sospetto che la linearizzazione (=il passaggio alla matrice Jacobiana) vada fatta sul problema continuo, PRIMA di discretizzare. Quindi, senza sapere né leggere né tantomeno scrivere, io credo che per prima cosa tu debba riscrivere il problema *continuo*

Quello che pensavo anche io, ma mi è stato detto che prima bisogna discretizzare, e poi differenziare.


Per quanto riguarda la mia $F(u)$, quella riportata nel primo post è errata (quell' $u^2$ non ha senso alcuno), ed è $F(u)=D_2 * \vecu - 4 - x^3/3 + (D_1 * \vecu )*(I*\vecu) - \vec{b}$.

Ora non ho il tempo per provare lo jacobiano con la definizione, appena avrò un momento stasera proverò e posterò ciò che mi risulta

[Raptorista]

Moderatore: Raptorista

Sposto da Analisi Superiore.

[feddy]

Come detto, mi sono messo giù per risolvere il limite che ho scritto nell'ultimo post, e mi sto trovando in seria difficoltà.

Questo diventa $ lim_(\epsilon -> 0) ((D(u+\epsilonv).*(u+\epsilonv)) - (Du).*u)/epsilon $ , e già qui ho difficoltà perchè non capisco il senso di moltiplicare un vettore per un altro vettore. So che "devono" venire fuori delle matrici, e pertanto da qualche parte dovrebbe venire fuori un $diag(...)$, ma non capisco veramente il senso.

Qualcuno di voi ha mai avuto a che fare con questo?

[dissonance]

No, e capisco le tue difficoltà perché sono le stesse obiezioni che dicevo io nel mio post precedente. Tanto per cominciare, metti dei vettori \(\mathbf{1}\) (questo simbolo indica il vettore con tutte le componenti uguali a \(1\)) nella definizione di \(F\). Per esempio invece di \(x^3\) credo dovresti scrivere \(x^3\mathbf{1}\). Ma il vero problema è il termine \(u(x)u'(x)\) che diventa una roba senza senso: tu lo hai scritto come prodotto scalare, ma allora è uno scalare?!? Non può essere, deve essere un vettore.

Secondo me il "giusto" prodotto è questo qui:
\[
u\bullet v = \begin{bmatrix} u_1v_1 \\ u_2 v_2 \\ \vdots \\ u_n v_n \end{bmatrix}. \]
Forse è utile osservare che
\[
u\bullet v = \mathrm{diag}(u\otimes v), \]
dove \(\otimes\) denota il prodotto tensoriale, che in questo caso è il prodotto colonna per riga:
\[
u\otimes v = [u_i v_j]_{1\le i, j\le n}.\]

[Raptorista]

Provo a rassicurare il buon dissonance
\(u\) è una funzione scalare di una variabile scalare, quindi \(u \cdot u'\) è uno scalare. Tutti questi vettori saltano fuori perché nel metodo delle differenze finite si scelgono \(n\) punti \(x^i\) appartenenti al dominio di \(u\), che in questo caso è \([1,3]\) e si definiscono i vettori \(u_h \in \mathbb R^n\) di elementi \(u_h^i = u_h(x^i)\).
A questo punto si scrive l'equazione differenziale in ogni punto \(x^i\) sostituendo le derivate continue con dei rapporti incrementali tipo
\[
u''(x^i) \approx \frac{u_h^{i+1}-2u_h^i+u_h^{i-1}}{x^{i+1}-x^{i-1}}
\]
e si ottiene un sistema nonlineare di dimensione \(n\) che qui è stato indicato come \(F(u) = 0\). Per rendere la vita più difficile agli studenti ci siamo messi tutti d'accordo di introdurre complicate matrici \(D\) che altro non fanno che mandare il vettore \(u_h\) nel vettore le cui componenti sono date dalla formula qui sopra [o analoghe per altre derivate], tenendo conto anche delle condizioni al bordo.

[feddy]

Grazie @dissonance per l'osservazione che hai fatto, è quello che cercavo e che ho visto fare a alezione, anche se nessuno ha parlato di tensori

Provo a risolvere sto benedetto limite, fattore per fattore: [NOTA: uso $"@"$ per indicare quello che tu scrivi \bullet, non so perché non lo legga correttamente]

1.
Il "pezzo" $(Du) \@ u = diag((Du)u)$

2.
Il pezzo $D(u+\epsilon v) \@ (u+\epsilon v)=(Du) \@ u + \epsilon (Du) \@ v + \epsilon (Dv) \@ u + \epsilon^2(Dv) \@ v$

Facendo il limite, mi risulta che il termine con $\epsilon^2$ va a $0$, mentre, raccogliendo $\epsilon$, il limite vale $(Du) \@ v + (Dv) \@ u= diag((Du) \@ v) + diag((Dv) \@ u)=2 diag((Du) \@ v)$.

Da cui dovrebbe seguire che $JF(u)=2 *diag(Du)$, che sarebbe una matrice diagonale che ha sulla diagonale le componenti del vettore $D \vec(u)$.

Non sono certo però della correttezza del risultato.

[feddy]

Ci tenevo a informarvi che ho risolto in due modi diversi: (il primo da solo, il secondo è quello esposto sopra ma risolto con l'aiuto del professore).

Modo 1.

Scrivo la generica riga i-esima del problema, che risulta essere $(u_{i+1} - u_{i-1}) / (2h) * u_i$ e derivo componente per componente, ottenendo così una matrice tridiagonale. Per situazioni più complesse però diventa alquanto orrendo (specie se devo inserire manualmente le componenti in Octave o MatLab).

Modo 2.

Scrivo in sintassi MatLab, che risulta essere più intuitivo (almeno al momento) per me, visto che di tensori non ne ho finora mai sentito parlare

Sviluppando il limite arrivo all'epressione

$(Du).*v + (Dv).*u = diag(Du).*v + diag(u)*Dv$

Pertanto, essendo l'ultima espressione pari a $JF(u)*v$, allora $JF(u)=diag(Du) + diag(u)*D$.

Ora dovrebbe essere tutto apposto