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

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

$$\partial_t u = \partial_{xx} u$$ 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 $x_k=\frac{k}{N+1}$ pour$k=1,\cdots,N$, alors j'ai$N$inconnues:$u_k=u(x_k)$pour$k=1,\cdots,N$, car$u$est 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+2$$u_0,\cdots,u_{N+1}$$x_{-1}$$x_{N+2}$

$$\frac{u_1-u_{-1}}{2 h} = 0 = \frac{u_{N+2}-u_N}{2 h}$$

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)$$ $$\partial_x u(0) = \partial_x u(1)$$

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

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=0$$x=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=1$$x=0$$N+1$$x_0$ $x_N$$x_N$ $x_0$

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

$$\frac{\partial}{\partial t} \left[\begin{array}{c} x_0 \\ x_1 \\ \vdots \\ x_N \end{array}\right] = \frac{1}{\Delta x^2} \left[\begin{array}{c} x_N - 2x_0 + x_1 \\ x_0 - 2x_1 + x_2 \\ \vdots \\ x_{N-1} - 2x_N + x_0 \end{array}\right]$$

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

$$\frac{\partial}{\partial t}\vec{x} = \frac{1}{\Delta x^2} \mathbf{A} \vec{x}$$ $$\mathbf{A} = \left[\begin{array}{c} -2 & 1 & 0 & \cdots & 0 & 1 \\ 1 & -2 & 1 & 0 & \cdots & 0 \\ & \ddots & \ddots & \ddots \\ && \ddots & \ddots & \ddots \\ 0 & \cdots & 0 & 1 & -2 & 1 \\ 1 & 0 & \cdots & 0 & 1 & -2 \end{array}\right].$$

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.

$$\partial_t u = \partial_{xx}u + b(t,x)$$ $x\in[-1,1)$$u_\text{Ref}(t,x) = \exp(-t)\cos(5\pi x)$$b(t,x) = (25\pi^2-1)\exp(-t)\cos(5\pi 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

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

$$\begin{equation} u(0, t) = u(1, t) \\ u_x(0, t) = u_x(1, t) \end{equation}$$

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

$$\begin{equation} \frac{u^{n+1}-u_{n}}{\Delta t}= \frac{u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i+1}}{\Delta x^2} \end{equation}$$

Résolution du système d'équations

$$\begin{equation} \begin{bmatrix} I-\frac{\Delta t}{\Delta x^2}A \end{bmatrix} \begin{bmatrix} u^{n+1}_1 \\ u^{n+1}_1 \\ \vdots \\ \\ \\ u^{n+1}_{N} \end{bmatrix} = \begin{bmatrix} u^n_1 \\ u^n_2 \\ \vdots \\ \\ \\ u^n_{N} \end{bmatrix} \end{equation}$$

$$\begin{equation} A= \begin{bmatrix} -2 & 1 & 0 & 0 & 0 & \ldots & 0 \\ 1 & -2 & 1 & 0 & 0 &\ldots & 0 \\ 0 & 1 & -2& 1 & 0 & \ldots & 0 \\ 0 & 0 & 1& -2 & 1 & \ldots & 0 \\ 0 & 0 & 0& 1 & -2 & \ldots & 0 \\ \vdots &\vdots &\vdots &\vdots & \ddots &\ddots &\vdots \\ 0 & 0 & 0 & \ldots &0 &1 &-2 \end{bmatrix} \end{equation}$$

$u_{0}$$u_{N+1}$

$$\begin{equation} u_1 - u_{N} = 0 \\ \frac{u_{2}-u_{0}}{2 \Delta x} - \frac{u_{N+1}-u_{N-1}}{2 \Delta x} = 0 \end{equation}$$

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

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

$$\begin{equation} \begin{bmatrix} 0 & 1 & 0 & 0 & \ldots & 0 & -1 & 0 \\ -1 & 0 & 1 & 0 & \ldots & 1 & 0 & -1 \\ 0 & & & & & & & 0\\ 0 & & & & & & & 0\\ 0 & & & I-\frac{\Delta t}{\Delta x^2}A & & & & 0\\ 0 & & & & & & & 0\\ 0 & & & & & & & 0\\ 0 & & & & & & & 0 \end{bmatrix} \begin{bmatrix} u^{n+1}_{0}\\ u^{n+1}_1 \\ u^{n+1}_2 \\ \\ \\ \\ \\ u^{n+1}_{N} \\ u^{n+1}_{N+1} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ u^n_1 \\ u^n_2 \\ \\ \\ \\ \\ u^n_{N} \end{bmatrix} \end{equation}$$

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.