Обычно при решении дифференциальных уравнений в частных производных (ДУЧП) требуют, чтобы начальные и граничные условия были согласованными, т.е. чтобы начальные условия удовлетворяли граничным условиям.
Что представляет собой решение ДУЧП в случае, когда начальные и граничные условия не согласованы? Как можно вычислить решение уравнения в этом случае?
В качестве примера рассмотрим однородное уравнение теплопроводности
с однородными граничными условиями
которые не согласуются с начальным условием
Подставляя данную функцию в уравнение и разделяя переменные, получим
или
Отсюда
Для нахождения постоянных Cn разложим начальную функцию в ряд Фурье:
Таким образом,
Построим график решения при помощи следующего кода в 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]
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 выдаёт предупреждение, что начальные и граничные условия не согласованы, но решение считает неверно.
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;
Решение вычисляется верно:
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]
Решение вычисляется правильно:
Отметим, что на границе вблизи начального момента времени вычисленный график показывает отрицательные значения.
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]
Что представляет собой решение ДУЧП в случае, когда начальные и граничные условия не согласованы? Как можно вычислить решение уравнения в этом случае?
В качестве примера рассмотрим однородное уравнение теплопроводности
с однородными граничными условиями
которые не согласуются с начальным условием
Точное решение
Найдем решение этого уравнения. Как обычно, ищем решение в видеПодставляя данную функцию в уравнение и разделяя переменные, получим
или
Отсюда
Общее решение уравнения имеет вид ряда
Таким образом,
Построим график решения при помощи следующего кода в 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:
![](http://2.bp.blogspot.com/-dDz5Cst1PUo/U-RHkFy3D_I/AAAAAAAAA0Q/9xEeO4shOwU/s1600/12.png)
Код в 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]
Решение вычисляется правильно:
Простейшая явная схема неправильно вычисляет решение, так как в явной схеме решение на следующем слое зависит только от значений в соседних узлах на предыдущем слое. Неявная схема правильно считает решение, так как в этом случае решение на следующем слое определяется как решение СЛАУ и зависит сразу от всего предыдущего слоя.
Wolfram Mathematica неправильно считает решение. MATLAB и FreeFem++ вычисляют решение правильно.
Решение практически постоянное. Это связано с тем, что вдали от граничных точек аппроксимация 2-й производной равна 0, так как начальное условие постоянно. В явной схеме значение на следующем временном слое зависит только от значений в соседних узлах на предыдущем слое. Поэтому на следующем слое значения остаются такими же, а должны убывать.
По всей видимости, условие устойчивости имеет место в случае согласованных начальных и граничных условий.
Неявная схема
Простейшая неявная схема имеет вид:или
где
![](http://4.bp.blogspot.com/-8ARDjD_1NRA/U-RKd41QmoI/AAAAAAAAA00/TmWJLrqRqG0/s1600/15.png)
Схема абсолютно устойчива.
Код в 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]
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
Спасибо анонимному комментатору, у меня была ошибка.
Код для явной схемы:
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]
График решения:Дополнение (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]
Комментарии
Do[
v[i, n + 1] = v[i, n] + tau / h^2 * (v[i - 1, n] - 2*v[i, n] + v[i + 1, n])
u[0, t] == t Exp[-t], Derivative[1, 0][u][1, t] == 0}
Добрый день. Не подскажете, как в Wolfram Mathematica модифицировать такое уравнение для его корректного решения?
(1) u_t = u_{xx},
(2) u|_{t=0} = sin(pi*x),
(3) u|_{x=0} = t*e^(-t), u_x|_{x=1} = 0.
Код для вычислений:
s = NDSolve[{D[u[x, t], t] == D[u[x, t], x, x],
u[x, 0] == Sin[\[Pi] x], u[0, t] == t Exp[-t],
Derivative[1, 0][u][1, t] == 0}, u, {x, 0, 1}, {t, 0, 100}];
usol[x_, t_] := Evaluate[u[x, t] /. s];
Plot3D[usol[x, t], {x, 0, 1}, {t, 0, 100}]
Wolfram Mathematica выдаёт предупреждение, что начальные и граничные условия не согласованы (начальная функция не удовлетворяет правому краевому условию).
Для проверки построим график решения при x = 0:
Plot[usol[0, t], {t, 0, 5}]
Видно, что левое краевое условие выполняется.
Построим график производной при x = 1:
Plot[Derivative[1, 0][usol][1, t], {t, 0, 100}]
Правое краевое условие не выполняется и равно -pi.
Построим график при t = 100:
Plot[usol[x, 100], {x, 0, 1}]
Получим функцию -pi*x.
В то же время решение задачи (1)-(3) должно стремиться к 0.
Если взять в качестве начальной функции, например, 1-cos(pi*x), то краевые условия выполняются, и решение, вероятно, посчитано верно.
Но если взять начальную функцию 1-cos(2*pi*x), то пакет почему-то выдаёт предупреждение и считает неверно, хотя начальная функция удовлетворяет краевым условиям.
s = NDSolve[{D[u[x, t], t] == D[u[x, t], x, x], u[x, 0] == Sin[Pi*x],
u[0, t] == t Exp[-t], Derivative[1, 0][u][1, t] == 0},
u, {x, 0, 1}, {t, 0, 10},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> False}];
При этом начальное условие как-то подстраивается под заданные краевые условия, и решение, которое выдаёт пакет, теперь удовлетворяет краевым условиям.
https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html
Также метод, который использует Mathematica, описан в статье Knapp R. A method of lines framework in Mathematica (https://goo.gl/m48pqZ)
Do[ u[i, j, 0] = v0[h*i, h*j], {i, 0, nn}, {j, 0, nn}]
Do[u[0, 0, k] = 0; u[nn, nn, k] = 0;
Do[u[i, j, k + 1] =
u[i, j, k] +
tau*(((u[i - 1, j, k] - 2*u[i, j, k] + u[i + 1, j, k])/
h^2) + ((u[i, j - 1, k] - 2*u[i, j, k] + u[i, j + 1, k])/
h^2)), {i, 1, nn - 1}, {j, 1, nn - 1}], {k, 0, mm - 1}]
не подскажите, как можно вывести результат? Пробовал через Annimate, пробовал через
z = ListInterpolation[
Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}, {k, 0, mm}], {{0, 1}, {0,
1}, {0, 1}}];
g[k_] := Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}];
p = Table[
ListPlot3D[g[k], PlotRange -> {{0, nn}, {0, mm}, {0, 6}}], {k, 0,
mm, 100}];
но не выводит. Буду признателен, если поможете.
Но весь процесс занимает очень много времени. Не подскажете, как можно ускорить вычисления?
Код такой:
LL = 1;
nn = 10;
mm = 1000;
h = 1/nn;
tau = LL/mm;
v0[x_, y_] := x;
Do[ u[i, j, 0] = v0[h*i, h*j], {i, 0, nn}, {j, 0, nn}]
Do[u[0, j, k] = 0; u[nn, j, k] = 0; u[i, 0, k] = 0; u[i, nn, k] = 0;
Do[u[i, j, k + 1] =
u[i, j, k] +
tau*(((u[i - 1, j, k] - 2*u[i, j, k] + u[i + 1, j, k])/
h^2) + ((u[i, j - 1, k] - 2*u[i, j, k] + u[i, j + 1, k])/
h^2)), {i, 1, nn - 1}, {j, 1, nn - 1}], {i, 1, nn - 1}, {j,
1, nn - 1}, {k, 0, mm - 1}]
g[k_] := Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}];
p = Table[
ListPlot3D[g[k], PlotRange -> {{0, nn}, {0, nn}, {0, 6}}], {k, 0,
mm, 100}]
вроде все делаю самым простым образом. Процесс идет, но крайне медленно.
даже с
LL = 1;
nn = 5;
mm = 100;
nn = 5;
mm = 200;
h = 1/nn;
tau = 1/mm;
v0[x_, y_] := 1;
Table[{x[i] = i*h, y[j] = j*h}, {i, 0, nn}, {j, 0, nn}];
Table[t[k] = k*t, {k, 0, mm}];
Table [u[i, j, 0] = v0[x[i], y[j]], {i, 0, nn}, {j, 0, nn}];
Table[u[0, j, k] = 0, {j, 0, nn}, {k, 0, mm}];
Table[u[nn, j, k] = 0, {j, 0, nn}, {k, 0, mm}];
Table[u[i, 0, k] = 0, {i, 1, nn - 1}, {j, 0, mm}];
Table[u[i, nn, k] = 0, {i, 1, nn - 1}, {k, 0, mm}];
Do[u[i, j, k + 1] =
u[i, j, k] +
tau*(((u[i - 1, j, k] - 2*u[i, j, k] + u[i + 1, j, k])/
h^2) + ((u[i, j - 1, k] - 2*u[i, j, k] + u[i, j + 1, k])/
h^2)), {i, 1, nn - 1}, {j, 1, nn - 1}, {k, 0, mm - 1}]
g[k_] := Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}];
p = Table[
ListPlot3D[g[k], PlotRange -> {{0, nn}, {0, nn}, {0, 6}}], {k, 0,
mm, 20}]
История та же.
Код:
d = 0.01;
a = b = 0.1;
nn = 10;
mm = 1000;
gamma = 0.35;
r = 1.5;
m = 4;
NPP0 = 65.766;
сn = 380;
tau = 19.2;
t = 1;
alpha = 720;
Nx = Ny = 100;
h = nn/Nx;
q = t/mm;
b0[x_, y_] := NPP0*tau;
Do[B[i, j, 0] = b0[h*i, h*j], {i, 0, Nx}, {j, 0, Ny}]
v0[x_, y_] := 461;
Do[с[i, j, 0] = v0[h*i, h*j], {i, 0, Nx}, {j, 0, Ny}]
Do[
Do[
с[Nx, j, k] = c[Nx - 1, j, k] + alpha*h;
с[1, j, k] = (c[0, j, k] + alpha*h), {j, 0, Ny}];
Do[
с[i, Ny, k] = (c[i, Ny - 1, k] + alpha*h);
с[i, 1, k] = (c[i, 0, k] + alpha*h), {i, 0, Nx}];
Do[
B[i, j, k + 1] =
B[i, j, k] +
q*(NPP0*(1 + gamma*Log[с[i, j, k]/сn]) - B[i, j, k]/tau);
с[i, j, k + 1] =
с[i, j, k] + q*(r*с[i, j, k] B[i, j, k])/(1 + m*с[i, j, k]) +
a*(с[i + 1, j, k] - с[i, j, k])/h +
b*(с[i, j + 1, k] - с[i, j, k])/h +
d*((с[i + 1, j, k] - 2*с[i, j, k] + с[i - 1, j, k])/(2*
h^2) + (с[i, j + 1, k] - 2*с[i, j, k] + с[i, j - 1, k])/(2*
h^2)), {i, 1, Nx - 1}, {j, 1, Ny - 1}], {k, 0, mm}]
g[k_] := Table[с[i, j, k], {i, 0, Nx}, {j, 0, Ny}];
p = Table[ListPlot3D[g[k], PlotRange -> All], {k, 0, mm, 100}]
Отправить комментарий