Решение дифференциальных уравнений в частных производных

В дифференциальном уравнении в частных производных (PDE) решаемая функция зависит от нескольких переменных, и дифференциальное уравнение может включать в себя частные производные, взятые относительно каждой из переменных. Дифференциальные уравнения в частных производных полезны для моделирования волн, теплового потока, дисперсии жидкости и других явлений с пространственным поведением, которое изменяется с течением времени.

Какие типы PDE можно решить с помощью MATLAB?

Решатель MATLAB ® PDEpdepe решает начально-граничные задачи для систем PDE в одной пространственной переменной x и времени t. Их можно рассматривать как ODE одной переменной, которые также изменяются по времени.

pdepe использует неформальную классификацию для 1-D уравнений, которые она решает:

Уравнения с производной по времени параболические. Примером может служить тепловое уравнение $\frac{}{} \frac{^{}}{^{∂u∂t=∂2u∂x2}}$ .
Уравнения без производной по времени являются эллиптическими. Примером может служить уравнение Лапласа $\frac{^{}}{^{}} ∂2u∂x2=0$ .

pdepe требует по меньшей мере одного параболического уравнения в системе. Другими словами, по меньшей мере одно уравнение в системе должно включать производную по времени.

pdepe также решает определенные задачи 2-D и 3-D, которые сводятся к проблемам 1-D из-за угловой симметрии (см. описание аргумента для константы симметрии m для получения дополнительной информации.

Дифференциальное уравнение Toolbox™ в частных производных расширяет эту функциональность до обобщенных задач в 2-D и 3-D с граничными условиями Дирихле и Неймана.

Решение 1-D PDE

1-D PDE включает в себя функцию u (x, t), которая зависит от времени t, и одну пространственную переменную x. pdepe решает системы 1-D параболического и овального PDEs формы

$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 удерживаются для t0 ≤ t ≤ tf и a ≤ x ≤ b.
Пространственный интервал [a, b] должен быть конечным.
m может быть 0, 1 или 2, что соответствует симметрии перекрытия, цилиндрической или сферической симметрии соответственно. Если m > 0, то ≥ 0 также должен быть удержан.
Коэффициент $f (\frac{}{} x,t,u,\partialu\partialx$ ) - член потока, а $s (\frac{}{}$ x,t,u,∂u∂x) - член источника.
Член потока должен зависеть от ∂u/∂x частных производных.

Связь частных производных по времени ограничивается умножением на диагональную матрицу $c (\frac{}{} x,t,u,\partialu\partialx$ ). Диагональные элементы этой матрицы являются либо нулевыми, либо положительными. Нулевой элемент соответствует эллиптическому уравнению, а любой другой элемент - параболическому уравнению. Должно быть хотя бы одно параболическое уравнение. Элемент c, который соответствует параболическому уравнению, может исчезнуть при изолированных значениях x, если они являются точками сетки (точками, где вычисляется решение). Разрывы в c и s из-за интерфейсов материала допускаются при условии, что точка сетки размещается на каждом интерфейсе.

Процесс решения

Решение PDE с помощью pdepe, необходимо определить коэффициенты уравнений для c, f и s, исходные условия, поведение решения на границах и сетку точек для вычисления решения. Вызов функции sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) использует эту информацию для вычисления решения по указанной сетке:

m - константа симметрии.
pdefun определяет решаемые уравнения.
icfun определяет начальные условия.
bcfun определяет граничные условия.
xmesh является вектором пространственных значений для x.
tspan - вектор значений времени для t.

Вместе, xmesh и tspan векторы образуют сетку 2-D, которая pdepe оценивает решение на.

Уравнения

Необходимо выразить PDE в стандартной форме, ожидаемой pdepe. Записанные в этой форме, можно считать значения коэффициентов c, f и s.

В MATLAB можно кодировать уравнения функцией формы.

function [c,f,s] = pdefun(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end

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

\frac{}{} \frac{^{}}{^{∂u∂t=∂2u∂x2}}

. Если существует несколько уравнений, то c, f и s являются векторами с каждым элементом, соответствующим одному уравнению.

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

В начальный момент времени t = t0 для всех x компоненты раствора удовлетворяют начальным условиям формы

$u (x,_{} t0)_{} = u0$ (x).

В MATLAB можно кодировать начальные условия с помощью функции формы.

function u0 = icfun(x)
u0 = 1;
end

В этом случае u0 = 1 определяет начальное условие u0 (x, t0) = 1. Если существует несколько уравнений, тоu0 - вектор с каждым элементом, определяющим начальное условие одного уравнения.

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

На границе x = a или x = b для всех t компоненты решения удовлетворяют граничным условиям формы

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

q (x, t) - диагональная матрица с элементами, которые либо равны нулю, либо никогда не равны нулю. Заметим, что граничные условия выражены в терминах потока f, а не в частных производных u относительно x. Кроме того, из двух коэффициентов p (x, t, u) и q (x, t) только p может зависеть от u.

В MATLAB граничные условия можно кодировать функцией формы.

function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
pL = uL;
qL = 0;
pR = uR - 1;
qR = 0;
end

pL и qL являются коэффициентами для левой границы, в то время как pR и qR - коэффициенты для правой границы. В этом случае bcfun определяет граничные условия

$\begin{array}{l} _{uL} (_{} xL, t) \\ _{=}_{} 0uR (xR \end{array}$ , t) = 1

Если существует несколько уравнений, то выходы pL, qL, pR, и qR - векторы с каждым элементом, определяющим граничное условие одного уравнения.

Варианты интеграции

Свойства интеграции по умолчанию в решателе MATLAB PDE выбираются для обработки распространенных проблем. В некоторых случаях можно повысить производительность решателя, переопределяя эти значения по умолчанию. Для этого используйте odeset для создания options структура. Затем передайте структуру в pdepe в качестве последнего входного аргумента:

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

Из опций базового решателя ОДУ ode15s, доступны только те, которые показаны в следующей таблице pdepe.

Категория	Имя опции
Управление ошибками	`RelTol`, `AbsTol`, `NormControl`
Размер шага	`InitialStep`, `MaxStep`
Ведение журнала событий	`Events`

Оценка решения

После решения уравнения с помощью pdepeMATLAB возвращает решение в виде массива 3-D sol, где sol(i,j,k) содержит k-й компонент раствора, оцененный при t(i) и x(j). В общем, можно извлечь kКомпонент решения th с командой u = sol(:,:,k).

Указанная временная сетка используется исключительно для вывода и не влияет на внутренние временные шаги, выполняемые решателем. Однако указанная пространственная сетка может влиять на качество и скорость решения. После решения уравнения можно использовать pdeval для оценки структуры решения, возвращенной pdepe с другой пространственной сеткой.

Пример: Тепловое уравнение

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

Примером параболического PDE является тепловое уравнение в одном измерении:

$\frac{}{} \frac{^{}}{^{∂u∂t=∂2u∂x2}} .$

Это уравнение описывает рассеяние тепла для $0≤x≤L$ и $t≥0$ . Целью является решение для температуры $u (x,$ t). Температура изначально ненулевая константа, поэтому начальное условие -

$u (x, 0)_{} =$ T0.

Кроме того, температура равна нулю на левой границе и ненулевой на правой границе, поэтому граничные условия

$u (0, t)$ = 0,

$u (L, t)$ = 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).

В этой форме уравнение тепла

$\frac{}{}^{} \frac{}{1⋅∂u∂t=x0∂∂x} (^{} \frac{}{} x0\partialu\partialx) +$ 0.

Значения коэффициентов следующие:

$m =$ 0
$c =$ 1
$\frac{}{f=∂u∂x}$
$s =$ 0

Значение $m$ передается в качестве аргумента pdepe, в то время как другие коэффициенты кодируются в функции для уравнения, которое

function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end

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

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

Функция начального условия для уравнения тепла назначает постоянное значение для $_{u0}$ . Эта функция должна принимать входные данные для $x$ , даже если они не используются.

function u0 = heatic(x)
u0 = 0.5;
end

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

Стандартная форма для граничных условий, ожидаемых pdepe решатель -

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

Записанные в этой форме граничные условия для этой задачи

$u (0, t) + (0\cdotf$ ) = 0,

$(u (L, t) - 1) + ($ 0⋅f) = 0.

Таким образом, значения для $p$ и $q$ являются

$_{pL} =_{} {uL, qL = 0.}_{}$
$_{pR} =_{} {uR-1, qR}_{=} 0$ .

Затем выполняется соответствующая функция.

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

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

Используйте пространственную сетку из 20 точек и временную сетку из 30 точек. Поскольку решение быстро достигает устойчивого состояния, моменты времени, близкие к $t =$ 0, расположены более близко друг к другу, чтобы зафиксировать это поведение на выходе.

L = 1;
x = linspace(0,L,20);
t = [linspace(0,0.05,20), linspace(0.5,5,10)];

Уравнение решения

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

m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);

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

Использовать imagesc для визуализации матрицы решения.

colormap hot
imagesc(x,t,sol)
colorbar
xlabel('Distance x','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')

$Figure contains an axes. The axes with title Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$ contains an object of type image.$

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

function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
function u0 = heatic(x)
u0 = 0.5;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end

Примеры и файлы PDE

Несколько доступных файлов примеров служат отличными отправными точками для большинства распространенных проблем 1-D PDE. Чтобы просмотреть и запустить примеры, используйте приложение Примеры дифференциальных уравнений. Чтобы запустить это приложение, введите

odeexamples

Чтобы открыть отдельный файл для редактирования, введите

edit exampleFileName.m

Для запуска примера введите

exampleFileName

Эта таблица содержит список доступных файлов примеров PDE.

Пример файла	Описание	Пример ссылки
`pdex1`	Простой PDE, иллюстрирующий формулировку, вычисление и построение графика решения.	Решение одного PDE
`pdex2`	Проблема, связанная с разрывами.	Решение PDE с прерыванием
`pdex3`	Проблема, требующая вычисления значений частных производных.	Решение PDE и вычисление частных производных
`pdex4`	Система из двух PDE, решение которых имеет граничные слои на обоих концах интервала и быстро изменяется для малых t.	Система решения ПДЭ
`pdex5`	Система PDE с пошаговыми функциями в качестве исходных условий.	Система решения PDE с пошаговыми функциями начального состояния

Ссылки

[1] Скил, Р. Д. и М. Берзиньш, «Метод пространственной дискретизации параболических уравнений в одной переменной пространства», SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp. 1-32.

См. также

bvp4c | ode45 | odeset | pdepe | pdeval

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