Condizioni di contorno periodiche

[5mrkv]

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\))

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.

[5mrkv]

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.

[5mrkv]

Bamp.