Problema differenziale (standard)

Messaggioda feddy » 13/11/2017, 17:54

Ciao a tutti, sono alle prese col seguente problema differenziale:

$ { ( y''(x) + (y'(x))/x =cos(y(x))),( y'(0)=0 ),( y(0)=1 ):} , x \in (0,1]$

Lo risolvo con differenze finite. Poiché ho condizioni di Neumann sul primo nodo e sempre sullo stesso ho una di Dirichlet, ho pensato di porre $u_1=1$, in quanto $y(0)=1$, e poi utilizzare un ghost node $u_0$ per discretizzare la derivata prima.
In sostanza, detto $h$ il passo di discretizzazione, ho $y'(0)=(u_2-u_0) / (2h) = 0$, da cui $u_0=u_2$.

Sostituendo nell'equazione, correggetemi se sbaglio, trovo che la prima riga del sistema deve risultare uguale a
$2/h^2 u_2 - 2/h^2 -cos(1)=0$
.

A tal fine, detta $D_1$ la matrice di discretizzazione della derivata prima, $D_2$ quella della seconda e $I$ l'identità, ho che il problema discretizzato equivale a risolvere l'equazione $F(u)= D_2 * u + 1/x * (D_1*u) - cos(u) =0$, dove in realtà, per non rendere pesante la notazione, qui $u=vecu$.

Evidentemente, c'è un problema con quel zero al denominatore: per aggirarlo ho considerato un linspace che parte da $0+eps$ e va fino a $1$.

Allego il codice del problema (in MatLab)

Codice:
m=100;
h=(1-eps)/(m-1);
x=linspace(0+eps,1,m)';

%matrici di discretizzazione + boundary conditions

D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
I=speye(m);


D2(1,1:2)=[0,2/(h^2)];
D1(1,2)=0;
I(1,1)=0;

b=zeros(m,1);
b(1,1)= 1 -2/h^2 - cos(1);

%funzione

Fun=@(u) D2*u + spdiags(1./x,0,m,m)*(D1*u) - cos(I*u) + b;

Jac=@(u) D2 + spdiags(1./x,0,m,m)*D1 + diag(sin(I*u));

%newton
u0=ones(m,1); %guess iniziale
tol=1e-10;

maxit=200;
u=u0;
delta=-Jac(u)\Fun(u);
i=0;

%Newton per sistemi NON-LINEARI
while(norm(delta,inf)>tol && i<maxit)
    i=i+1;
    u=u+delta;
    delta=-Jac(u) \ Fun(u);
end
u=u+delta;


plot(x,u,'*')


Il plot che ne risulta è il seguente: dalle condizioni al contorno mi aspetterei tangente orizzontale, invece qui è proprio il contrario. Sbaglio a applicare differenze finite oppure c'è un errore nell'implementazione?

Immagine


Grazie a chiunque mi possa illuminare :-D
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1557 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Problema differenziale (standard)

Messaggioda Raptorista » 14/11/2017, 09:54

feddy ha scritto:Poiché ho condizioni di Neumann sul primo nodo e sempre sullo stesso ho una di Dirichlet,

Per questo tipo di problemi si chiamano condizioni di Cauchy.
feddy ha scritto: ho pensato di porre $u_1=1$, in quanto $y(0)=1$, e poi utilizzare un ghost node $u_0$ per discretizzare la derivata prima.
In sostanza, detto $h$ il passo di discretizzazione, ho $y'(0)=(u_2-u_0) / (2h) = 0$, da cui $u_0=u_2$.

E \(u_0\) poi a cosa è uguale? Non mi sembra risolto così il problema.

Non esiste un solo schema a differenze finite, anzi ne esistono una miriade. Comincia a scrivere chiaramente quale stai usando e come discretizza le derivate prima e seconda.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 4683 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Problema differenziale (standard)

Messaggioda feddy » 14/11/2017, 10:48

Ciao,
lo schema alle differenze finite che sto usando è quello che discretizza la derivata prima e la derivata seconda mediante le matrici $D1$ e $D2$, rispettivamente scritte come sotto ( $m$ è il numero di nodi, mentre $h$ è il passo della discretizzazione )

Esplicitamente, la derivata prima la discretizzo così $(u_{i+1}- u_{i-1})/(2h $, mentre la derivata seconda con $(u_{i+1} - 2u_i + u_{i-1})/h^2$


Codice:
D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));




Detto questo, passiamo al problema:
Visto che lo $0$ è escluso dalle $x$ del problema, ho considerato un linspace che va da $0+ \epsilon$ fino a $1$, con $epsilon =eps$


$u_0$ è un nodo fantasma, o ghost node, e serve per discretizzare la derivata prima nell'estremo sinistro.
In pratica, per scrivere la condizione sul primo nodo $y'(0)=0$, dicretizzo $y'(0)$ scrivendo
$(u_2-u_0)/(2h)= 0$
da cui ricavo $u_2=u_0$.

Grazie a questa condizione, scrivo la prima riga del problema, cioè:
$(u_0 - 2u_1 + u_2)/h^2 + 1/x*(u_2- u_0)/2h - cos(u_1)=0$


Sfruttando $u_2=u_0$,e il fatto che l'altra condizione mi dice che $u_1=1$, la prima riga del sistema deve risultare uguale a:
$2/h^2 u_2 -2/h^2 -cos(1)=0$


Affinché questo accada, ho dovuto correggere l'equazione con un termine noto $b$, che nel codice che ho allegato è un vettore di zeri che in prima posizione ha il valore $b(1,1)= 1 -2/h^2 - cos(1)$.


La funzione da azzerare è perciò
Codice:
F(u)= D2*u + spdiags(1./x,0,m,m)*(D1*u) - cos(I*u) + b
(scritto in sintassi MatLab).
Evidentemente è un problema non-lineare e va risolto, per esempio, con Newton.
Lo jacobiano risulta essere
Codice:
Jac(u)= D2 + spdiags(1./x,0,m,m)*D1 + diag(sin(I*u))


Poi ho applicato Newton e ho plottato la soluzione, che come vedi non pare proprio essere quella giusta. Mi aspetterei tangente orizzontale in $x=0$, mentre sembra quasi che ci sia un asintoto dal mio plot. Evidentemente è perché c'è quella $x=0$, e non so come andare avanti.


Spero di essere stato abbastanza chiaro
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1558 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Problema differenziale (standard)

Messaggioda Raptorista » 14/11/2017, 12:35

Se tu aggiungi una variabile \(u_0\) ed una equazione \(u_0 = u_2\) hai solo ingrandito il sistema di una riga e una colonna, ma non hai risolto il tuo problema. Io credo che tu debba assegnare \(u_0 = u_1 = 1\) visto che la derivata prima è nulla.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 4686 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Problema differenziale (standard)

Messaggioda feddy » 14/11/2017, 20:26

Scusa ma non mi pare di aver aggiunto nulla. Infatti la prima riga del sistema non dipende da $u_0$. Ad ogni modo, ho posto $u_1=1$. Questo viene dalla condizione $u_1=1$. Sono abbastanza certo di questo procedimento
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1560 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Problema differenziale (standard)

Messaggioda Raptorista » 14/11/2017, 22:32

Negli episodi precedenti...
feddy ha scritto:ho pensato di porre $u_1=1$, in quanto $y(0)=1$, e poi utilizzare un ghost node $u_0$ per discretizzare la derivata prima.


And now...
feddy ha scritto:Scusa ma non mi pare di aver aggiunto nulla.


Immagine


Comunque,
feddy ha scritto:Ad ogni modo, ho posto $ u_1=1 $. Questo viene dalla condizione $ u_1=1 $. Sono abbastanza certo di questo procedimento

Non fa una piega.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 4689 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Problema differenziale (standard)

Messaggioda feddy » 14/11/2017, 22:58

Raptorista ha scritto:Comunque,

feddy ha scritto: Ad ogni modo, ho posto u1=1. Questo viene dalla condizione u1=1. Sono abbastanza certo di questo procedimento



Non fa una piega.
Scusa ma ho scritto questo post dopo una giornata intera di studio e evidentemente sono fuso. :lol:

Ricapitolo:
-condizione $y(0)=1$, segue subito $u_1=1$.
Sul "ghost node" la tecnica che ho descritto l'abbiamo sempre usata per problemi ai limiti, e non ho mai avuto alcun tipo di problema. Provo a spiegarmi meglio: dalla condizione sulla derivata prima $y'(0)=0$, trovo che $(u_2 -u_0)/(2h)=0$, da cui l'uguaglianza $u_0=u_2$, e scrivendo la prima riga del problema ho che questa riga deve essere
$2/h^2 u_2 -2/h^2 -cos(1)=0$


Ora si tratta solamente di forzare il sistema ad avere la prima riga fatta proprio come scritto sopra, e lo faccio aggiungendo un termine noto. Su questo non credo ci siano dubbi, è una tecnica che ho sempre visto usare e che mi pare sensata.

Non riesco, sinceramente, a capire quello che mi contesti. Ti chiedo scusa
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1561 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Problema differenziale (standard)

Messaggioda Raptorista » 15/11/2017, 10:03

Non sto contestando niente XD
Converrai però che il nodo fantasma è una cosa che si aggiunge a posteriori, non è parte del problema iniziale. Quando poi dici che "segue subito \(u_1 = 1\)" credo che ci sia un piccolo imbroglio: è vero che \(u_1 = 1\) ma sospetto che il procedimento corretto è di mettere uno o più nodi fantasma a sinistra [il numero dipende dal grado di differenze finite in uso] e poi usare lo schema per dedurre quelli di destra. In questo caso lo schema è centrato di ordine 2 e quindi il risultato è comunque giusto, ma non so se la cosa succederebbe anche con altri schemi.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 4692 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Problema differenziale (standard)

Messaggioda feddy » 23/11/2017, 15:48

Mi sono preso del tempo per pensarci, e finalmente ora mi torna tutto ! Grazie Raptorista :)
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 1582 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Problema differenziale (standard)

Messaggioda Raptorista » 24/11/2017, 12:05

Good to know! Buona fortuna!
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 4705 di 9616
Iscritto il: 28/09/2008, 19:58


Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite