Condizioni di contorno periodiche
Data, con \(x \in [a,b], \ t \in [t_{0},T]\), l'equazione
\begin{cases}
u_{t}+au_{x}&=0 \\
u(x,0)&=f(x) \\
\end{cases}
ho le condizioni di contorno
\begin{split}
&u(a,t)&=u(b,t) \\
&u(a+h,t)&=u(b+h,t) \\
\end{split}
Il problema è che la soluzione che mi viene non assomiglia a quella esatta. Ho fatto così: dato che conosco il valore della soluzione al tempo \(t=0\) discretizzando il dominio spaziale ottengo un vettore soluzione al tempo \(t=0\). Quello che devo fare quindi è discretizzare il dominio temporale e calcolare il vettore soluzione al tempo \(t=1\), e così via fino al tempo finale. Quello che le condizioni di contorno mi dicono è che (ad esempio per vettore tempo iniziale \(j=1\))
Nel programma calcolo la soluzioni con il metodo di Euler in avanti nei punti del dominio discretizzato \(x=a+h,...,b\) e non in \(x=a\) in cui viene sovra scritta la condizione di contorno. In più viene aggiunto un punto in più a destra del vettore soluzione, sempre attraverso le condizioni di contorno, per quei metodi del tipo \(u(i,j+1)=f(u(i+1,j))\) che vogliono questo punto in più per calcolare la soluzione all'estrema destra del vettore soluzione. Metodo usato:
\begin{split}
u_{i,j+1}=u_{i,j}-akh^{-1}(u_{i+1,j}-u_{i,j})
\end{split}
In rosso è la soluzione esatta mentre il verde quella approssimata, per un confronto al tempo \(t=10\) (vedi commento blocco verifica)
script.m
eulerAv.m
f1.m (cond. iniziale)
f2.m (sol.)
All'interno del programma c'è una seconda opzione per una diversa condizione di contorno.
\begin{cases}
u_{t}+au_{x}&=0 \\
u(x,0)&=f(x) \\
\end{cases}
ho le condizioni di contorno
\begin{split}
&u(a,t)&=u(b,t) \\
&u(a+h,t)&=u(b+h,t) \\
\end{split}
Il problema è che la soluzione che mi viene non assomiglia a quella esatta. Ho fatto così: dato che conosco il valore della soluzione al tempo \(t=0\) discretizzando il dominio spaziale ottengo un vettore soluzione al tempo \(t=0\). Quello che devo fare quindi è discretizzare il dominio temporale e calcolare il vettore soluzione al tempo \(t=1\), e così via fino al tempo finale. Quello che le condizioni di contorno mi dicono è che (ad esempio per vettore tempo iniziale \(j=1\))
intx=[-2, 2 ]; intt=[ 0, 0.5]; h=0.05; k=0.05; x=intx(1):h:intx(2); t=intt(1):k:intt(2); n=length(x); m=length(t); %u(i,j) u(1,1)=u(n,1); u(n+1,1)=u(2,1);
Nel programma calcolo la soluzioni con il metodo di Euler in avanti nei punti del dominio discretizzato \(x=a+h,...,b\) e non in \(x=a\) in cui viene sovra scritta la condizione di contorno. In più viene aggiunto un punto in più a destra del vettore soluzione, sempre attraverso le condizioni di contorno, per quei metodi del tipo \(u(i,j+1)=f(u(i+1,j))\) che vogliono questo punto in più per calcolare la soluzione all'estrema destra del vettore soluzione. Metodo usato:
\begin{split}
u_{i,j+1}=u_{i,j}-akh^{-1}(u_{i+1,j}-u_{i,j})
\end{split}
In rosso è la soluzione esatta mentre il verde quella approssimata, per un confronto al tempo \(t=10\) (vedi commento blocco verifica)
script.m
%EQUAZIONE DI TRASPORTO intx=[-2, 2 ]; intt=[ 0, 0.5]; p=0; %p=1; eulerAv(@f1,@f2,intx,intt,p); %eulerIn(@f1,@f2,intx,intt,p); %laxFr (@f1,@f2,intx,intt,p);
eulerAv.m
%METODO DI EULERO ALL'AVANTI function eulerAv(f1,f2,intx,intt,p) h=0.05; k=0.05; x=intx(1):h:intx(2); t=intt(1):k:intt(2); n=length(x); m=length(t); %inizializzo il vettore for i=1:n, x1=x(i); u(i,1)=f1(x1); end if p u(1,1)=u(2,1); u(n+1,1)=u(n,1); elseif not(p) u(1,1)=u(n,1); u(n+1,1)=u(2,1); end %eulero all'avanti for j=1:m-1, for i=2:n, x1=intx(1)*(h^-1)*k; x2=(u(i,j)-u(i+1,j))*x1; u(i,j+1)=u(i,j)+x2; %#ok<*AGROW> end if p u(1,j)=u(2,j); u(n+1,j)=u(n,j); elseif not(p) u(1,j)=u(n,j); u(n+1,j)=u(2,j); end end %blocco verifica T=ones(1,n+1); T=T*t(m-1); x=[x,x(n)+h]; v=f2(x,T); %plot(x,u(:,1) ,'*k '); hold on plot(x,u(:,m-1),'g-.'); hold on plot(x,v ,'r--');
f1.m (cond. iniziale)
function a=f1(x) a=exp(-5*x.^2);
f2.m (sol.)
function a=f2(x,t) a=exp(-5*(x+2*t).^2);
All'interno del programma c'è una seconda opzione per una diversa condizione di contorno.
Risposte
Risolto. Ora ho un problema con il metodo di Euler all'indietro. Conosco il vettore soluzione al tempo iniziale \(u(:,1)\). Lo schema mi dice:
\begin{split}
& \alpha=akh^{-1} \\
& (1-\alpha)u_{i,j+1}+\alpha u_{i+1,n+1}=u_{i,n} \\
\end{split}
Ma non riesco ad impostare il sistema.
\begin{split}
& \alpha=akh^{-1} \\
& (1-\alpha)u_{i,j+1}+\alpha u_{i+1,n+1}=u_{i,n} \\
\end{split}
Ma non riesco ad impostare il sistema.
Bamp.