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

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

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

MATLAB^® УЧП pdepe решает задачи начального-краевого значения для систем PDE в одной пространственной переменной x и времени t. Можно думать о них как об ОДУ одной переменной, которые также изменяются относительно времени.

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

Уравнения с производной по времени являются parabolic. Примером является тепловое уравнение $\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}}$ .
Уравнения без производной по времени являются elliptic. Примером является уравнение Лапласа $\frac{\partial^{2} u}{\partial x^{2}} = 0$ .

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

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

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

Решение 1-D PDE

1-D УЧП включает u функции (x, t), которая зависит от временных t и одной пространственной переменной x. Решатель MATLAB PDE pdepe решает системы 1-D параболических и эллиптических PDE вида

$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}) .$

Уравнение имеет свойства:

УЧП держится для <reservedrangesplaceholder6> 0 <reservedrangesplaceholder5> <reservedrangesplaceholder4> <reservedrangesplaceholder3> и <reservedrangesplaceholder2> ≤ <reservedrangesplaceholder1> ≤ <reservedrangesplaceholder0>.
Пространственный интервал [a, b] должен быть конечным.
m может быть 0, 1 или 2, соответствующий плите, цилиндрической или сферической симметрии, соответственно. Если m > 0, то a ≥ 0 также должен удерживаться.
Коэффициент $f (x, t, u, \frac{\partial u}{\partial x})$ является термином потока и $s (x, t, u, \frac{\partial u}{\partial x})$ - исходный термин.
Член потока должен зависеть от частной производной ∂ u/ ∂ x.

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

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

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

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{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}}

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

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

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

$u (x, t_{0}) = u_{0} (x) .$

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

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

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

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

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

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial 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} u_{L} (x_{L}, t) = 0 \\ u_{R} (x_{R}, t) = 1 \end{array}$

Если существует несколько уравнений, то выходы 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 возвращает решение как трехмерный массив sol, где sol(i,j,k) содержит kпервый компонент решения, оцениваемый в t(i) и x(j). В целом можно извлечь kкомпонент решения с командой u = sol(:,:,k).

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

Пример: Уравнение тепла

Открыть Live Script

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

$\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}} .$

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

$u (x, 0) = T_{0} .$

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

$u (0, t) = 0,$

$u (L, t) = 1 .$

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

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

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

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

$1 \cdot \frac{\partial u}{\partial t} = x^{0} \frac{\partial}{\partial x} (x^{0} \frac{\partial u}{\partial x}) + 0 .$

Таким образом, значения коэффициентов следующие:

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

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

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

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

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

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

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

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

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

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

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

$u (0, t) + (0 \cdot f) = 0,$

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

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

$p_{L} = u_{L}, q_{L} = 0 .$
$p_{R} = u_{R} - 1, q_{R} = 0$ .

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

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

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

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

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

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

Наконец, решите уравнение, используя симметрию $m$ , уравнение УЧП, начальное условие, граничные условия и сетки для $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 УЧП. Чтобы исследовать и запустить примеры, используйте Примеры дифференциальных уравнений, приложение. Чтобы запустить это приложение, введите

odeexamples

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

edit exampleFileName.m

Чтобы запустить пример, введите

exampleFileName

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

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

Ссылки

[1] Skeel, R. D. and M. Berzins, «A Method for the Spatial Discretization of Parabolic Equations in One Space Variable», SIAM Journal on Scientific and Statistical Computing, Vol., 1990, pp. 1-32.

См. также

bvp4c | ode45 | odeset | pdepe | pdeval

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