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
Ma figurati, sei stato gentilissimo, prenditi pure il tempo che ti serve per fare i tuoi esami e più avanti se avrai voglia potrai ripassare per questa discussione. In ogni caso mi hai già aiutato tanto e per questo ti ringrazio!
Purtroppo il tuo ultimo commento ammetto di non averlo capito benissimo. In particolare, lo stesso problema ma svolto in modo lineare, quindi con $ k $ costante nota, riesco a farlo funzionare e forse per questo pensavo fosse relativamente semplice estenderlo al caso non lineare utilizzando il metodo di Newton.
Giusto per completezza ti riporto il codice del caso lineare, che magari guarderai più avanti quando avrai più tempo:
Lo avevo scritto prima che tu mi spiegassi come scrivere le matrici agilmente, quindi ci sono ancora i cicli for.
Buona fortuna per i tuoi esami!
Purtroppo il tuo ultimo commento ammetto di non averlo capito benissimo. In particolare, lo stesso problema ma svolto in modo lineare, quindi con $ k $ costante nota, riesco a farlo funzionare e forse per questo pensavo fosse relativamente semplice estenderlo al caso non lineare utilizzando il metodo di Newton.
Giusto per completezza ti riporto il codice del caso lineare, che magari guarderai più avanti quando avrai più tempo:
clear clc %CODICE CAMPO DI T CON ASSEGNATA TEMPERATURA E FLUSSO IN INGRESSO L=10; N=100; %numero elementi dx=L/N; %SPECIFICA LA TEMPERATURA ALL'INGRESSO E IL FLUSSO TERMICO T0=100; q0=-500; %COEFFICIENTI DELL'EQUAZIONE: generazione di potenza e conduttività termica g=250; k=50; %conduttività termica A=zeros(N,N); b=zeros(N,1); A(1,1)=1; %prima riga A(2,2)=+1; %seconda riga b=[T0 (T0-dx/k*q0-1/2*dx^2*g/k)]'; for j=3:N A(j,j-2)=1; A(j,j-1)=-2; A(j,j)=1; b(j,1)=-g/k*dx^2; end T=A\b; x=linspace(0,L,N); plot(x,T)
Lo avevo scritto prima che tu mi spiegassi come scrivere le matrici agilmente, quindi ci sono ancora i cicli for.
Buona fortuna per i tuoi esami!
Buona fortuna per i tuoi esami!
Grazie !
Ad ogni modo ti ripeto che per risolvere ODE, che è quello che vuoi fare, devi ricondurre l'equazione diff. (che era al secondo ordine) in un sistema al prim'ordine e usare MatLab. Il fatto che plottasse qualcosa non vuol dire niente in realtà. Ad esempio, tu hai messo esplicitamente una condizione sul bordo a $x=L$. Quindi, faccio finta che sulle ascisse ci sia il tempo e ti scrivo 5 righe di come andrebbe fatto uno stupidissimo Eulero Esplicito per esempio col tuo problema
\( \begin{cases} (d^2 u)/dx^2 +g/k=0 \\ u(0)=u_0\\ \frac{du}{dx}(0)=-\frac{q_0}{(k(u_0))} \end{cases} \)
Ho scritto il problema al prim'ordine introducendo $u_{1}'=u_2$ e risolto un sistema differenziale con condizioni iniziali che ho scelto io, visto che tu precedentemente non le davi (ho solo assunto $u(0)=u_0=0$ poiché suppongo che la temperatura all'istante iniziale sia nulla, comunque non cambia niente, basta che sistemi $y_0$). Ancora, visto che il risultato comparirà sulla tua tesi, se vuoi garantirti un risultato migliore, senza preoccuparti della tanta teoria che ci sta dietro, puoi usare il comando ode di MatLab, che mi pare tu conosca già.
Ad ogni modo, i tuoi problemi sono più che altro sulla teoria, e continuare a scrivere pezzi di codice di problemi diversi rende questa discussione veramente "enorme" e si rischia di perdere di vista gli obiettivi. Ti ripeto che secondo me ti conviene darti una letta a del materiale sull'argomento, almeno saper distinguere una ODE da un problema ai limiti, altrimenti diventa impossibile andare avanti
clear all close all tf=5; #tempo finale ts=1000; #time steps k=tf/ts; #passo temporale #costanti T0=5; q0=0; g=0; k0=100; %primo termine della polinomiale per la conduttività termica: k=k0+aT k1=20; %coefficiente della T f=@(y) [y(2);-g./(k0 + k1*y(1))]; df=@(y) [1,0;0,(g*k1)/((k0+k1*y(1)).^2)]; #Eulero esplicito y0=[0;2]; #DATO INIZIALE: corrisponde a u(0) e u'(0)=-g/k. y=NaN(2,ts+1); y(:,1)=y0; for n=1:ts y(:,n+1)=y(:,n)+k*f(y(:,n)); end t=linspace(0,tf,ts+1); plot(t,y(1,:))
Ah quindi mi stai dicendo che anche la soluzione che ho trovato nel caso lineare fosse sbagliata? Eh si che sembrava tanto carina
Non so cosa debba venire, ma se non hai condizioni a un bordo non te le puoi certo inventare.
Hai ragione, probabilmente ho bisogno di capire prima un po' meglio la teoria e cercherò di leggere bene il testo che mi hai passato. A proposito, se hai altri testi, in particolare con esempi di codici matlab, sono ben accetti.
Non capisco proprio a cosa ti riferisci quando dici che io mi sia inventato una condizione a un bordo: io ho impostato la funzione al nodo $1$ e la funzione al nodo $2$ (nota la derivata prima in $0$) e la soluzione trovata nel caso lineare viene corretta (l'ho confrontata con la soluzione analitica e con altri risolutori, quindi sono sicuro che sia corretta).
Non so cosa debba venire, ma se non hai condizioni a un bordo non te le puoi certo inventare
Non capisco proprio a cosa ti riferisci quando dici che io mi sia inventato una condizione a un bordo: io ho impostato la funzione al nodo $1$ e la funzione al nodo $2$ (nota la derivata prima in $0$) e la soluzione trovata nel caso lineare viene corretta (l'ho confrontata con la soluzione analitica e con altri risolutori, quindi sono sicuro che sia corretta).
Un testo recente è "Modellistica numerica per problemi differenziali - Quarteroni".
Per quanto riguarda l'ultima cosa: quello che intendo è che hai risolto un problema ai valori iniziali con un metodo per problemi ai limiti, visto che non conoscevi la condizione al bordo destro. Questo è sbagliato, che venga "giusta" poi è un altro conto. Oltretutto, quello che hai fatto te è stato mettere la condizione di Dirichlet in $0$ e nel nodo successivo quella di Neumann... ma converrai che così non hai risolto veramente il problema assegnato, visto che non hai modificato lo stesso nodo. Insomma, a scanso di equivoci: per ODE si usano metodi per ODE, per BVP differenze finite
Per quanto riguarda l'ultima cosa: quello che intendo è che hai risolto un problema ai valori iniziali con un metodo per problemi ai limiti, visto che non conoscevi la condizione al bordo destro. Questo è sbagliato, che venga "giusta" poi è un altro conto. Oltretutto, quello che hai fatto te è stato mettere la condizione di Dirichlet in $0$ e nel nodo successivo quella di Neumann... ma converrai che così non hai risolto veramente il problema assegnato, visto che non hai modificato lo stesso nodo. Insomma, a scanso di equivoci: per ODE si usano metodi per ODE, per BVP differenze finite
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo