Condition aux limites périodique pour l'équation de la chaleur en] 0,1 [

13

Considérons une condition initiale lisse et l'équation de chaleur dans une dimension:

tu=xxu
dans l'intervalle ouvert ]0,1[ , et supposons que nous voulons le résoudre numériquement avec des différences finies.

Je sais que pour que mon problème soit bien posé, je dois le doter de conditions aux limites à x=0 et x=1 . Je sais que Dirichlet ou Neumann fonctionnent bien.

Si j'ai dans le premier cas N points intérieurs xk=kN+1 pourk=1,,N, alors j'aiNinconnues:uk=u(xk)pourk=1,,N, caruest prescrit aux limites.

Dans le deuxième cas, j'ai vraiment inconnues u 0 , , u N + 1 , et je sais comment utiliser le (homogène) Neumann BC pour discrétiser le laplacien à la frontière, par exemple avec l'adjonction de deux points fictifs x - 1 et x N + 2 et les égalités:N+2u0,,uN+1x1xN+2

u1u12h=0=uN+2uN2h

Ma question concerne la Colombie-Britannique périodique. J'ai le sentiment que je pourrais utiliser une équation, à savoir mais peut-être deux, puis j'utiliserais x u ( 0 ) = x u ( 1 )

u(0)=u(1)
xu(0)=xu(1)

mais je ne suis pas sûr. Je ne sais pas non plus combien d'inconnues je devrais avoir. Est-ce ?N+1

bela83
la source
Avez-vous des conditions aux limites de Dirichlet ou Neumann? Le nombre de cellules fantômes dépend de l'ordre d'approximation des conditions aux limites de Neumann.
ilciavo
@ilciavo, la question porte sur les conditions aux limites périodiques.
Bill Barth

Réponses:

8

La meilleure façon de le faire est (comme vous l'avez dit) d'utiliser simplement la définition des conditions aux limites périodiques et de configurer correctement vos équations dès le départ en utilisant le fait que . En fait, encore plus fortement, les conditions aux limites périodiques identifient x = 0 avec x = 1 . Pour cette raison, vous ne devez avoir qu'un seul de ces points dans votre domaine de solution. Un intervalle ouvert n'a pas de sens lors de l'utilisation de conditions aux limites périodiques car il n'y a pas de limite .u(0)=u(1)x=0x=1

Ce fait signifie que vous ne devez pas placer un point à car il est identique à x = 0 . Discrétisant avec N + 1 points, vous utilisez ensuite le fait que par définition, le point à gauche de x 0 est x N et le point à droite de x N est x 0 .x=1x=0N+1x0 xNxN x0

schéma d'une grille périodique

Votre PDE peut ensuite être discrétisé dans l'espace comme

t[x0x1xN]=1Δx2[xN2x0+x1x02x1+x2xN12xN+x0]

Cela peut être écrit sous forme de matrice comme A=[ - 2 1 0 0 1 1 - 2 1 0 0

tx=1Δx2Ax
A=[21001121000012110012].

Bien sûr, il n'est pas nécessaire de créer ou de stocker cette matrice. Les différences finies doivent être calculées à la volée, en prenant soin de traiter les premier et dernier points spécialement si nécessaire.

tu=xxu+b(t,x)
x[1,1)uRef(t,x)=exp(t)cos(5πx)b(t,x)=(25π21)exp(t)cos(5πx)
clear

% Solve: u_t = u_xx + b
% with periodic boundary conditions

% analytical solution:
uRef = @(t,x) exp(-t)*cos(5*pi*x);
b = @(t,x) (25*pi^2-1)*exp(-t)*cos(5*pi*x);

% grid
N = 30;
x(:,1) = linspace(-1,1,N+1);

% leave off 1 point so initial condition is periodic
% (doesn't have a duplicate point)
x(end) = [];
uWithMatrix = uRef(0,x);
uNoMatrix = uRef(0,x);

dx = diff(x(1:2));
dt = dx.^2/2;

%Iteration matrix:
e = ones(N,1);
A = spdiags([e -2*e e], -1:1, N, N);
A(N,1) = 1;
A(1,N) = 1;
A = A/dx^2;

%indices (left, center, right) for second order centered difference
iLeft = [numel(x), 1:numel(x)-1]';
iCenter = (1:numel(x))';
iRight = [2:numel(x), 1]';

%plot
figure(1)
clf
hold on
h0=plot(x,uRef(0,x),'k--','linewidth',2);
h1=plot(x,uWithMatrix);
h2=plot(x,uNoMatrix,'o');
ylim([-1.2, 1.2])
legend('Analytical solution','Matrix solution','Matrix-free solution')
ht = title(sprintf('Time t = %0.2f',0));
xlabel('x')
ylabel('u')
drawnow

for t = 0:dt:1
    uWithMatrix = uWithMatrix + dt*( A*uWithMatrix + b(t,x) );
    uNoMatrix = uNoMatrix + dt*(  ( uNoMatrix(iLeft) ...
                                - 2*uNoMatrix(iCenter) ...
                                  + uNoMatrix(iRight) )/dx^2 ...
                                + b(t,x) );
    set(h0,'ydata',uRef(t,x))
    set(h1,'ydata',uWithMatrix)
    set(h2,'ydata',uNoMatrix)
    set(ht,'String',sprintf('Time t = %0.2f',t))
    drawnow
end

Tracé de l'état initial

Tracé de solution à t = 0,5

Tracé de solution à t = 1,0

Tracé de solution à t = 2,0

Doug Lipinski
la source
1
Grande et simple solution !! au cas où quelqu'un en aurait besoin ici les implémentations en Python
ilciavo
X
@ bela83 Vous avez raison, il n'est pas nécessaire de spécifier autre chose que la condition initiale. Cela entraînerait un système surdéterminé. Tout ce que vous devez faire est d'être un peu prudent près des points finaux de l'intervalle pour être sûr d'enrouler correctement les choses régulièrement. Il existe de nombreuses façons valables de procéder.
Doug Lipinski
-1

Selon cela, vous devez imposer des conditions aux limites périodiques comme:

u(0,t)=u(1,t)uX(0,t)=uX(1,t)

Une façon de discrétiser l'équation de chaleur en utilisant implicitement Euler vers l'arrière est

un+1-unΔt=uje+1n+1-2ujen+1+uje+1n+1ΔX2

Résolution du système d'équations

[je-ΔtΔX2UNE][u1n+1u1n+1uNn+1]=[u1nu2nuNn]

UNE=[-2100001-2100001-2100001-2100001-2000001-2]

u0uN+1

u1-uN=0u2-u02ΔX-uN+1-uN-12ΔX=0

Selon la section 2.11 LeVeque, cela vous donne une précision de second ordre pouruX

Enfin, votre système d'équations ressemblera à:

[01000-10-101010-100000je-ΔtΔX2UNE0000000][u0n+1u1n+1u2n+1uNn+1uN+1n+1]=[00u1nu2nuNn]

Ce qui vous donne N + 2 équations et N + 2 inconnues.

Vous pouvez également vous débarrasser du premier des équations et des cellules fantômes et arriver à un système de N équations et N inconnues.

ilciavo
la source
I don't understand the statement : "add two additional cells u1 and uN" because uN is already a point in ]0,1[. I have in mind xk=kN+1 so that xN=NN+1.
bela83
It's just a matter of indexing. You start with N cells (or points) from u0 to uN1 and you add two cells u1 and uN. If you have cells going from u1 to uN then you need to add cells at u0 and uN+1
ilciavo
OK then I don't understand the two more equations ! The first should be (with the notation from the question): u0=uN+1 right ? But I don't understand the second one: why not choose a centered difference approximation ? Finally, that makes N+1 unknowns, if first and last values are equal. Please compare to the situation with (homogeneous) Neumann BC in the question.
bela83
I've change the notation. It depends on the order of approximation for ux. The first comes from u(0,t)=u(1,t) and the second from ux(0,t)=ux(1,t)
ilciavo
1
u1=uN adds an additional restriction(equation u1−uN=0) to your system
ilciavo