Equazione differenziale con metodo di Newton

[feddy]

Sui metodi numerici per le PDE c'è un macello di roba, e dipende sempre dal tipo di equazione che si ha davanti (paraboliche, iperboliche,.ecc.). Io sono al terzo anno di Matematica e quello che ho visto su metodi numerici per alcuni tipi di PDE è ancora poco, dovresti chiedere al tuo relatore o magari aspettare l'intervento di qualcuno più ferrato in materia che passi di qui.

O magari, chiedere su computational science

[luca.milano]

Nel frattempo rimaniamo sulle ODE, ho provato a scrivere questo codice che sostanzialmente risolve la stessa equazione di prima ma con $ T(0) $ e $ T(L) $ che possono essere posti a piacimento, il problema è che mi da errore e non capisco bene cosa sto sbagliando. Questo è il codice (è molto simile al tuo):

Codice:

L=1; %lunghezza
m=50; %numero nodi
h=L/m; %intervallo spaziale

T0=0;  %temperatura in x=0
TL=0;  %temperatura in x=L

g=50;  %termine di generazione
k0=10; %primo termine della polinomiale per la conduttività termica: k=k0+aT
a=200;  %coefficiente della T
save indiretto

%MATRICE CHE DISCRETIZZA LA DERIVATA SECONDA
A=toeplitz(sparse([1,1],[1,2],[-2,1],1,m));
I=speye(m);

%BOUNDARY CONDITIONS
A(1,1:2)=[1,0]; %temperatura ingresso nota
A(m,m-1:m)=[0,1]; %temperatura uscita nota

b=@(u) sparse([1,m],[1,1],[g/(k0+a*u)+T0 , g/(k0+a*u)+TL],m,1);
 F=@(u) A*u + g./(k0+a*I*u)-b(u);
JF=@(u) A - g*a*spdiags([0;1./((k0+a*u(2:m-1)).^2);0],0,m,m);

%Newton:

tol=h^2; %diff finite secondo ordine
u0=ones(m,1);
res=-JF(u0)\F(u0);
iter=0;
maxit=50;
while (norm(res,inf)>tol) & (iter<maxit)
   iter=iter+1;
   u0=u0+res;
   res=-JF(u0)\F(u0);
end
u=u0+res;

x=linspace(0,L,m);
plot(x,u0,'b-o')

L'errore è nella funzione F: dice che devono avere la stessa dimensione ma a me sembra proprio che abbiano la stessa dimensione. Te come lo modificheresti?

EDIT
mi ero accorto di un errore e ho modificato il codice. Purtroppo ancora non gira.

[feddy]

Più che ODE, questi sono BVP (boundary value problems / problemi ai limiti).

Non mi sembra sensato definire il termine noto come funzione di $u$, non bastava cambiare solo la prima e ultima componente, i.e. $b(0) , b(m)$, dove nelle rispettive espressioni compare $T(0)$ e $T(L)$ ?

[luca.milano]

Giusto, BVP
Effettivamente si fa prima a modificare b, a me l'unica soluzione che viene in mente è questa:

Codice:

b=zeros(m,1);
b(1)=-g/(k0+a*T0)-T0;
b(m)=-g/(k0+a*TL)-TL;

F=@(u) A*u + g./(k0+a*I*u) + b;


Eppure non funziona

[feddy]

Con il mio codice, e la tua modifica su $b(1),b(m)$ appena riportata, le condizioni al bordo vengono corrette. Ho provato più volte, con diversi valori. Prova a ricopiare il mio codice e a mettere il tuo vettore $b$, funziona.

[luca.milano]

Se invece volessi risolvere un problema con condizione di Robin di questo tipo:
$ { ( (partial^2 u)/(partial x^2) +g/(k(u))=0),( -k(u)*(partial u)/(partial x) (0)= h(T-T0) ),( T(L)=T_L ):} $

In cui $ T_0 $, $ T_L $, $ h $ e $ g $ sono dei termini noti, mentre: $ k(u)= k0+k1*u+k2*u^2 $
Chiaramente $ k0 $, $ k1 $ e $ k2 $ sono termini noti.

Ci sto riflettendo ma non saprei proprio come fare (senza dover ricorrere a dei cicli per creare le matrici F e JF)!

[feddy]

La condizione di Robin è con $-k(u(0))$, giusto? Intendo, si sa quanto vale $k(u(0))$? Bisogna scrivere la prima riga discretizzata e utilizzare un "nodo fantasma" al bordo $x=0$ per la derivata prima. Tutto quello che alla fine risulta cambiato è solamente la prima riga.
Ora non ho tempo di mettermici, ti conviene guardare su qualche testo o dispensa in rete per problemi ai limiti.

[luca.milano]

Si ci hai preso, era $ k(u(0)) $ e non si sa quanto vale proprio perchè non si sa quanto vale $ u(0) $ . E' proprio la classica condizione di Robin. Ho guardato qualcosa in rete, ma si trovano solo casi lineari e io sto avendo qualche problema con una parte di codice da te riportata, che funziona, ma non la capisco:

Codice:

F=@(u) A*u + g./(k0+a*I*u) + b;

Alla matrice $ A*u $ ci sommi prima il termine diagonale $ (g)/(k0+a*I*u) $ (è diagonale, se ho capito bene, perché così vado a sommare solo i termine del nodo $ i-esimo $ ). Poi però ci sommi $ b $ e qui mi perdo perché mi sarei aspettato anche qui un termine di tipo diagonale, in modo da andare a toccare solo i termini $ (i,i) $ della matrice che stiamo costruendo.

Per ultimo, nel momento in cui viene dato questo comando

Codice:

u0=ones(m,1);
res=-JF(u0)\F(u0);

$ F(u0) $ si trasforma inspiegabilmente (per me) in un vettore.

Io mi sarei aspettato di creare $ F $ direttamente come vettore così:

Codice:

F=@(u) A*ones(m,1)*u

per poi aggiungerci il termine non lineare e il termine correttivo, sempre pensati come vettori colonna.

Scusa se ti ho fatto un po' di domande tutte in una volta, abbi pazienza

[feddy]

Strano che su internet non si trovi nulla... se ho un po' di tempo in questi giorni mi ci metto, però l'approccio che ho scritto sopra non dovrebbe cambiare.

Rispondo ora alle tue domande:
Innazitutto:

Alla matrice A*u [...]

non è una matrice , è un vettore, e infatti il termine $(g)/(k0+a*I*u)$ è un vettore . Segue dunque, che $b$ è la correzione apportata alla prima e ultima riga di $F$. Questa $F$, per inciso, è e deve essere una funzione da $RR^m \rarr RR^m$ e quindi ovviamente $F(u)$ un vettore colonna.
Questo modo di procedere, per apportare la "correzione" alle righe destinate alle condizioni ai bordi non mi piace moltissimo, ma è dovuto al fatto che MatLab non permette una scrittura diversa della $F$, cosa che invece la sua alternativa free Octave fa.

Avrei potuto scrivere, e ciò sarebbe stato chiarissimo:

Codice:

F=@(u) [u(1)-T0; (A*u + g./(k0+a*I*u))(2:m-1);u(m)-TL];

E' evidente che la soluzione $\bar(u)$ di $F(u)=0$, calcolata mediante Newton, sarà tale che $u(1) - T0=0$, da cui $u(1)=T0$. Analogamente per $u(m)$ si vede che le condizioni ai bordi sono verificate.

Io mi sarei aspettato di creare F direttamente come vettore così:

Codice:

 F=@(u) A*ones(m,1)*u

In tal caso la soluzione sarebbe $\vec{u}=\vec{0}$... mancano tutti gli altri termini.

[luca.milano]

Ok, ora è tutto molto più chiaro e il codice è partito e ora funziona bene.
Ho provato a cimentarmi con un problema simile:
$ { ( (d^2 u)/dx^2 +g/k=0),( u(0)=u_0 ),( (d u)/dx(0)=-q_0/(k(u_0)) ):} $

$ q_0 $ è una costante (fisicamente è il flusso termico in ingresso del corpo). $q(x)=-k*(du)/(x)$ dove $u$ è la temperatura (ma questo ci interessa meno).

Tutti i termini noti e il $ k $ polinomiale sono sempre gli stessi come li ho definiti nelle risposte precedenti. Per quella questione del nodo fantasma ho usato un approccio diverso ma anche usando la definizione di nodo fantasma trovata in internet arrivo allo stesso risultato. Tramite lo sviluppo in serie di Taylor al secondo ordine posso scrivere il nodo $2$ in questo modo:

$ u_2 = u_1+(du)/(dx)(0)*Delta x +1/2*(d^2u)/(dx^2)(0)*Delta x^2 $

Sostituendo:

$ { ( (du)/(dx)(0)=-(q(0))/(k(u0)) ),( (d^2u)/(dx^2)(0)=-g/(k(u_0)) ):} $

La prima delle due viene dalle BC, mentre la seconda viene dall'equazione differenziale. Fatto questo si può scrivere per il nodo $2$:

$ u_2=u_0-q_0/(k(u(0)))*Delta x-1/2*g/(k(u_0))*Delta x^2 $

Ho quindi trovato una relazione che mi permette di calcolare $u_2$. Per quanto riguarda il codice ho tenuto il tuo schema, facendo qualche piccola modifica sulle condizioni al contorno, che ovviamente dovranno essere imposte sui nodi $ 1 $ e $ 2 $ invece che sui nodi $ 1 $ ed $ m $.
Questo è il codice:

Codice:

%CODICE LASTRA PIANA 1D, GENERAZIONE,
%TEMPERATURA IMPOSTA IN INGRESSO E FLUSSO IMPOSTO IN INGRESSO
%EQUAZIONE: d2T/dx2 +g/k=0 con k=k0+k1*T
clear
clc

L=1;   %lunghezza
m=50; %numero nodi
h=L/m; %intervallo spaziale

T0=5; %temperatura in x=0
q0=0;   %flusso termico in x=0

g=0;  %termine di generazione
k0=100; %primo termine della polinomiale per la conduttività termica: k=k0+aT
k1=20;  %coefficiente della T

%MATRICE CHE DISCRETIZZA LA DERIVATA SECONDA
A=toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
I=speye(m);

%BOUNDARY CONDITIONS
A(1,1:2)=[1 0];
A(2,1:3)=[-1 1 0];
I(1,1)=0; I(2,2)=0;

b=zeros(m,1);
b(1)= -g/k0 - T0;
b(2)= -g/k0 + q0/(k0+k1*T0)*h + 1/2*g/(k0+k1*T0)*h^2;

F=@(u) A*u + g./(k0+k1.*I*u) + b;

JF=@(u) A-k1.*g.*spdiags([0;0;1./((k0+k1.*u(3:m))).^2],0,m,m);

%Newton:

tol=(h^2); %diff finite secondo ordine
u0=ones(m,1);
res=-JF(u0)\F(u0);
iter=0;
maxit=50;
while (norm(res,inf)>tol) & (iter<maxit)
   iter=iter+1;
   u0=u0+res;
   res=-JF(u0)\F(u0);
end
u=u0+res;

x=linspace(0,L,m);
plot(x,u0)



A me sembra uguale identico a quello scritto con le condizioni di prima (che funziona), mentre questo gira (mi sembra che calcoli la funzione al nodo $ 2 $ in modo corretto) ma il resto del calcolo è sicuramente sbagliato.

Così di primo impatto mi verrebbe da pensare che c'è qualche errore nella funzione o nello Jacobiano, ma non vedo davvero niente di strano (ed è tutto il giorno che ci lavoro ). Te noti niente di sbagliato?