Sistema equazioni differenziali in matlab

[luca.milano]

Ciao a tutti, sto cercando di risolvere questo sistema di equazioni differenziali, ma la soluzione che trovo non è per niente fisica, quindi molto probabilmente sto sbagliando a scrivere il codice. Questo è il sistema in questione:

$ { ( (dTf)/(dx)=h*(Ts-Tf)/(rhouCp) ),( (dq)/(dx)=-beta*qr ),( (d^2Ts)/(dx^2)=1/(ks)*((-beta_*qr)-h*(Ts-Tf)) ):} $

Le variabili dell'equazione sono Tf, Ts e qr. Le altre sono tutte costanti
Il codice che ho scritto è il seguente, per quanto riguarda il main (tralascio ovviamente la definizione di tutte le costanti):

save odeprova45
y0 = [Tf1,50000,q_sun,0];

tspan = 0:0.01:L;

[t,y] = ode45('fun',tspan,y0);

Per quanto riguarda invece la function:

function yp=fun(t,y)
load odeprova45

Tf=y(1);
Ts=y(2);
qr=y(3);
dTs=y(4);

yp=[h*(Ts-Tf)/(rho_u*cp);
dTs;
-beta_*qr
1/ks*(-beta_*qr)-h*(Ts-Tf) %d2Ts/dx2=d2Ts
];

è il comando giusto? Anche usando ode15i mi vengono risultati un po' strani (anche se ode15i ha una sintassi completamente diversa).

Spero in un vostro aiuto, perchè con Matlab faccio davvero fatica

[feddy]

Benvenuto nel forum,

first of all: impara a scrivere le formule/codici

non vedo il motivo per cui devi fare una function per risolvere l'equazione. Prova a seguire in ordine questo che ti scrivo, e che non dovresti avere problemi a riprodurre. L'unica cosa che cambia veramente è il modo (meno laborioso) di scrivere la function.

- Definisci le costanti

- Scriviti il sistema come un sistema del primo ordine (cosa che mi pare tu abbia fatto, visto che il termine noto ha 4 componenti).
Per fare questo basta definire una funzione. Ad esempio

Codice:

fun=@(t,y) [y(2)^2; sin(y(1))]

definisce il termine destro $f$ del corrispondente sistema differenziale $\mathbf{y}'=f(t,\mathbf{y})$.

- Scriviti il termine noto $y_0$, cosa che hai già fatto.

- Nel vettore con i tempi non serve che definisca tu il passo, è sufficiente scrivere

Codice:

tspan=[0 L]

, poi il metodo numerico che andrai ad applicare utilizzerà un passo variabile. A meno che tu non voglia esplicitamente che segua il passo da te fornito.

- Scrivi infine

Codice:

[t,y]=ode45(fun,tspan,y_0)

Nel caso non riuscissi posta il codice, con costanti ecc.

Anche usando ode15i mi vengono risultati un po' strani (anche se ode15i ha una sintassi completamente diversa).

E' un comando abbastanza particolare, questo solver è stato creato per sistemi differenziali di un certo tipo. Non credo abbia molto senso usarlo in un contesto come questo. Se frequenterai un corso di metodi numerici per equazioni differenziali sarai sicuramente in grado di scrivere un algoritmo in grado di risolvere il sistema, senza doverti affidare a odexyz.

[luca.milano]

Innanzitutto ti ringrazio dell'aiuto: sono riuscito a implementare ma la soluzione resta uguale alla precedente.
Per cui comincio a pensare che io debba porre delle condizioni al contorno diverse. Il comando ode permette di porre le condizioni al contorno solo nel punto iniziale del dominio, io invece vorrei mettere delle condizioni al contorno in altri punti del dominio. In particolare vorrei porre queste condizioni:
$ { ( T(0)=Tf1 ),( qr(0)=qsun ),( (dTs)/dx (0)=0 ),( (dTs)/dx(xrarr infty) =0):} $


PS
Per quanto riguarda l'equazione:
$ (dq)/(dx)=-beta*qr $

Che ha come condizione al contorno:
$ qr(0)=qsun $

Questa è facilmente integrabile con soluzione:
$ qr(x)=qsun*e^(-beta*x) $

C'è un modo per inserirla direttamente in forma algebrica o è necessario porla come equazione differenziale + condizione al contorno?

[feddy]

C'è un modo per inserirla direttamente in forma algebrica o è necessario porla come equazione differenziale + condizione al contorno?

Certamente, la definisci come una vera e propria funzione $qr(x)$, e verrà definita tramite

Codice:

qr =@(x) qsun*exp(-beta*x);

Per cui comincio a pensare che io debba porre delle condizioni al contorno diverse.

Il fatto che tu voglia mettere la condizione al bordo $x=\+infty$ cambia le cose, perché non si tratta più di un problema ai valori iniziali, ma è più un problema ai limiti e vanno usate altre tecniche. Forse dovresti fornire il testo completo dell'esercizio, il contesto, o comunque da dove l'hai preso. Perché non sapendo di cosa si parla non è semplice per me aiutarti

[luca.milano]

Grazie per la pazienza! Io provo a porti tutto il problema, considera che prima ti stavo proponendo una versione semplificata che stavo iniziando a fare per cominciare a prendere confidenza con Matlab. Si tratta di uno scambiatore di calore con un termine sorgente dato dall'irraggiamento solare. Le equazioni sono giuste poichè sono prese da un articolo scientifico. In ogni caso, il sistema di equazioni in questione è questo:

$ { (d/dx(rho*c_p*u*Tf)=d/dx(kf*d/dx(Tf))+h*(Ts-Tf)) ,(0=d/dx(ks*d/dx(Ts))+h*(Tf-Ts)-d/dx(q_r)):} $

dove le incognite sono:

$ Tf $ e $ Ts $ mentre $x$ va da $0$ ad $L$ (quest'ultimo è una costante).

E il termine $ q_(i n) $ vale:

$ q_r=q_(i n)*e^(-beta*x) $, in cui $q_(i n)$ dipende in realtà da $Ts$ come definirò sotto.

Tutti i termini presenti non sono davvero delle costanti come ti avevo detto nel messaggio iniziale, proprio perchè stavo iniziando a risolvere il problema in modo semplificato, con l'intenzione di complicarlo poco alla volta.

$ rho=(p)/(R*Tf) $ in cui $p$ e $R$ sono delle opportune costanti,

$ c_p=1060+0.449Tf+1.14*10^-3*Tf^2 $

$ u=m/(rho*A) $ in cui $m$ e $A$ sono anch'esse delle costanti,

$kf$ è un'altra polinomiale in funzione di Tf;

$h$ e $ks$ sono delle costanti

Ora le condizioni al contorno:

$ { ( Tf(0)=400 ),( (dTf)/(dx)(xrarr infty)=0 ),( (dTs)/(dx)(0)=0 ),((dTs)/(dx)(xrarr infty)=0),(q_(i n)=q_(sun)-alpha*(Ts(0)^4-Tamb^4)):} $


Il problema è questo, come mi consigli di procedere?

[feddy]

Mmm così su due piedi due cose:

1. Se hai risolto un problema diverso, mi pare giusto che la tua soluzione non sia "fisica"
2. Guardando le condizioni iniziali effettivamente non è un problema ai valori iniziali, ma un problema ai limiti. Potresti usare uno schema alle differenze finite, cioè discretizzare $[0,L]$ e discretizzare ogni derivata tramite con rapporti incrementali. Se non le hai mai viste o sentite però rischia di essere un bel problema, a meno che non esistano risolutori già pronti (che io di sicuro non conosco).

Tocca aspettare l'intervento di qualcuno più esperto, oppure chiedere delucidazioni al tuo prof

[luca.milano]

Cercando "problemi ai limiti" ho trovato il comando bvp per Matlab:

https://it.mathworks.com/help/matlab/ref/bvp4c.html

sai dirmi qualcosa a riguardo o non c'entra niente col mio problema? è un comando affidabile?

[feddy]

Certo, bvp è l'acronimo di boundary value problem, aka "problema ai limiti".
Nella descrizione dell'algoritmo (ultime righe) c'è esplicitamente scritto che usa uno schema alle differenze finite... quindi dovrebbe funzionare. Non ho mai usato questo comando, quindi non so come funzioni di preciso e non posso aiutarti, ma ci sono sicuramente tutorial ed esempi. Devi capire come passargli le condizioni ai bordi.

[luca.milano]

Per quanto riguarda l'ultima condizione al contorno:

$q_(i n)=q_(sun)-alpha*(Ts(0)^4-T_(amb)^4)$

Ha la particolarità di essere una condizione che dipende dalla soluzione stessa. Cercando un po' ho trovato che questo tipo di condizione è detta "Condizione di Robin" da contrapporsi alle condizioni di Dirichlet e Neumann.
Se fosse un problema ai valori iniziali, sapresti dirmi come implementarla col comando ode45? Può essere che una volta capito come funziona il comanda bvp potrei riuscire a trasferire la sintassi!
Questo perchè se non ho letto male sul mathworks, il comando bvp dovrebbe avere la stessa sintassi del comando ode.

[feddy]

Allora, condizioni di Dirichlet sono condizioni come la tua $T_f(0)=400$. Condizioni di Neumann sono quelle sulla/e derivata/e. Di Robin sono quelle di tipo misto.
Ma ti ho già detto come si lancia ode45, e hai detto di essere riuscito a lanciarlo. Non credo abbia senso quello che stia cercando di fare, insomma stai cercando di sviare il problema, e poi "trasferire la sintassi" non credo funzioni visto che ode45 non prende dati iniziali ai "bordi".
Hai provato a guardare gli esempi sulla documentazione del comando bvp?