Граничные условия для уравнения переноса, дискретизированного методом конечных разностей

Я пытаюсь найти некоторые ресурсы, которые помогут объяснить, как выбирать граничные условия при использовании методов конечных разностей для решения PDE.

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

Общие правила, регулирующие стабильность при наличии границ, слишком сложны для вводного текста; они требуют сложных математических механизмов

(А. Изерлес. Первый курс по численному анализу дифференциальных уравнений)

Например, при попытке реализовать двухэтапный метод чехарды для уравнения адвекции:

$u_i^{n+1} = u_i^{n-1} + \mu (u_{i+1}^n - u_{i-1}^n)$

используя MATLAB

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);
    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

Решение ведет себя хорошо, пока не достигнет границы, когда оно внезапно начинает вести себя плохо.

Где я могу узнать, как обрабатывать граничные условия, как это?

Ответ Слоеда очень тщательный и правильный. Я просто хотел добавить несколько моментов, чтобы было легче понять.

По сути, любое волновое уравнение имеет собственную скорость и направление волны. Для одномерного волнового уравнения: скорость волны - это постоянная a, которая определяет не только скорость, с которой информация распространяется в области, но и ее направление. Если a > 0 , информация идет слева направо, а если a < 0, это наоборот.

$$u_t + a u_x = 0$$ $a$$a>0$$a<0$

Для метода скачкообразной лягушки, когда вы дискретизируете уравнения, вы получаете: или: u n i =u n - 2 i +μ(u n - 1 i + 1 -u n - 1 i - 1 ),гдеμ=-aΔt/Δx. В вашем случаеμ>

$$\frac{u_i^{n} - u_i^{n-2}}{2 \Delta t} + a \frac{u_{i+1}^{n-1} - u_{i-1}^{n-1}}{2 \Delta x} = 0$$ $$u_i^{n} = u_i^{n-2} + \mu (u_{i+1}^{n-1} - u_{i-1}^{n-1})$$ $\mu = -a \Delta t/\Delta x$$\mu > 0$что переводится как волна, идущая влево. Теперь, если вы подумаете об этом, для волны, которая движется влево, понадобится только граничное условие на правой границе, поскольку все значения слева обновляются через своих правых соседей. Фактически, указание любого значения на левой границе не согласуется с природой проблемы. В некоторых методах, например при простом ветре, об этом заботятся автоматически, поскольку в схеме также используются только правильные соседи. В других методах, таких как скачкообразная лягушка, вы должны указать «правильное» значение.

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

Общий ответ
Ваша проблема в том, что вы вообще не устанавливаете (или даже не указываете) граничные условия - ваша численная задача плохо определена.

Как правило, существует два возможных способа задания граничных условий:

1. $u_{0}$$u_{101}$
2. Измените числовой трафарет, чтобы он использовал только внутреннюю информацию на границе.

Какой путь вы выберете, зависит от физики вашей проблемы. Для задач типа волнового уравнения обычно определяют собственные значения потока Якобиана, чтобы решить, нужны ли внешние граничные условия, или следует использовать внутреннее решение (этот метод обычно называют «намоткой»).

---

$u_{i-1}^n$$u_{i+1}^n$$i$$n+1$$i = 1$$u_{0}^n$$u_{100}^{n+1}$$u_{101}^n$

$u_{1}^n$$u_{100}^n$

Вы можете найти измененную версию вашего исходного кода ниже:

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    %u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);

    % Apply the numerical stencil to all interior points
    for j = 2:M-1
        u(j,i) = u(j,i-2) + mu*(u(j+1,i-1) - u(j-1,i-1));
    end

    % Set the boundary values by interpolating linearly from the interior
    u(1,i) = 2*u(2,i) - u(3,i);
    u(M,i) = 2*u(M-1,i) - u(M-2,i);

    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

Итак, я рассмотрел это более подробно, и кажется, что это (по крайней мере, в основных случаях, с которыми я работаю) зависит от групповой скорости метода.

Метод чехарды (например):

$u_{i}^{n+1} = u_{i}^{n-1} + \mu (u_{i+1}^n - u_{i-1}^n)$

$u_k^n = e^{i( \zeta k \Delta x + \omega(\zeta) n \Delta t)}$

$e^{2 i \omega \Delta t} = 1 + \mu e^{i\omega \Delta t} ( e^{i \zeta \Delta x} - e^{-i \zeta \Delta x} )$

$\sin(\omega \Delta t) = \mu \sin (\zeta \Delta x)$

$\frac{d\omega}{d\zeta} = \frac{\cos (\zeta \Delta x)}{\sqrt{1-\mu^2 sin^2(\zeta \Delta x)}} \in [-1,1]$

Теперь нам нужно выяснить групповую скорость граничных условий:

Мой метод :$u_1^{n+2} = u_1^n + \mu u_2^{n+1}$

Мы можем вычислить граничную групповую скорость следующим образом:

$2 i \sin(\omega \Delta t) = \mu e^{i \zeta \Delta x}$

поэтому, чтобы найти некоторые групповые скорости, которые позволяют границы, нам нужно найти:

$\omega = c \zeta$

$\cos(\zeta \Delta x) = 0, \mu \sin(\zeta \Delta x) = 2 \sin (\zeta c \Delta t)$

$\zeta = \frac{\pi}{2 \Delta x}$ would give $\mu = 2 \sin (\frac{c \mu \pi}{2})$ for which a solution for $c \in [-1,1]$ will exist. (For most choices of $\mu$ at least)

---

The solution which I've found in the literature is to take $u_0^{n+1} = u_1^{n}$ since this has a boundary wave number which lies outside $[-1,1]$.

I've still quite quite a bit more to read up about this before I understand it completely. I think the key words I'm looking for are GKS theory.

Source for all this A Iserles Part III notes

---

A clearer calculation of what I've done can be found here: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Guys I am very new to this site. Maybe this isn't the place to ask, but please forgive me as I am very new here :) I am having an extremely similar problem, the only difference being the starting function which, in my case, is a cosine wave. My code is this: clear all; clc; close all;

M = 1000; N = 2100;

mu = 0.5;

c = [mu 0 -mu]; f = @(x)1- cos(20*pi*x-0.025).^2; u = zeros (M, N); x=0:(1/M):0.05; u(1:length(x),1) = f(x); u(1:length(x),2) = f(x - mu/(M)); x=linspace(0,1,M);

for i = 3:N hold off;

% Apply the numerical stencil to all interior points
for j = 2:M-1
    u(j,i) = u(j,i-2) - mu*(u(j+1,i-1) - u(j-1,i-1));
end

% Set the boundary values by interpolating linearly from the interior
u(M,i) =  2*u(M-1,i-1) - u(M-2,i-1);

plot(x, u(:,i)); axis( [ 0 1.5 -0.5 2] ) drawnow; %pause end

There is already this code here, but for some reason, probably related to the cosine wave, my code fails :/ any help would be appreciated :) thanks!