Решающая система PDE

Открыть Live Script

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

Рассмотрите систему PDE

$\frac{\partial u_{1}}{\partial t} = 0.024 \frac{\partial^{2} u_{1}}{\partial x^{2}} - F (u_{1} - u_{2}),$

$\frac{\partial u_{2}}{\partial t} = 0.170 \frac{\partial^{2} u_{2}}{\partial x^{2}} + F (u_{1} - u_{2}) .$

(Функция $F (y) = e^{5.73 y} - e^{- 11.46 y}$ используется в качестве стенограммы.)

Уравнение удерживается на интервале $0 \leq x \leq 1$ для раз $t \geq 0$ . Начальные условия

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

$u_{2} (x, 0) = 0 .$

Граничные условия

$\begin{array}{l} \frac{\partial}{\partial x} u_{1} (0, t) = 0, \\ u_{2} (0, t) = 0, \\ \frac{\partial}{\partial x} u_{2} (1, t) = 0, \\ u_{1} (1, t) = 1 . \end{array}$

Чтобы решить это уравнение в MATLAB, необходимо кодировать уравнение, начальные условия и граничные условия, затем выбрать подходящий mesh решения перед вызовом решателя pdepe. Можно либо включить необходимые функции в качестве локальных функций в конце файла (как это сделано здесь), либо сохранить их как отдельные, именованные файлы в директории по пути MATLAB.

Кодовое уравнение

Прежде чем вы сможете кодировать уравнение, необходимо убедиться, что это в том виде, в котором pdepe решатель ожидает:

$c (x, t, u, \frac{\partial u}{\partial x}) \frac{\partial u}{\partial t} = x^{- m} \frac{\partial}{\partial x} (x^{m} f (x, t, u, \frac{\partial u}{\partial x})) + s (x, t, u, \frac{\partial u}{\partial x}) .$

В этой форме коэффициенты УЧП являются матричными, и уравнение становится

$[\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}] \frac{\partial}{\partial t} [\begin{matrix} u_{1} \\ u_{2} \end{matrix}] = \frac{\partial}{\partial x} [\begin{matrix} 0.024 \frac{\partial u_{1}}{\partial x} \\ 0.170 \frac{\partial u_{2}}{\partial x} \end{matrix}] + [\begin{matrix} - F (u_{1} - u_{2}) \\ F (u_{1} - u_{2}) \end{matrix}] .$

Таким образом, значения коэффициентов в уравнении

$m = 0$

$c (x, t, u, \frac{\partial u}{\partial x}) = [\begin{matrix} 1 \\ 1 \end{matrix}]$ (только по диагонали)

$f (x, t, u, \frac{\partial u}{\partial x}) = [\begin{matrix} 0.024 \frac{\partial u_{1}}{\partial x} \\ 0.170 \frac{\partial u_{2}}{\partial x} \end{matrix}]$

$s (x, t, u, \frac{\partial u}{\partial x}) = [\begin{matrix} - F (u_{1} - u_{2}) \\ F (u_{1} - u_{2}) \end{matrix}]$

Теперь можно создать функцию, чтобы кодировать уравнение. Функция должна иметь подпись [c,f,s] = pdefun(x,t,u,dudx):

x - независимая пространственная переменная.
t - независимая временная переменная.
u - зависимая переменная, дифференцируемая по отношению к x и t. Это двухэлементный вектор, где u(1) является $u_{1} (x, t)$ и u(2) является $u_{2} (x, t)$ .
dudx - частная пространственная производная $\partial u / \partial x$ . Это двухэлементный вектор, где dudx(1) является $\partial u_{1} / \partial x$ и dudx(2) является $\partial u_{2} / \partial x$ .
Выходные выходы c, f, и s соответствуют коэффициентам в стандартной форме уравнения УЧП, ожидаемым по pdepe.

В результате уравнения в этом примере могут быть представлены функцией:

function [c,f,s] = pdefun(x,t,u,dudx)
c = [1; 1];
f = [0.024; 0.17] .* dudx;
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];
end

(Примечание. Все функции включены в качестве локальных функций в конце примера.)

Начальные условия кода

Затем напишите функцию, которая возвращает начальное условие. Начальное условие применяется в первый раз и предоставляет значение $u (x, t_{0})$ для любого значения x. Количество начальных условий должно равняться количеству уравнений, поэтому для этой задачи существует два начальных условия. Используйте сигнатуру функции u0 = pdeic(x) для записи функции.

Начальные условия

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

$u_{2} (x, 0) = 0 .$

Соответствующая функция является

function u0 = pdeic(x)
u0 = [1; 0];
end

Граничные условия кода

Теперь напишите функцию, которая оценивает граничные условия

$\begin{array}{l} \frac{\partial}{\partial x} u_{1} (0, t) = 0, \\ u_{2} (0, t) = 0, \\ \frac{\partial}{\partial x} u_{2} (1, t) = 0, \\ u_{1} (1, t) = 1 . \end{array}$

Для задач, поставленных на интервале $a \leq x \leq b$ , граничные условия применяются для всех $t$ и либо $x = a$ или $x = b$ . Стандартная форма для граничных условий, ожидаемых решателем,

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

Записанные в этой форме граничные условия для частных производных $u$ должны быть выражены в терминах потока $f (x, t, u, \frac{\partial u}{\partial x})$ . Таким образом, граничные условия для этой задачи

Для $x = 0$ , уравнение

$[\begin{matrix} 0 \\ u_{2} \end{matrix}] + [\begin{matrix} 1 \\ 0 \end{matrix}] \cdot [\begin{matrix} 0.024 \frac{\partial u_{1}}{\partial x} \\ 0.170 \frac{\partial u_{2}}{\partial x} \end{matrix}] = 0 .$

Коэффициенты:

$p_{L} (x, t, u) = [\begin{matrix} 0 \\ u_{2} \end{matrix}],$

$q_{L} (x, t) = [\begin{matrix} 1 \\ 0 \end{matrix}] .$

Аналогично, для $x = 1$ уравнение

$[\begin{matrix} u_{1} - 1 \\ 0 \end{matrix}] + [\begin{matrix} 0 \\ 1 \end{matrix}] \cdot [\begin{matrix} 0.024 \frac{\partial u_{1}}{\partial x} \\ 0.170 \frac{\partial u_{2}}{\partial x} \end{matrix}] = 0 .$

Коэффициенты:

$p_{R} (x, t, u) = [\begin{matrix} u_{1} - 1 \\ 0 \end{matrix}],$

$q_{R} (x, t) = [\begin{matrix} 0 \\ 1 \end{matrix}] .$

Граничная функция должна использовать сигнатуру функции [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t):

Входы xl и ul соответствуют $u$ и $x$ для левого контура.
Входы xr и ur соответствуют $u$ и $x$ для правого контура.
t - независимая временная переменная.
Выходные выходы pl и ql соответствуют $p_{L} (x, t, u)$ и $q_{L} (x, t)$ для левого контура ( $x = 0$ для этой задачи).
Выходные выходы pr и qr соответствуют $p_{R} (x, t, u)$ и $q_{R} (x, t)$ для правого контура ( $x = 1$ для этой задачи).

Граничные условия в этом примере представлены функцией:

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];
end

Выберите Mesh решения

Решение этой проблемы изменяется быстро, когда $t$ маленькая. Хотя pdepe выбирает временной шаг, подходящий для разрешения резких изменений, чтобы увидеть поведение на выходных графиках, необходимо выбрать соответствующее время выхода. Для пространственного mesh в решении есть граничные слои на обоих концах $0 \leq x \leq 1$ , поэтому необходимо задать точки mesh там, чтобы разрешить резкие изменения.

x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

Решение уравнения

Наконец, решите уравнение, используя симметрию $m$ , уравнение УЧП, начальные условия, граничные условия и сетки для $x$ и $t$ .

m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);

pdepe возвращает решение в трехмерный массив sol, где sol(i,j,k) аппроксимирует kпервый компонент решения $u_{k}$ оценивается в t(i) и x(j). Извлеките каждый компонент решения в отдельную переменную.

u1 = sol(:,:,1);
u2 = sol(:,:,2);

Решение для построения графика

Создайте объемные поверхностные диаграммы решений для $u_{1}$ и $u_{2}$ нанесенный на график в выбранных точках mesh для $x$ и $t$ .

surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

Figure contains an axes. The axes with title u_1(x,t) contains an object of type surface.

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

Figure contains an axes. The axes with title u_2(x,t) contains an object of type surface.

Локальные функции

Здесь перечислены локальные вспомогательные функции, которые решатель PDE pdepe вызывается для вычисления решения. Также можно сохранить эти функции как собственные файлы в директории по пути MATLAB.

function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve
c = [1; 1];
f = [0.024; 0.17] .* dudx;
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];
end
% ---------------------------------------------
function u0 = pdeic(x) % Initial Conditions
u0 = [1; 0];
end
% ---------------------------------------------
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];
end
% ---------------------------------------------

См. также

pdepe

Документация