Решите систему УЧП

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

Рассмотрите систему УЧП

$\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$ .
Выходные параметры cF, и 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

Выберите Solution 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компонент th решения $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')

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

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

Перечисленный здесь локальные функции помощника что решатель УЧП 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

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