Diffusione attraverso due layer

[vinxp1]

Salve, sono nuovo e mi sono appena presentato. Chiedo cortesemente il vostro aiuto per la risoluzione della seguente equazione differenziale alle derivate parziali (corredata delle condizioni a contorno):

$(delc)/(delt)$=D/r^2*$(del)/(delr)$(r^2*$(delc)/(delr)$)

-Il primo layer va da 0 a un raggio Ri, in cui il coefficiente di diffusione è D=Df e le condizioni a contorno sono le seguenti:
$\{(t=0, c=cf),(r=0, (∂cf)/(∂r)=0),(r=Ri, (∂cf)/(∂r)=(∂cs)/(∂r))}$
L'ultima condizione è quella all'interfaccia che può anche essere scritta come: r=Ri, cf=10*cs

-Il secondo layer va dal raggio Ri al raggio Rs, in cui il coefficiente di diffusione è D=Ds e le condizioni a contorno sono le seguenti:
$\{(t=0, c=0),(r=Ri, (∂cf)/(∂r)=(∂cs)/(∂r)),(r=Rs, (∂cs)/(∂r)=0)}$

Dispongo di Matlab, e vorrei realizzare un qualche script con il metodo delle differenze finite ma non riesco proprio (altri metodi sono ben accetti ). Grazie anticipatamente.

[Raptorista1]

Ciao, hai mai scritto uno script a differenze finite in matlab? Puoi partire da una cosa semplice, per esempio risolvere l'equazione con uno strato solo, per poi aggiungere il secondo strato. Sapresti fare almeno questo?

[vinxp1]

Ciao Raptorista, grazie per la risposta. Allego lo script che ho realizzato, sperando che sia scritto giusto. Ho però due problemi:
-il primo è come si inserisce il tempo nello script;
-il secondo è cosa scrivere nel command per far funzionare lo script(avevo pensato ad fsolve... )

Lo script è il seguente:

function f=diff(c)

Df=10^-5;
Ds=4*10^-3;
dr=1;
dt=1;
cf=10;

%CONDIZIONE INIZIALE
for i=1:Ri

f(i,1)=c(i,1)-cf;

end


%BOUNDARY E CONDIZIONI AI NODI
for j=1:10


f(1,j)=(c(2,j)-c(1,j))/dr;

for i=2:Ri-1


f(i,j)=(c(i,j-1)-c(i,j))/dt+(Df/(i*dr)^2)*(2*i*dr*(c(i+1,j)-c(i,j))/dr+(i*dr^2)*(c(i+1,j)-2*c(i,j)+(c(i-1,j)))/dr^2);

end

f(Ri,j)=c(Ri,j)-10*c(Ri+1,j);
end

[Raptorista1]

"vinxp":Allego lo script che ho realizzato, sperando che sia scritto giusto.

Questo è uno script per cosa? Il problema completo con due strati?
Potresti metterlo tra i tag di codice forniti dal forum? E indentarlo come si deve in modo che sia leggibile.

"vinxp":
-il primo è come si inserisce il tempo nello script;

Puoi tenerne traccia in una variabile, o contare il numero di passi da eseguire se il passo temporale è fissato.

"vinxp":
-il secondo è cosa scrivere nel command per far funzionare lo script(avevo pensato ad fsolve... )

Non capisco XD

[vinxp1]

Scusami ma per me è tutto un dubbio...facciamo così, questo è lo script completo per due strati dove i è lo spazio e j è il tempo (anche se non so se è giusto così)...e provo a riscriverlo in modo più ordinato:

function f=diffusion(c)

$Df=10^-5$;
$Ds=4*10^-3$;
$dr=1$;
$dt=1$;
$c_f=10$;

%CONDIZIONE INIZIALE
for i=1:Ri

$f_(i,1)=c_(i,1)-c_f$;

end

%BOUNDARY E CONDIZIONI AI NODI
for j=1:10

$f_(1,j)=(c_(2,j)-c_(1,j))/(dr)$;

for i=2:(Ri-1)

$f_(i,j)=(c_(i,j-1)-c_(i,j))/(dt)+((Df)/(i*dr^2))*(2*i*dr*(c_(i+1,j)-c_(i,j))/(dr)+(i*dr^2)*(c_(i+1,j)-2*c_(i,j)+c_(i-1,j))/(dr^2))$; questa è l'equazione differenziale discretizzata

end

$f_(Ri,j)=c_(Ri,j)-10*c_(Ri+1,j)$;
end

Cosa ne pensi?

[Raptorista1]

Apprezzo il tentativo, ma quello che dovevi fare era di scrivere il codice così:

function f=diff(c)

Df=10^-5;
Ds=4*10^-3;
dr=1;
dt=1;
cf=10;

%CONDIZIONE INIZIALE
for i=1:Ri
    f(i,1)=c(i,1)-cf;
end


%BOUNDARY E CONDIZIONI AI NODI
for j=1:10
    f(1,j)=(c(2,j)-c(1,j))/dr;
   
    for i=2:Ri-1
        f(i,j)=(c(i,j-1)-c(i,j))/dt+(Df/(i*dr)^2)*(2*i*dr*(c(i+1,j)-c(i,j))/dr+(i*dr^2)*(c(i+1,j)-2*c(i,j)+(c(i-1,j)))/dr^2);
    end

    f(Ri,j)=c(Ri,j)-10*c(Ri+1,j);
end

Mi confermi che questo è il codice che hai per le mani?

Poi, altre cose non mi sono chiare.
\(cf\) è una cosa unica o è \(c\cdot f\)? E perché nelle condizioni al bordo c'è una derivata rispetto a \(cf\)?
Stesse domande per \(cs\).
Che cos'è [inline]f[/inline] nel codice?
Servono più commenti, che cosa fa il ciclo [inline]for[/inline] innestato con la formula lunga?
Sai qualcosa della buona posizione del problema?