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

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

Какие типы УЧП можно решить с MATLAB?

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

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

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

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

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

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

Решение 1D УЧП

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

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

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

УЧП содержат для t ₀ ≤ t ≤ t _f и a ≤ x ≤ b.
Пространственный интервал [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, если они - точки mesh (точки, где решение оценено). Разрывы в c и s из-за материальных интерфейсов разрешены при условии, что точка mesh помещается в каждый интерфейс.

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

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

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

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

Уравнения

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

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

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

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

В этом случае u0 = 1 задает начальное условие u ₀ (x, t ₀) = 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 выбраны, чтобы решить типичные проблемы. В некоторых случаях можно улучшать производительность решателя путем переопределения этих значений по умолчанию. Для этого использовать odeset создать options структура. Затем передайте структуру pdepe как последний входной параметр:

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

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

Категория	Имя опции
Контроль ошибок	`RelTol'AbsTol'` `'NormControl'`
Неродной размер	`InitialStep'MaxStep'`
Регистрация событий	`Events`

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

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

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

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

Скрипт Open 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

Выберите Solution Mesh

Используйте пространственную сетку 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 object. The axes object with title Heat Equation for 0 less equals x less equals 1 and 0 less equals t less equals 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

Примеры УЧП и файлы

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

odeexamples

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

edit exampleFileName.m

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

exampleFileName

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

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

Ссылки

[1] Skeel, R. D. и М. Берзиньш, "Метод для Пространственной Дискретизации Параболических уравнений в Одной Пространственной переменной", SIAM Journal на Научном и Статистическом Вычислении, Издании 11, 1990, стр 1–32.

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