Решение ДУЧП с несогласованными начальными и граничными условиями

Обычно при решении дифференциальных уравнений в частных производных (ДУЧП) требуют, чтобы начальные и граничные условия были согласованными, т.е. чтобы начальные условия удовлетворяли граничным условиям.
Что представляет собой решение ДУЧП в случае, когда начальные и граничные условия не согласованы? Как можно вычислить решение уравнения в этом случае?

В качестве примера рассмотрим однородное уравнение теплопроводности

с однородными граничными условиями

которые не согласуются с начальным условием

Точное решение

Найдем решение этого уравнения. Как обычно, ищем решение в виде

Подставляя данную функцию в уравнение и разделяя переменные, получим

или

Отсюда

 

 Общее решение уравнения имеет вид ряда

Для нахождения постоянных C_n разложим начальную функцию в ряд Фурье:

Таким образом,

Построим график решения при помощи следующего кода в Wolfram Mathematica (также можно использовать Облако Wolfram):

tt = 0.1;

usol[x_, t_] :=
 Sum[2/(Pi*n)*(1 - (-1)^n)*Sin[Pi*n*x]*Exp[-Pi^2*n^2*t], {n, 1, 100}]

Plot3D[usol[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

Приближённое решение

Wolfram Mathematica

Воспользуемся стандартным решателем Wolfram Mathematica:

tt = 1;

s = NDSolve[{D[u[x, t], t] == D[u[x, t], x, x],
   u[x, 0] == 1,
   u[0, t] == 0, u[1, t] == 0},
  u, {t, 0, tt}, {x, 0, 1}];

usol[x_, t_] := Evaluate[u[x, t] /. s];

Plot3D[usol[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

Mathematica выдаёт предупреждение, что начальные и граничные условия не согласованы, но решение считает неверно.

MATLAB

Попробуем применить MATLAB для решения данного уравнения.

function a

tt = 0.1;

m = 0;
x = linspace(0,1,20);
t = linspace(0,tt,20);

sol = pdepe(m,@apde,@aic,@abc,x,t);
u = sol(:,:,1);

figure;
surf(x,t,u);
title('Numerical solution');
xlabel('Distance x');
ylabel('Time t');

figure;
plot(x,u(end,:));
title('Solution at t = 0.1');
xlabel('Distance x');
ylabel('u(x,0.1)');

% ----------------------------------------------------------------

function [c,f,s] = apde(x,t,u,DuDx)
c = 1;
f = DuDx;
s = 0;

% ----------------------------------------------------------------

function u0 = aic(x)
u0 = 1;
   
% ----------------------------------------------------------------

function [pl,ql,pr,qr] = abc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur;
qr = 0;

Решение вычисляется верно:

FreeFem++

Теперь применим FreeFem++:

func u0 = 1;
real T = 0.1;
int tn = 100;
real dt = T/(tn-1);

real h = 0.2;
int npts = 100;

mesh Th = square(100,10,[x,h*y]);

fespace Vh(Th,P1);

Vh u=u0, v, uold;
problem thermic(u,v) =
  int2d(Th)(u*v/dt + (dx(u) * dx(v) + dy(u) * dy(v)))
  - int2d(Th)(uold*v/dt) + on(2,4,u=0);

ofstream gnu("output.txt");
gnu.precision(10);

for (int i = 1; i <= npts; i++) {
  real xi = 1.0*(i-1)/(npts-1);
  gnu << u(xi,h/2) << endl;
}

real t = 0;
for (int tnum = 2; tnum <= tn; tnum++) {
  uold = u;
  thermic;
  t += dt;

  for (int i = 1; i <= npts; i++) {
    real xi = 1.0*(i-1)/(npts-1);
    gnu << u(xi,h/2) << endl;
  }
}

plot(u,fill=true,wait=true);
//plot(u,fill=true,wait=true,value=1);

Для отображения результатов вычислений используем следующий код в Wolfram Mathematica:

tt = 0.1;

tnn = 100;

npts = 100;

str = OpenRead["C:/ff_new/0/output.txt"];

Do[
 Do[inp[i, j] = Read[str, Number], {i, npts}];
 , {j, tnn}]

u = ListInterpolation[
   Table[inp[i, j], {i, npts}, {j, tnn}], {{0, 1}, {0, tt}},
   Method -> "Spline"];

Plot3D[u[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

Решение вычисляется правильно:

Отметим, что на границе вблизи начального момента времени вычисленный график показывает отрицательные значения.

Но при этом все числа, которые выдаёт FreeFem, неотрицательны. Такое поведение графика вызвано погрешностью интерполяции.

Программирование разностных схем

Теперь выясним, как вычисляют решение стандартные разностные схемы [1].

Явная схема

Простейшая явная схема имеет вид:

Схема условно устойчива, т.е. устойчива при  (условие Куранта).
Код в Wolfram Mathematica:

tt = 0.1;

tnn = 100;

xnn = 20;

tau = tt/tnn;

h = 1/xnn;

u0[x_] := 1;

Do[v[i, 0] = u0[h*i], {i, 0, xnn}]

Do[
 v[0, n] = 0; v[xnn, n] = 0;
 Do[
  v[i, n + 1] = v[i, n] + tau*(v[i - 1, n] - 2*v[i, n] + v[i + 1, n])
  , {i, 1, xnn - 1}]
 , {n, 0, tnn - 1}]

u = ListInterpolation[
   Table[v[i, j], {i, 0, xnn}, {j, 0, tnn}], {{0, 1}, {0, tt}}];

Plot3D[u[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

Шаги сетки выбраны так, чтобы выполнялось условие Куранта. Тем не менее, схема неправильно считает решение.

Решение практически постоянное. Это связано с тем, что вдали от граничных точек аппроксимация 2-й производной равна 0, так как начальное условие постоянно. В явной схеме значение на следующем временном слое зависит только от значений в соседних узлах на предыдущем слое. Поэтому на следующем слое значения остаются такими же, а должны убывать.
По всей видимости, условие устойчивости имеет место в случае согласованных начальных и граничных условий.

 

Неявная схема

Простейшая неявная схема имеет вид:

или

где .
Схема абсолютно устойчива.
Код в Wolfram Mathematica:

tt = 0.1;

tnn = 100;

xnn = 100;

tau = tt/tnn;

h = 1/xnn;

r = tau/h^2;

u0[x_] := 1;

Do[v[i, 0] = u0[h*i], {i, 0, xnn}]

Do[
 v[0, n + 1] = 0; v[xnn, n + 1] = 0;
 Do[a[i, j] = 0, {i, 1, xnn - 1}, {j, 1, xnn - 1}];
 Do[a[i, i] = 1 + 2*r, {i, 1, xnn - 1}];
 Do[a[i, i - 1] = -r, {i, 2, xnn - 1}];
 Do[a[i, i + 1] = -r, {i, 1, xnn - 2}];
 Do[b[i] = v[i, n], {i, 1, xnn - 1}];
 aa = Table[a[i, j], {i, 1, xnn - 1}, {j, 1, xnn - 1}];
 bb = Table[b[i], {i, 1, xnn - 1}];
 l = LinearSolve[N[aa], N[bb]];
 Do[v[i, n + 1] = l[[i]], {i, 1, xnn - 1}];
 , {n, 0, tnn - 1}]

u = ListInterpolation[
   Table[v[i, j], {i, 0, xnn}, {j, 0, tnn}], {{0, 1}, {0, tt}}];

Plot3D[u[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

Решение вычисляется правильно:

Выводы

Если начальные и граничные условия не согласованы, то решение всё равно можно найти в виде ряда Фурье, хотя оно не будет непрерывно в "угловых" точках пространственно-временного цилиндра (x=0,t=0), (x=1,t=0). Во всех остальных точках решение непрерывно.
Простейшая явная схема неправильно вычисляет решение, так как в явной схеме решение на следующем слое зависит только от значений в соседних узлах на предыдущем слое. Неявная схема правильно считает решение, так как в этом случае решение на следующем слое определяется как решение СЛАУ и зависит сразу от всего предыдущего слоя.
Wolfram Mathematica неправильно считает решение. MATLAB и FreeFem++ вычисляют решение правильно.

Литература

[1] Алексеев Г.В. Введение в численные методы решения дифференциальных уравнений. Уч. пос. Владивосток: Дальрыбвтуз (ТУ), 2002. - 140 с.

Дополнение (11.09.2014)

Wolfram Mathematica может решать уравнения с несогласованными начальными и граничными условиями, просто в граничных точках нужно начальное условие приравнять к граничным.

Код:

tt = 1;

s = NDSolve[{D[u[x, t], t] == D[u[x, t], x, x],
u[x, 0] == If[x == 0 || x == 1, 0, 1],
u[0, t] == 0,
u[1, t] == 0}, u, {t, 0, tt}, {x, 0, 1}];

usol[x_, t_] := Evaluate[u[x, t] /. s];

Plot3D[usol[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

Источник - Русскоязычная поддержка Wolfram Mathematica: http://vk.com/wall-1172233_16828

Дополнение (21.11.2015)

Спасибо анонимному комментатору, у меня была ошибка.

Код для явной схемы:

tt = 0.1;

tnn = 100;

xnn = 20;

tau = tt/tnn;

h = 1/xnn;

u0[x_] := 1;

Do[v[i, 0] = u0[h*i], {i, 0, xnn}]

Do[v[0, n] = 0; v[xnn, n] = 0;
 Do[v[i, n + 1] =
   v[i, n] + tau/h^2*(v[i - 1, n] - 2*v[i, n] + v[i + 1, n]), {i, 1,
   xnn - 1}], {n, 0, tnn - 1}]

u = ListInterpolation[
   Table[v[i, j], {i, 0, xnn}, {j, 0, tnn}], {{0, 1}, {0, tt}}];

Plot3D[u[x, t], {x, 0, 1}, {t, 0, tt}, PlotRange -> All]

График решения: