Система решения ПДЭ

Открыть сценарий в реальном времени

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

Рассмотрим систему ПДЭ

$\frac{_{}}{} \frac{^{}_{}}{^{}} ∂u1∂t=0.024∂2u1∂x2-F (_{}_{} u1-u2$ ),

$\frac{_{}}{} \frac{^{}_{}}{^{}} ∂u2∂t=0.170∂2u2∂x2+F (_{}_{} u1-u2$ ).

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

Уравнение удерживается на интервале, $0≤x≤1$ для времени $t≥0$ . Начальные условия:

$_{u1} (x, 0)$ = 1,

$_{u2} (x, 0)$ = 0.

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

$\begin{array}{l} \frac{}{}_{∂∂xu1} (0, t) \\ {= 0, u2}_{} (0, t) \\ \frac{}{}_{} =0,∂∂xu2 (1, t \\ {) = 0,}_{} u1 (1, t) \end{array}$ = 1.

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

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

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

$c (\frac{}{} x,t,u,\partialu\partialx \frac{)}{}^{} \frac{}{} {∂u∂t=x-m∂∂x}^{} (xmf (\frac{}{} x,t,u,∂u∂x)) + s \frac{(}{}$ x,t,u,∂u∂x).

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

$[\begin{matrix} \end{matrix} \frac{}{} \begin{matrix} _{} \\ _{} \end{matrix} \frac{}{} \begin{matrix} \frac{_{}}{} \\ \frac{_{}}{} \end{matrix} 1001]\partial\partialt[u1u2]=\partial\partialx[0.024\partialu1\partialx0.170\partialu2\partialx]+ \begin{matrix} [- F_{(}_{} \\ u1-u2)_{F} (_{} \end{matrix}$ u1-u2)].

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

$m =$ 0

$c (\frac{}{} x,t,u,\partialu\partialx) \begin{matrix} = \end{matrix} [$ 11] (только диагональные значения)

$f (\frac{}{} x,t,u,\partialu\partialx) \begin{matrix} \frac{_{}}{} \\ \frac{_{}}{} \end{matrix}$ =[0.024∂u1∂x0.170∂u2∂x]

$s (\frac{}{} x,t,u,\partialu\partialx) \begin{matrix} = [-_{F} (_{} \\ {u1-u2}_{)} F_{} ( \end{matrix}$ u1-u2)]

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

x является независимой пространственной переменной.
t - независимая переменная времени.
u является зависимой переменной, дифференцируемой по отношению к x и t. Это двухэлементный вектор, где u(1) является $_{u1} (x,$ t) иu(2) является $_{u2} (x,$ t).
dudx - частная пространственная производная $∂u/∂x$ . Это двухэлементный вектор, где dudx(1) является $_{} ∂u1/∂x$ и dudx(2) является $_{} ∂u2/∂x$ .
Продукция c, f, и s соответствуют коэффициентам в стандартной форме уравнения PDE, ожидаемым 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,_{}$ t0) для любого значения x. Число начальных условий должно равняться числу уравнений, поэтому для этой задачи существует два начальных условия. Использовать подпись функцииu0 = pdeic(x) для записи функции.

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

$_{u1} (x, 0)$ = 1,

$_{u2} (x, 0)$ = 0.

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

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

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

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

$\begin{array}{l} \frac{}{}_{∂∂xu1} (0, t) \\ {= 0, u2}_{} (0, t) \\ \frac{}{}_{} =0,∂∂xu2 (1, t \\ {) = 0,}_{} u1 (1, t) \end{array}$ = 1.

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

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

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

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

$[\begin{matrix} _{} \end{matrix} \begin{matrix} \end{matrix} \begin{matrix} \frac{_{}}{} \\ \frac{_{}}{} \end{matrix} 0u2]+[10]⋅[0.024∂u1∂x0.170∂u2∂x]=0 .$

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

$_{pL} (x, t, u \begin{matrix} ) \\ _{=} \end{matrix} [$ 0u2],

$_{qL} (x, t) \begin{matrix} = \end{matrix} [$ 10].

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

$[\begin{matrix} _{} \end{matrix} \begin{matrix} \end{matrix} \begin{matrix} \frac{_{}}{} \\ \frac{_{}}{} \end{matrix} u1-10]+[01]⋅[0.024∂u1∂x0.170∂u2∂x]=0 .$

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

$_{pR} (x, t, u \begin{matrix} )_{} = \\ [ \end{matrix}$ u1-10],

$_{qR} (x, t) \begin{matrix} = \end{matrix} [$ 01].

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

Исходные данные xl и ul соответствуют $u$ и $x$ для левой границы.
Исходные данные xr и ur соответствуют $u$ и $x$ для правой границы.
t - независимая переменная времени.
Продукция pl и ql соответствуют $_{pL} (x, t,$ u) $_{и} qL ($ x, t) для левой $границы$ (x = 0 для этой задачи).
Продукция pr и qr соответствуют $_{pR} (x, t,$ u) $_{и} qR ($ 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

Выбрать сетку решения

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

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$ , уравнение PDE, начальные условия, граничные условия и сетки для $x$ и $t$ .

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

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

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

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

Создайте графики поверхности решений для $_{u1}$ и $_{u2}$ , выводимых на печать в выбранных точках сетки для $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

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