Periodische Randbedingung für die Wärmegleichung in] 0,1 [

Betrachten wir eine glatte Ausgangsbedingung und die Wärmegleichung in einer Dimension:

$$\partial_t u = \partial_{xx} u$$ im offenen Intervall $]0,1[$ , und nehmen wir an, dass wir sie numerisch mit endlichen Differenzen lösen wollen.

Ich weiß, dass ich mein Problem mit Randbedingungen bei $x=0$ und ausstatten muss, damit es gut gestellt ist $x=1$. Ich weiß, dass Dirichlet oder Neumann gut funktionieren.

Wenn ich im ersten Fall $N$ innere Punkte habe, ist fürk=1,⋯,N, dann habe ichNUnbekannte:uk=u(xk)fürk=1,⋯,N, weil$x_k=\frac{k}{N+1}$$k=1,\cdots,N$$N$$u_k=u(x_k)$$k=1,\cdots,N$ an den Grenzen vorgeschrieben ist.$u$

Im zweiten Fall habe ich wirklich Unbekannte u 0 , ⋯ , u N + 1 , und ich weiß, wie man den (homogenen) Neumann BC verwendet, um den Laplace an der Grenze zu diskretisieren, zum Beispiel mit dem Zusatz von zwei fiktiven Punkten x - 1 und x N + 2 und die Gleichungen:$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}$$

Meine Frage betrifft das periodische BC. Ich habe das Gefühl, dass ich eine Gleichung verwenden könnte, nämlich aber vielleicht zwei, und dann würde ich ∂ x u ( 0 ) = ∂ x u

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

aber ich bin mir nicht sicher. Ich weiß auch nicht, wie viele Unbekannte ich haben sollte. Ist es $N+1$ ?

Der beste Weg, dies zu tun, ist (wie Sie sagten), einfach die Definition der periodischen Randbedingungen zu verwenden und Ihre Gleichungen von Anfang an korrekt unter Verwendung der Tatsache aufzustellen, dass . Tatsächlich identifizieren periodische Randbedingungen noch stärker x = 0 mit x = 1 . Aus diesem Grund sollten Sie nur einen dieser Punkte in Ihrer Lösungsdomäne haben. Ein offenes Intervall ist bei periodischen Randbedingungen nicht sinnvoll, da es keine Grenze gibt .$u(0)=u(1)$$x=0$$x=1$

Diese Tatsache bedeutet, dass Sie keinen Punkt bei da er mit x = 0 identisch ist . Diskretisierung mit N + 1 Punkten, verwenden Sie dann die Tatsache , dass der Punkt links von definitions x 0 heißt x N und der Punkt rechts von x N ist$x=1$$x=0$$N+1$$x_0$ $x_N$$x_N$ $x_0$ .

Ihre PDE kann dann im Raum als ∂ diskretisiert werden

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

Dies kann in Matrixform als wobei 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].$$

Natürlich muss diese Matrix nicht erstellt oder gespeichert werden. Die endlichen Differenzen sollten im laufenden Betrieb berechnet werden, wobei darauf zu achten ist, dass der erste und der letzte Punkt nach Bedarf behandelt werden.

Als einfaches Beispiel löst das folgende MATLAB-Skript mit periodischen Randbedingungen auf der Domäne x ∈ [ - 1 , 1 ) . Die hergestellte Lösung U Ref ( t , x ) = exp ( - t ) cos ( 5 π x ) = ( 25 π 2

$$\partial_t u = \partial_{xx}u + b(t,x)$$ $x\in[-1,1)$$u_\text{Ref}(t,x) = \exp(-t)\cos(5\pi x)$ verwendet wird, was bedeutet , . Ich habe der Einfachheit halber die Euler-Vorwärts-Zeitdiskretisierung verwendet und die Lösung sowohl mit als auch ohne Bildung der Matrix berechnet. Die Ergebnisse sind unten gezeigt.$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

Nach dieser sollten Sie periodische Randbedingungen wie zu verhängen:

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

Eine Möglichkeit, die Wärmegleichung implizit mit dem Rückwärts-Euler zu diskretisieren, ist

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

Lösen des Gleichungssystems

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

Wo

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

Sie können Ihre periodischen Randbedingungen einbeziehen, indem Sie zwei weitere Gleichungen und zwei zusätzliche (Geister-) Zellen hinzufügen $u_{0}$ und $u_{N+1}$ so dass:

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