Equazione differenziale con metodo di Newton
Buongiorno a tutti, sto cercando di risolvere un'equazione differenziale non lineare discretizzando in N nodi e poi risolvendo il sistema risultante con il metodo di Newton. Il codice mi sembra che giri ma è abbastanza lento e vorrei capire se c'è qualche errore o se semplicemente questi sono i limiti del metodo di Newton.
Questa è l'equazione:
$ { ( (partial^2 T)/(partial x^2) + g/k=0 ),( T(0)=T_0 ),( T(L)=T_L ):} $
dove il termine $ k $ è pari a:
$ k=k0+a*T $
Vi riporto lo script principale:
La function per il vettore delle equazioni è la seguente:
Lo Jacobiano è invece questo:
Vedete qualche errore? Ci sono metodi migliori per risolvere sistemi non lineari?
Questa è l'equazione:
$ { ( (partial^2 T)/(partial x^2) + g/k=0 ),( T(0)=T_0 ),( T(L)=T_L ):} $
dove il termine $ k $ è pari a:
$ k=k0+a*T $
Vi riporto lo script principale:
L=1; %lunghezza
N=50; %numero nodi
dx=L/N; %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=10; %coefficiente della T
save indiretto
%METODO DI NEWTON
T=zeros(N,1);
x1=ones(N,1);
n=0;
while norm(T-x1)>1e-5
n=n+1;
x1=T;
r=jac2(T)\sist2(T);
T=T-r;
end
T;
n;La function per il vettore delle equazioni è la seguente:
function F=sist2(T)
%devo creare il vettore colonna del sistema: F=0
load indiretto
F(1)=T(1)-T0;
F(N)=T(L)-TL;
for j=2:N-1
k(j)=k0+(10^-1)*T(j);
F(j)=T(j+1)-2*T(j)+T(j-1)+g/k(j)*dx^2;
end
F=F';
Lo Jacobiano è invece questo:
function J=jac2(T)
load indiretto
J=zeros(N,N);
J(1,1)=1; J(N,N)=1;
for j=2:N-1
k(j)=k0+10^-1*T(j);
J(j,j-1)=1;
J(j,j)=-2-a*g/(k0+a*T(j))^2;
J(j,j+1)=1;
endVedete qualche errore? Ci sono metodi migliori per risolvere sistemi non lineari?
Risposte
Ciao,
deduco che tu stia usando uno schema alle differenze finite equispaziate. E' un semplice problema ai limiti con condizioni al contorno di Dirichlet.
Probabilmente evitare di scrivere due diverse funzioni solo per scrivere la $F$ e lo jacobiano (con cicli for !) è meglio. Questo perché MatLab è "ottimizzato" per questo genere di operazioni, vedi anche qui.
Per quel che mi riguarda, preferisco sempre discretizzare gli operatori differenziali, in questo caso $\frac{d}{dx^2}$, in modo "matriciale", ossia creando una matrice triangolare, ma non cambia nulla a livello concettuale.
Il tuo metodo di Newton mi pare corretto, anche se le dichiarazione $x_1=T$ non vedo che senso possa avere.
deduco che tu stia usando uno schema alle differenze finite equispaziate. E' un semplice problema ai limiti con condizioni al contorno di Dirichlet.
Probabilmente evitare di scrivere due diverse funzioni solo per scrivere la $F$ e lo jacobiano (con cicli for !) è meglio. Questo perché MatLab è "ottimizzato" per questo genere di operazioni, vedi anche qui.
Per quel che mi riguarda, preferisco sempre discretizzare gli operatori differenziali, in questo caso $\frac{d}{dx^2}$, in modo "matriciale", ossia creando una matrice triangolare, ma non cambia nulla a livello concettuale.
Il tuo metodo di Newton mi pare corretto, anche se le dichiarazione $x_1=T$ non vedo che senso possa avere.
Si hai ragione, è lineare, è che il mio problema è un altro e volevo provare una versione più semplificata per fare pratica 
.
Il problema completo sarebbe questo:
$ k(T)*(partial^2 T)/(partial x^2) + a(T)*(partial T)/(partial x) + g=0 $
Poi a sua volta questa equazione devo inserirla in un sistema di 5 equazioni non lineari (se possibile in geometria bidimensionale, altrimenti la tengo monodimensionale).
Il problema è che se già la versione semplificata che ho proposto gira con fatica, posso solo immaginare quanto potrà metterci a risolvere il problema completo.
.Il problema completo sarebbe questo:
$ k(T)*(partial^2 T)/(partial x^2) + a(T)*(partial T)/(partial x) + g=0 $
Poi a sua volta questa equazione devo inserirla in un sistema di 5 equazioni non lineari (se possibile in geometria bidimensionale, altrimenti la tengo monodimensionale).
Il problema è che se già la versione semplificata che ho proposto gira con fatica, posso solo immaginare quanto potrà metterci a risolvere il problema completo.
N.B.: codice editato.
L'output prodotto è:
- - - - - - - - - - - - -
Chiaramente è un approccio diverso dal tuo, dunque subito potrebbe non essere chiaro.
Un paio di osservazioni sul mio codice:
clear all close all L=1; %lunghezza m=50; %numero nodi h=L/m;%passo spaziale x=linspace(0,L,m); %Dati al bordo + costanti del problema T0=0;TL=0; g=50;k0=10;a=10; %Matrice che discretizza 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(m,m-1:m)=[0,1]; I(1,1)=0; I(m,m)=0; b=zeros(m,1); b(1)=-g/k0; b(m)=-g/k0; F=@(u) A*u + g./(k0+a*I*u) + b; 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++; u0+=res; res=-JF(u0)\F(u0); end u0+=res; plot(x,u0,'b-o')
L'output prodotto è:
- - - - - - - - - - - - -
Chiaramente è un approccio diverso dal tuo, dunque subito potrebbe non essere chiaro.
Un paio di osservazioni sul mio codice:
- 1. E' scritto in Octave, quindi se usi MatLab incrementi del tipo i++, u0+=res vanno scritti rispettivamente i=i+1 e u0=u0 +res
2. Ho scritto la funzione $F$, a valori vettoriali, da azzerare. Ho scelto di aggiunere un vettore "correttivo" del termine noto per imporre le condizioni ai bordi per cercare di fare il più simile possibile a te.
3. Ora ho dovuto effettivamente utilizzare Newton, quindi ho calcolato la jacobiano $JF$, identico al tuo [/list:u:2apwfbvz]
Editato
Beh, se $k$ e $a$ sono funzioni di $T$ allora anche il secondo problema da te proposto è ancora lineare...
Ti ringrazio per la risposta, sei stato gentilissimo. Innanzitutto mi hai fatto vedere un metodo per creare le matrici più agilmente e cercherò di studiarlo e impararlo .
L'unico problema riguardo il tuo codice è che la soluzione non dipende dal parametro $ a $ ! Se provo a cambiarlo rimane tutto uguale. E' possibile che tu abbia dimenticato un pezzo?
.L'unico problema riguardo il tuo codice è che la soluzione non dipende dal parametro $ a $ ! Se provo a cambiarlo rimane tutto uguale. E' possibile che tu abbia dimenticato un pezzo?
Sì hai ragione, leggendo in fretta non avevo visto anche $k0$... ho editato il messaggio sopra
Grazie ancora! Mi confermi che posso estendere il metodo anche a sistemi di equazioni differenziali alle derivate parziali? Le soluzioni che trovo sono accurate?
Eh questo credo sia un po' troppo... in che contesto ti è richiesto di risolvere questo?
Tesi di laurea magistrale in ingegneria. Altrimenti che strumenti bisogna adoperare?
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
O magari, chiedere su computational science
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):
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.
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.
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)$ ?
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)$ ?
Giusto, BVP 
Effettivamente si fa prima a modificare b, a me l'unica soluzione che viene in mente è questa:
Eppure non funziona
Effettivamente si fa prima a modificare b, a me l'unica soluzione che viene in mente è questa:
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
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.
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)!
$ { ( (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)!
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.
Ora non ho tempo di mettermici, ti conviene guardare su qualche testo o dispensa in rete per problemi ai limiti.
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:
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
$ F(u0) $ si trasforma inspiegabilmente (per me) in un vettore.
Io mi sarei aspettato di creare $ F $ direttamente come vettore così:
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
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
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ì:
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
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:
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:
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.
In tal caso la soluzione sarebbe $\vec{u}=\vec{0}$... mancano tutti gli altri termini.
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:
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ì:
F=@(u) A*ones(m,1)*u
In tal caso la soluzione sarebbe $\vec{u}=\vec{0}$... mancano tutti gli altri termini.
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:
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?
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 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?
Sbagli dal principio, quella è una ODE, hai solo i dati a $x=0$, non hai nessuna condizione sull'estremo destro. Non puoi discretizzare l'intervallo se un intervallo non ce l'hai. Questa è roba base base, ti conviene prenderti ([size=55]scaricare[/size]) un libro sulla risoluzione numerica di equazioni differenziali. O magari qualche dispensa tipo questa, ma è solo un esempio.
Devi ricondurti a un sistema del prim'ordine introducendo delle variabili ausiliarie e poi applicare uno schema numerico opportuno (Eulero implicito, trapezi, Runge Kutta). Se non vuoi scriverti il codice del metodo basta prendere MatLab e usare odexyz e trovarti la soluzione.
Quello che segue era la mia prima risposta a questo problema. Non mi ero accorto che fosse una ODE e mi sono concentrato solo sul risponderti sulla condizione di Robin del BVP.
[size=150]EDIT[/size]
Ora per me sono giorni un po' pesanti (come lo sono per te immagino), ma appena ho un po' di tempo do' meglio una occhiata al tuo codice.
Per la condizione di Robin a $x=0$ provo a dirti come farei io, senza passare dallo sviluppo di Taylor, per comodità:
Questa è $\frac{d}{dx}u(0) + \frac{q_0}{k_0 + k_1 u+k_2 u^2}=0$
Discretizzo questa riga
$\frac{u_0 -2u_1 +u_2}{h^2} + \frac{g}{k} =0$, ($\square$) dove con $u_0$ indico il valore della soluzione nel nodo fantasma, mentre $h$ è il passo spaziale. Chiaramente il problema ora è la ricerca del valore della soluzione al nodo fantasma, e questo viene ricavato dalla condizione di Robin stessa, che ora scrivo per esteso:
Ora da qui ricavi $u_0$, lo sostituisci in $\square$ e ricavi l'espressione della prima riga, che poi modificherai di conseguenza
[size=150]Fine EDIT[/size]
Devi ricondurti a un sistema del prim'ordine introducendo delle variabili ausiliarie e poi applicare uno schema numerico opportuno (Eulero implicito, trapezi, Runge Kutta). Se non vuoi scriverti il codice del metodo basta prendere MatLab e usare odexyz e trovarti la soluzione.
Quello che segue era la mia prima risposta a questo problema. Non mi ero accorto che fosse una ODE e mi sono concentrato solo sul risponderti sulla condizione di Robin del BVP.
[size=150]EDIT[/size]
Tutto quello scritto sotto vale per la condizione di Robin nel caso di problemi ai limiti, ma non risolve l'ultimo problema dell'OP.
Ora per me sono giorni un po' pesanti (come lo sono per te immagino), ma appena ho un po' di tempo do' meglio una occhiata al tuo codice.
Per la condizione di Robin a $x=0$ provo a dirti come farei io, senza passare dallo sviluppo di Taylor, per comodità:
Questa è $\frac{d}{dx}u(0) + \frac{q_0}{k_0 + k_1 u+k_2 u^2}=0$
Discretizzo questa riga
$\frac{u_0 -2u_1 +u_2}{h^2} + \frac{g}{k} =0$, ($\square$) dove con $u_0$ indico il valore della soluzione nel nodo fantasma, mentre $h$ è il passo spaziale. Chiaramente il problema ora è la ricerca del valore della soluzione al nodo fantasma, e questo viene ricavato dalla condizione di Robin stessa, che ora scrivo per esteso:
$\frac{u_2 - u_0}{2h} + \frac{q_0}{k_0 + k_1 u_1 + k_2 * u_{1}^{2}}=0$
Ora da qui ricavi $u_0$, lo sostituisci in $\square$ e ricavi l'espressione della prima riga, che poi modificherai di conseguenza
[size=150]Fine EDIT[/size]
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo