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.

Risposte
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.

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.