Condición límite periódica para la ecuación de calor en] 0,1 [

Consideremos una condición inicial suave y la ecuación de calor en una dimensión:

$$\partial_t u = \partial_{xx} u$$ en el intervalo abierto $]0,1[$ , y supongamos que queremos resolverlo numéricamente con diferencias finitas.

Yo sé que para mi problema a ser bien planteado, que necesito para dotarlo de las condiciones de contorno en $x=0$ y $x=1$ . Sé que Dirichlet o Neumann funcionan bien.

Si tengo en el primer caso $N$ puntos interiores $x_k=\frac{k}{N+1}$ para$k=1,\cdots,N$, entonces tengo$N$incógnitas:$u_k=u(x_k)$para$k=1,\cdots,N$, porque$u$está prescrito en los límites.

En el segundo caso, realmente tengo incógnitas de u 0 , ⋯ , u N + 1 , y sé cómo usar el Neumann BC (homogéneo) para discretizar el laplaciano en el borde, por ejemplo con la unión de dos puntos ficticios x - 1 y x N + 2 y las igualdades:$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}$$

Mi pregunta es sobre BC periódica. Tengo la sensación de que podría usar una ecuación, es decir, pero tal vez dos, y luego usaría ∂ x u ( 0 ) = ∂ x u ( 1

$$u(0) = u(1)$$ $$\partial_x u(0) = \partial_x u(1)$$

pero no estoy seguro. Tampoco sé cuántas incógnitas debería tener. ¿Es ?$N+1$

La mejor manera de hacer esto es (como dijiste) simplemente usar la definición de condiciones de contorno periódicas y configurar tus ecuaciones correctamente desde el principio usando el hecho de que . De hecho, aún más fuertemente, las condiciones de contorno periódicas identifican x = 0 con x = 1 . Por esta razón, solo debe tener uno de estos puntos en su dominio de solución. Un intervalo abierto no tiene sentido cuando se usan condiciones de límite periódicas ya que no hay límite .$u(0)=u(1)$$x=0$$x=1$

Este hecho significa que no debe colocar un punto en ya que es lo mismo que x = 0 . Discretizando con N + 1 puntos, entonces usa el hecho de que, por definición, el punto a la izquierda de x 0 es x N y el punto a la derecha de x N es$x=1$$x=0$$N+1$$x_0$ $x_N$$x_N$ $x_0$ .

Su PDE se puede discretizar en el espacio como

$$\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]$$

Esto se puede escribir en forma de matriz como donde A=[ - 2 1 0 ⋯ 0 1 1 - 2

$$\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].$$

Por supuesto, no hay necesidad de crear o almacenar esta matriz. Las diferencias finitas deben calcularse sobre la marcha, teniendo cuidado de manejar el primer y el último punto especialmente según sea necesario.

Como ejemplo simple, el siguiente script de MATLAB resuelve con condiciones de contorno periódicas en el dominio x ∈ [ - 1 , 1 ) . Se utiliza la solución fabricada u Ref ( t , x ) = exp ( - t ) cos ( 5 π x ) , lo que significa b ( t

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

De acuerdo con esto , debe imponer condiciones de contorno periódicas como:

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

One way of discretising the Heat Equation implicitly using backward Euler is

$$\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}$$

Solving the system of equations

$$\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}$$

Where

$$\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}$$

Your periodic boundary conditions can be included by adding two more equations and two additional(ghost) cells $u_{0}$ and $u_{N+1}$ such that:

$$\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}$$

According to Section 2.11 LeVeque this gives you a 2nd order accuracy for $u_x$

Finally your system of equations will look like:

$$\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}$$

Which gives you N+2 equations and N+2 unknowns.

You can also get rid of the first to equations and the ghost cells and arrive at a system of N equations and N unknowns.