Периодическое граничное условие для уравнения теплопроводности в] 0,1 [

Давайте рассмотрим гладкое начальное условие и уравнение теплопроводности в одном измерении: в открытом интервале и предположим, что мы хотим решить его численно с конечными разностями.

\partial_{T} U знак равно \partial_{Икс Икс} U

$\partial_t u = \partial_{xx} u$

] 0, 1 [

$]0,1[$

Я знаю, что для того, чтобы моя задача была правильно поставлена, мне нужно наделить ее граничными условиями при и . Я знаю, что Дирихле или Нейман работают хорошо. $x=0$ $x=1$

Если у меня в первом случае внутренних точек для , то у меня есть неизвестных: для , потому что прописан на границах. $N$ $x_k=\frac{k}{N+1}$ $k=1,\cdots,N$ $N$ $u_k=u(x_k)$ $k=1,\cdots,N$ $u$

Во втором случае у меня действительно есть $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}

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

Мой вопрос о периодической до нашей эры. У меня есть ощущение, что я мог бы использовать одно уравнение, а именно но, может быть, два, и тогда я бы использовал

u (0) = u (1)

$u(0) = u(1)$

\partial_{x} u (0) = \partial_{x} u (1)

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

но я не уверен. Я тоже не знаю, сколько неизвестных мне нужно иметь. Это ? $N+1$

pde boundary-conditions parabolic-pde bela83
источник

У вас есть граничные условия Дирихле или Неймана? Количество клеток-призраков зависит от порядка аппроксимации для вас граничных условий Неймана.

ilciavo

@ilciavo, вопрос о периодических граничных условиях.

Билл Барт

Ответы:

Лучший способ сделать это - (как вы сказали) просто использовать определение периодических граничных условий и правильно настроить ваши уравнения с самого начала, используя тот факт, что . Фактически, еще сильнее периодические граничные условия отождествляют с . По этой причине у вас должна быть только одна из этих точек в вашем домене решений. Открытый интервал не имеет смысла при использовании периодических граничных условий, так как нет границы . $u(0)=u(1)$ $x=0$ $x=1$

Этот факт означает, что вам не следует помещать точку в поскольку она совпадает с . Дискретизируя с точками, вы затем используете тот факт, что по определению точка слева от равна а точка справа от равна . $x=1$ $x=0$ $N+1$ $x_0$ $x_N$ $x_N$ $x_0$

схема периодической сетки

Ваш PDE может быть затем дискретизирован в пространстве как

\frac{\partial}{\partial t} [\begin{matrix} x_{0} \\ x_{1} \\ ⋮ \\ x_{N} \end{matrix}] = \frac{1}{Δ x^{2}} [\begin{matrix} x_{N} - 2 x_{0} + x_{1} \\ x_{0} - 2 x_{1} + x_{2} \\ ⋮ \\ x_{N - 1} - 2 x_{N} + x_{0} \end{matrix}]

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

Это можно записать в матричной форме как , где

\frac{\partial}{\partial t} \vec{x} = \frac{1}{Δ x^{2}} A \vec{x}

$\frac{\partial}{\partial t}\vec{x} = \frac{1}{\Delta x^2} \mathbf{A} \vec{x}$

A = [\begin{matrix} - 2 & 1 & 0 & \dots & 0 & 1 \\ 1 & - 2 & 1 & 0 & \dots & 0 \\ ⋱ & ⋱ & ⋱ \\ ⋱ & ⋱ & ⋱ \\ 0 & \dots & 0 & 1 & - 2 & 1 \\ 1 & 0 & \dots & 0 & 1 & - 2 \end{matrix}] .

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

Конечно, нет необходимости создавать или хранить эту матрицу. Конечные различия должны быть вычислены на лету, обращаясь к первым и последним точкам, особенно по мере необходимости.

В качестве простого примера, следующие MATLAB скрипт решает с периодическими граничными условиями на области . Изготовленный раствор

\partial_{t} u = \partial_{x x} u + b (t, x)

$\partial_t u = \partial_{xx}u + b(t,x)$

x \in [- 1, 1)

$x\in[-1,1)$

, что означает

u_{Ref} (t, x) = \exp (- t) \cos (5 π x)

$u_\text{Ref}(t,x) = \exp(-t)\cos(5\pi x)$

. Я использовал прямую дискретизацию Эйлера по времени для простоты и вычислил решение как с матрицей, так и без нее. Результаты показаны ниже.

b (t, x) = (25 π^{2} - 1) \exp (- t) \cos (5 π 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

Участок начального состояния

Участок раствора при t = 0,5

Участок решения при t = 1,0

Участок решения при t = 2,0

Дуг Липински
источник

Отличное и простое решение !! в случае, если кому-то это нужно, здесь реализации в Python

ilciavo

x

$x$

@ bela83 Вы правы, что нет необходимости указывать что-либо большее, чем начальное условие. Это может привести к переопределению системы. Все, что вам нужно сделать, - это быть немного осторожнее рядом с конечными точками интервала, чтобы быть уверенным, что вы периодически правильно оборачиваете вещи. Есть много способов сделать это.

Дуг Липински

-1

В соответствии с этим вы должны наложить периодические граничные условия как:

U (0, t) = u (1, t) u_{x} (0, t) = u_{x} (1, t)

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

Одним из способов неявного выражения уравнения теплопроводности с использованием обратного Эйлера является

\frac{u^{n + 1} - u_{n}}{Δ t} = \frac{u_{i + 1}^{n + 1} - 2 u_{i}^{n + 1} + u_{i + 1}^{n + 1}}{Δ x^{2}}

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

Решение системы уравнений

[\begin{matrix} I - \frac{Δ t}{Δ x^{2}} A \end{matrix}] [\begin{matrix} u_{1}^{n + 1} \\ u_{1}^{n + 1} \\ ⋮ \\ u_{N}^{n + 1} \end{matrix}] = [\begin{matrix} u_{1}^{n} \\ u_{2}^{n} \\ ⋮ \\ u_{N}^{n} \end{matrix}]

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

A = [\begin{matrix} - 2 & 1 & 0 & 0 & 0 & \dots & 0 \\ 1 & - 2 & 1 & 0 & 0 & \dots & 0 \\ 0 & 1 & - 2 & 1 & 0 & \dots & 0 \\ 0 & 0 & 1 & - 2 & 1 & \dots & 0 \\ 0 & 0 & 0 & 1 & - 2 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & 0 & 1 & - 2 \end{matrix}]

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

u_{1} - u_{N} = 0 \frac{u_{2} - u_{0}}{2 Δ x} - \frac{u_{N + 1} - u_{N - 1}}{2 Δ x} = 0

$\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{matrix} 0 & 1 & 0 & 0 & \dots & 0 & - 1 & 0 \\ - 1 & 0 & 1 & 0 & \dots & 1 & 0 & - 1 \\ 0 & 0 \\ 0 & 0 \\ 0 & I - \frac{Δ t}{Δ x^{2}} A & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{matrix}] [\begin{matrix} u_{0}^{n + 1} \\ u_{1}^{n + 1} \\ u_{2}^{n + 1} \\ u_{N}^{n + 1} \\ u_{N + 1}^{n + 1} \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ u_{1}^{n} \\ u_{2}^{n} \\ u_{N}^{n} \end{matrix}]

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

ilciavo
источник

I don't understand the statement : "add two additional cells

u_{- 1}

$u_{-1}$ and

u_{N}

$u_N$ " because

u_{N}

$u_N$ is already a point in

] 0, 1 [

$]0,1[$ . I have in mind

x_{k} = \frac{k}{N + 1}

$x_k=\frac{k}{N+1}$ so that

x_{N} = \frac{N}{N + 1}

$x_N=\frac{N}{N+1}$ .

bela83

It's just a matter of indexing. You start with

N

$N$ cells (or points) from

u_{0}

$u_0$ to

u_{N - 1}

$u_{N-1}$ and you add two cells

u_{- 1}

$u_{-1}$ and

u_{N}

$u_{N}$ . If you have cells going from

u_{1}

$u_1$ to

u_{N}

$u_N$ then you need to add cells at

u_{0}

$u_0$ and

u_{N + 1}

$u_{N+1}$

ilciavo

OK then I don't understand the two more equations ! The first should be (with the notation from the question):

u_{0} = u_{N + 1}

$u_0=u_{N+1}$ right ? But I don't understand the second one: why not choose a centered difference approximation ? Finally, that makes

N + 1

$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

u_{x}

$u_x$ . The first comes from

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

$u(0,t)=u(1,t)$ and the second from

u_{x} (0, t) = u_{x} (1, t)

$u_x(0,t)=u_x(1,t)$

ilciavo

u1=uN adds an additional restriction(equation u1−uN=0) to your system

ilciavo