pdepe

Решение 1-D параболических и эллиптических PDE

Синтаксис

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

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

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

Описание

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) решает систему параболических и эллиптических PDE с одной пространственной переменной x и временной t. По крайней мере, одно уравнение должно быть параболическим. Скалярная m представляет симметрию задачи (плита, цилиндрическая или сферическая). Решаемые уравнения кодируются в pdefunначальное значение кодируется в icfun, и граничные условия закодированы в bcfun. Обыкновенные дифференциальные уравнения (ОДУ), следующие из дискретизации в пространстве, интегрируются, чтобы получить приблизительные решения в моменты времени, заданные в tspan. pdepe функция возвращает значения решения на сетке mesh, представленной в xmesh.

пример

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) также использует настройки интегрирования, заданные как options, который создается с помощью odeset функция. Для примера используйте AbsTol и RelTol опции для задания абсолютных и относительная погрешность допусков.

пример

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) также находит, где функции (t, u (x, t)), называемые функциями события, равны нулю. В выходах te - время события, sole является решением во время события и ie - индекс инициируемого события. tsol - вектор-столбец времени, заданный в tspan, до первого терминального события.

Для каждой функции события задайте, должно ли интегрирование завершаться в нуле и имеет ли значение направление пересечения нулем. Сделайте это, установив 'Events' опция odeset к функции, такой как @myEventFcn, и создание соответствующей функции: [value, isterminal, direction] = myEventFcn(m, t, xmesh, umesh). The xmesh вход содержит пространственный mesh и umesh - решение в точках mesh.

Примеры

свернуть все

Решение уравнения тепла в цилиндрических координатах

Открыть Live Script

Решить тепловое уравнение в цилиндрических координатах можно используя pdepe, и постройте график решения.

В цилиндрических координатах с угловой симметрией тепловое уравнение

$\frac{\partial u}{\partial t} = \frac{1}{x} \frac{\partial}{\partial x} (x \frac{\partial u}{\partial x}) .$

Уравнение задано для $0 \leq x \leq 1$ время от времени $t \geq 0$ . Начальное условие определяется терминами функции бесселя $J_{0} (x)$ и его первый нуль $n = 2.404825557695773$ как

$u (x, 0) = J_{0} (nx) .$

Поскольку эта задача находится в цилиндрических координатах (m = 1), pdepe автоматически применяет условие симметрии в $x = 0$ . Правое граничное условие

$u (1, t) = J_{0} (n) e^{- n^{2} t} .$

Начальные и граничные условия выбраны так, чтобы соответствовать аналитическому решению задачи,

$u (x, t) = J_{0} (nx) e^{- n^{2} t} .$

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

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

Прежде чем вы сможете кодировать уравнение, необходимо переписать его в форме, которую pdepe решатель ожидает. Стандартная форма, которая 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}) .$

Написанный в этой форме, УЧП становится

$\frac{\partial u}{\partial t} = x^{- 1} \frac{\partial}{\partial x} (x^{1} \frac{\partial u}{\partial x})$

С помощью уравнения в надлежащей форме можно считать релевантные условия:

$m = 1$
$c (x, t, u, \frac{\partial u}{\partial x}) = 1$
$f (x, t, u, \frac{\partial u}{\partial x}) = \frac{\partial u}{\partial x}$
$s (x, t, u, \frac{\partial u}{\partial x}) = 0$

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

x - независимая пространственная переменная.
t - независимая временная переменная.
u - зависимая переменная, дифференцируемая по отношению к x и t.
dudx - частная пространственная производная $\partial u / \partial x$ .
Выходные выходы c, f, и s соответствуют коэффициентам в стандартной форме уравнения УЧП, ожидаемым по pdepe. Эти коэффициенты кодируются в терминах входных переменных x, t, u, и dudx.

В результате уравнение в этом примере может быть представлено функцией:

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

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

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

Затем напишите функцию, которая возвращает начальное условие. Начальное условие применяется при первом временном значении tspan(1). Функция должна иметь подпись u0 = heatic(x).

Соответствующая функция для $u (x, 0) = J_{0} (nx)$ является

function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end

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

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

$u (1, t) = J_{0} (n) e^{- n^{2} t} .$

Поскольку эта задача находится в цилиндрических координатах (m = 1), pdepe автоматически применяет условие симметрии в $x = 0$ , поэтому вам не нужно задавать левое граничное условие.

Для задач, поставленных на интервале $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})$ . Таким образом, правильное граничное условие для этой задачи является

$u (1, t) - J_{0} (n) e^{- n^{2} t} + 0 \cdot f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

Граничная функция должна иметь сигнатуру функции [pl,ql,pr,qr] = heatbc(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] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end

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

Перед решением уравнения необходимо задать точки mesh $(t, x)$ при котором вы хотите pdepe для оценки решения. Задайте точки как векторы t и x. Векторы t и x играть различные роли в решателе. В частности, стоимость и точность решения сильно зависят от длины вектора x. Однако расчет намного менее чувствительно к значениям в векторе t.

Для этой задачи используйте mesh с 25 равномерно разнесенными точками в пространственном интервале [0,1] для обоих x и t.

x = linspace(0,1,25);
t = linspace(0,1,25);

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

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

m = 1;
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

pdepe возвращает решение в трехмерный массив sol, где sol(i,j,k) аппроксимирует kпервый компонент решения $u_{k}$ оценивается в t(i) и x(j). Размер sol является length(t)-by- length(x)-by- length(u0), с u0 задает начальное условие для каждого компонента решения. Для этой задачи u имеет только один компонент, так что sol является матрицей 25 на 25, но в целом можно извлечь kкомпонент решения с командой u = sol(:,:,k).

Извлеките первый компонент решения из sol.

u = sol(:,:,1);

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

Создать объемную поверхностную диаграмму решения. Поскольку задача поставлена с цилиндрическими координатами на диске, x значения показывают температуру на диске на некотором расстоянии от центра и t значения показывают, как температура в конкретном месте изменяется с течением времени.

surf(x,t,u)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])

Figure contains an axes. The axes contains an object of type surface.

Постройте график изменения температуры в центре ( $x = 0$ ) диска.

plot(t,sol(:,1))
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at center of disc')

Figure contains an axes. The axes with title Temperature change at center of disc contains an object of type line.

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

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

function [c,f,s] = heatcyl(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
%----------------------------------------------
function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end
%----------------------------------------------
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end
%----------------------------------------------

Решение колебательного УЧП с логгированием событий

Открыть Live Script

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

Рассмотрим уравнение

$\frac{1}{x} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} (\frac{1}{t} u) .$

Уравнение задано для $0 \leq x \leq 1$ и $t \geq 0.1$ . Начальное условие является

$u (x, 0.1) = 1 .$

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

$u (0, t) = 1,$

$u (1, t) = \cos (π t) .$

Кроме того, представляют интерес пересечения решения с пересечениями нулем.

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

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

Уравнение УЧП уже в таком виде:

$\frac{1}{x} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} (\frac{1}{t} u) .$

Так можно прочитать соответствующие условия:

$m = 0$
$c (x, t, u, \frac{\partial u}{\partial x}) = \frac{1}{x}$
$f (x, t, u, \frac{\partial u}{\partial x}) = \frac{1}{t} u$
$s (x, t, u, \frac{\partial u}{\partial x}) = 0$

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

x - независимая пространственная переменная.
t - независимая временная переменная.
u - зависимая переменная, дифференцируемая по отношению к x и t.
dudx - частная пространственная производная $\partial u / \partial x$ .
Выходные выходы c, f, и s соответствуют коэффициентам в стандартной форме уравнения УЧП, ожидаемым по pdepe. Эти коэффициенты кодируются в терминах входных переменных x, t, u, и dudx.

В результате уравнение в этом примере может быть представлено функцией:

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

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

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

Затем напишите функцию, которая возвращает начальное условие. Начальное условие применяется при первом временном значении tspan(1). Функция должна иметь подпись u0 = oscic(x).

Соответствующая функция для $u (x, 0.1) = 1$ является

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

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

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

$u (0, t) = 1,$

$u (1, t) = \cos (π t) .$

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

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

$u (0, t) - 1 + 0 \cdot f (x, t, u, \frac{\partial u}{\partial x}) = 0,$

$u (1, t) - \cos (π t) + 0 \cdot f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

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

Входы xl и ul соответствуют $x$ и $u$ для левого контура.
Входы xr и ur соответствуют $x$ и $u$ для правого контура.
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] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end

Функция события кода

Используйте функцию события для регистрации пересечений нулем решения в интегрировании. Функция события имеет сигнатуру функции [value,isterminal,direction] = pdevents(m,t,xmesh,umesh):

m - симметрия координат, заданная в качестве первого входа pdepe.
t - текущее время (скаляр).
xmesh - пространственный mesh.
umesh содержит решение в точках mesh.
value является интересующим уравнением, обычно выраженным в терминах решения umesh. Когда value равен 0, происходит событие.
isterminal определяет, останавливает ли событие интегрирование. Если isterminal 0, события регистрируются, но интегрирование не останавливается. Если isterminal равен 1, тогда интегрирование останавливается, когда происходит событие.
direction задает направление пересечения нуля. Если 1, событие инициируется только пересечениями нулем с положительным наклоном. Если -1, пересечения нулем должны иметь отрицательный уклон. Если 0, то любое пересечение нулем запускает событие.

В каждый момент интегрирования решатель вызывает функцию события для проверки пересечений нуля. Чтобы логгировать все пересечения нуля, value должен искать изменения в знаке в векторе решения umesh. Задайте isterminal и direction как нулевые векторы с таким же размером, поскольку события в этом примере не являются терминальными, и пересечения нулем могут происходить с любым наклоном.

Функция события для этой задачи

function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end

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

Перед решением уравнения необходимо задать точки mesh $(t, x)$ при котором вы хотите pdepe для оценки решения. Для этой задачи используйте мелкий mesh 50 точек в интервалах $0 \leq x \leq 1$ и $0.1 \leq t \leq 1$ . Мелкий mesh даёт хорошее разрешение колебательного решения.

x = linspace(0,1,50);
t = linspace(0.1,pi,50);

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

Наконец, решить уравнение используя симметрию m, уравнение УЧП, начальные условия, граничные условия, функция события и сетки для x и t. Использование odeset чтобы создать структуру опций, которая ссылается на функцию событий и передается в структуре как последний входной параметр pdepe. Задайте пять выходных аргументов, чтобы вернуть информацию как от функции события, так и от решателя:

sol - решение вычисляется pdepe.
tsol является вектором времени перед терминальным событием. tsol равно t когда никакие события не являются терминалами.
sole является решением во время каждого события.
te - время каждого события.
ie - индекс каждого события. Начиная с values = umesh в функции события, ie приводит индексы umesh который инициировал событие на каждом временном шаге.

m = 0;
options = odeset('Events',@pdevents);
[sol,tsol,sole,te,ie] = pdepe(m,@oscpde,@oscic,@oscbc,x,t,options);

Извлеките решение как матрицу u.

u = sol(:,:,1);

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

Создать объемную поверхностную диаграмму решения и просмотреть график сверху.

surf(x,t,u)
view(2)

Figure contains an axes. The axes contains an object of type surface.

Постройте график точек, где произошли события, с поверхностью $f (x, t) = 0$ для ссылки. Выходной вектор индекса ie полезно выбрать местоположения событий. Выражение x(ie)' задает значения x, где произошли события, и выражение sole(x==x(ie)') задает соответствующие значения решения.

view([39 30])
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
hold on
plot3(x(ie)',te,sole(x==x(ie)'),'r*')
surf(x,t,zeros(size(u)),'EdgeColor','flat')
hold off

Figure contains an axes. The axes contains 3 objects of type surface, line.

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

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

function [c,f,s] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end
%----------------------------------------------
function u0 = oscic(x)
u0 = 1;
end
%----------------------------------------------
function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end
%----------------------------------------------
function [value, isterminal, direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
%----------------------------------------------

Входные параметры

свернуть все

`m` - Константа симметрии
`0` | `1` | `2`

Константа симметрии, заданная как одно из значений в этой таблице. m задает тип задачи в уравнениях, заданных pdefun. Как только вы переписываете уравнения в форме, ожидаемой pdepe, можно считать значение m из уравнения.

Значение Описание Лаплакиан потока

Значение	Описание	Лаплакиан потока
`0`	1-D Декартовы координаты без симметрии	$Δ f = \frac{\partial^{2} f}{\partial x^{2}}$
`1`	2-D цилиндрические координаты с азимутальной угловой симметрией	$Δ f = \frac{1}{ρ} \frac{\partial}{\partial ρ} (ρ \frac{\partial f}{\partial ρ}) + \frac{1}{ρ^{2}} \frac{\partial^{2} f}{\partial φ^{2}}$ , где $\frac{\partial f}{\partial φ} = 0$ (азимутальная симметрия).
`2`	3-D Сферические координаты с азимутальной и зенитной угловой симметриями	$Δ f = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial f}{\partial r}) + \frac{1}{r^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial f}{\partial θ}) + \frac{1}{r^{2} \sin^{2} θ} \frac{\partial^{2} f}{\partial φ^{2}}$ , где $\frac{\partial f}{\partial θ} = \frac{\partial f}{\partial φ} = 0$ (азимутальные и зенитные угловые симметрии).

0

1-D Декартовы координаты без симметрии

Δ f = \frac{\partial^{2} f}{\partial x^{2}}

1

2-D цилиндрические координаты с азимутальной угловой симметрией

$Δ f = \frac{1}{ρ} \frac{\partial}{\partial ρ} (ρ \frac{\partial f}{\partial ρ}) + \frac{1}{ρ^{2}} \frac{\partial^{2} f}{\partial φ^{2}}$ ,

где $\frac{\partial f}{\partial φ} = 0$ (азимутальная симметрия).

2

3-D Сферические координаты с азимутальной и зенитной угловой симметриями

$Δ f = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial f}{\partial r}) + \frac{1}{r^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial f}{\partial θ}) + \frac{1}{r^{2} \sin^{2} θ} \frac{\partial^{2} f}{\partial φ^{2}}$ ,

где $\frac{\partial f}{\partial θ} = \frac{\partial f}{\partial φ} = 0$ (азимутальные и зенитные угловые симметрии).

Когда m > 0, pdepe требует, чтобы левое пространственное граничное a было ≥ 0, и это автоматически накладывает левое граничное условие для расчета разрыва координат в источник. В этом случае, pdepe игнорирует все заданные левые граничные условия в bcfun.

`pdefun` - функция PDE для уравнений
указатель на функцию

Функция PDE для уравнений, заданная как указатель на функцию, который задает коэффициенты УЧП.

pdefun кодирует коэффициенты для PDE вида

Условия в уравнении:

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

УЧП держится для <reservedrangesplaceholder6> 0 <reservedrangesplaceholder5> <reservedrangesplaceholder4> и <reservedrangesplaceholder3> ≤ <reservedrangesplaceholder2> ≤ <reservedrangesplaceholder1>. Значения tspan(1) и tspan(end) соответствуют t 0 и _tf, в то время как xmesh(1) и xmesh(end) соответствуют a и b. Интервал [a, b] должен быть конечным. m может быть 0, 1 или 2, соответствующий плите, цилиндрической или сферической симметрии, соответственно. Если m > 0, то a должно быть ≥ 0.

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

Напишите функцию pdefun который вычисляет значения коэффициентов c, f и s. Используйте сигнатуру функции

[c,f,s] = pdefun(x,t,u,dudx)

Входные параметры в pdefun являются скалярами x и t и векторы u и dudx. Векторная u аппроксимирует u решения и dudx аппроксимирует его частную производную относительно x. Выходные аргументы c, f, и s Векторы-столбцы с количеством элементов, равным количеству уравнений. c сохраняет диагональные элементы матричной c.

Пример

Тепловое уравнение $\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}}$ может быть переписан в форме, ожидаемой решателем как

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

В этой форме Вы можете прочитать значение коэффициентов как m = 0, c = 1, f =  <reservedrangesplaceholder2> / <reservedrangesplaceholder1> и s = 0. Функция для этого уравнения:

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

Типы данных: function_handle

`icfun` - Функция начальных условий
указатель на функцию

Функция начального условия, заданная как указатель на функцию, который задает начальное условие.

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

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

Кодирование начального условия

Напишите функцию icfun который определяет начальные условия. Используйте сигнатуру функции:

function u0 = icfun(x)

При вызове с аргументом x, функция icfun оценивает и возвращает начальные значения компонентов решения в x в векторе-столбце u0. Количество элементов в u0 равен количеству уравнений.

Пример

Постоянное начальное условие $u (x, 0) = 1 / 2$ соответствует этой функции:

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

Типы данных: function_handle

`bcfun` - Функция граничных условий
указатель на функцию

Функция граничного условия, заданная как указатель на функцию, который задает граничные условия.

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

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

Элементы q либо равны нулю, либо никогда не равны нулю. Обратите внимание, что граничные условия выражены с точки зрения потока f, а не  <reservedrangesplaceholder3> /  <reservedrangesplaceholder2>. Кроме того, из двух коэффициентов только p могут зависеть от u.

Кодирование граничных условий

Напишите функцию bcfun который определяет условия p и q границ. Используйте сигнатуру функции

function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

ul - приблизительное решение на левом контуре xl = a, и ur является приблизительным решением на правом контуре xr = b.
pl и ql Векторы-столбцы, соответствующие p и q, оцениваются в xl. Точно так же pr и qr соответствуют xr.
Количество элементов в выходах pl, ql, pr, и qr равен количеству уравнений.

Когда m > 0 и a = 0, граница решения около x = 0 требует, чтобы поток f исчезнуть при a = 0. pdepe автоматически накладывает это граничное условие и игнорирует значения, возвращенные в pl и ql.

Пример

Учитывайте граничные условия $u (0, t) = 0, u (1, t) = 1$ на интервале $0 \leq x \leq 1$ . В форме, ожидаемой решателем, граничные условия

$\begin{array}{l} u (0, t) + 0 \cdot f (0, t, u, \frac{\partial u}{\partial x}) = 0, \\ u (1, t) - 1 + 0 \cdot f (0, t, u, \frac{\partial u}{\partial x}) = 0. \end{array}$

В этой форме понятно, что оба q условия равны нулю. Условия p ненулевые, чтобы отражать значения u. Соответствующая функция является

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

Типы данных: function_handle

`xmesh` - Пространственный mesh
вектор

Пространственный mesh, заданная как вектор [x0 x1 ... xn] определение точек, в которых запрашивается числовое решение для каждого значения в tspan. Приближения решения второго порядка выполняются на mesh, указанной в xmesh. pdepe не выбирает mesh в x автоматически. Необходимо предоставить соответствующую фиксированный mesh в xmesh.

Элементы xmesh должен удовлетворять x0 < x1 < ... < xn. Как правило, лучше всего использовать близко расположенные точки mesh, где решение изменяется быстро. Длина xmesh должен быть >=3, и вычислительная стоимость решения сильно зависит от длины xmesh.

Для задач с разрывами следует поместить точки mesh в разрывы так, чтобы задача была гладкой в каждом подынтервале. Однако, когда m > 0, не обязательно использовать мелкий mesh около x = 0, чтобы принять во внимание сингулярность координат.

Пример: xmesh = linspace(0,5,25) использует пространственный mesh 25 точек между 0 и 5.

Типы данных: single | double

`tspan` - Временной промежуток интегрирования
вектор

Временной промежуток интегрирования, заданный как вектор [t0 t1 ... tf] определение точек, в которых запрашивается решение для каждого значения в xmesh. pdepe функция выполняет интегрирование во времени с решателем ОДУ, который выбирает и временной шаг, и формулу динамически. Элементы tspan просто укажите, где вы хотите ответы, и поэтому вычислительная стоимость решения зависит только слабо от длины tspan.

Элементы tspan должен удовлетворять t0 < t1 < ... < tf. Длина tspan должен быть >=3.

Пример: tspan = linspace(0,5,5) использует пять временных точек от 0 до 5.

Типы данных: single | double

`options` - Структура опций
массив структур

Структура опции, заданная как массив структур. Используйте odeset функция для создания или изменения структуры опций. pdepe поддерживает следующие опции:

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

Пример: options = odeset('RelTol',1e-5) задает относительную погрешность допуск 1e-5.

Типы данных: struct

Выходные аргументы

свернуть все

`sol` - Массив решений
трехмерный массив

Массив решения, возвращенный как трехмерный массив.

pdepe возвращает решение как многомерный массив. _ui = ui = sol(:,:, i) является приближением к iпервый компонент векторной u решения. Элемент ui(j, k) = sol(j, k, i) аппроксимирует _ui в (t, x ) = (tspan(j), xmesh(k)).

ui = sol(j,:, i) аппроксимирует компонентные i решения во время tspan(j) и mesh точек xmesh(:). Использовать pdeval вычислить приближение и его _{частную производную <reservedrangesplaceholder2>} / <reservedrangesplaceholder1> в точках, не включенных в xmesh. Посмотрите pdeval для получения дополнительной информации.

`tsol` - Время решения
Вектор-столбец

Время решения, возвращенное как вектор-столбец времени, заданное в tspan которые находятся до первого терминального события.

`sole` - Решение во время событий
массив

Решение во время событий, возвращаемое как массив. Время событий в te соответствуют решениям, возвращаемым в sole, и ie определяет, какое событие произошло.

`te` - Время событий
Вектор-столбец

Время событий, возвращаемое как вектор-столбец. Время событий в te соответствуют решениям, возвращаемым в sole, и ie определяет, какое событие произошло.

`ie` - Индекс функции инициируемого события
Вектор-столбец

Индекс функции инициируемого события, возвращаемый как вектор-столбец. Время событий в te соответствуют решениям, возвращаемым в sole, и ie определяет, какое событие произошло.

Совет

Если uji = sol(j,:,i) аппроксимирует i компонента решения во время tspan(j) и mesh xmesh, затем pdeval оценивает приближение и его частную производную ∂ _ui/ ∂ x в массиве точек xout и возвращает их в uout и duoutdx: [uout,duoutdx] = pdeval(m,xmesh,uji,xout). pdeval функция оценивает частную производную ∂ _ui/ ∂ x а не поток. Поток непрерывен, но при взаимодействии материала частная производная может иметь переход.

Алгоритмы

Интегрирование во времени осуществляется с ode15s решатель. pdepe использует возможности ode15s для решения дифференциально-алгебраических уравнений, возникающих, когда УЧП содержит эллиптические уравнения, и для обращения с якобианцами с заданным шаблоном разреженности.

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

Ссылки

[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, ppp.1-32.

См. также

ode15s | odeget | odeset | pdeval

Темы

Представлено до R2006a

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

pdepe

Синтаксис

Описание

Примеры

Решение уравнения тепла в цилиндрических координатах

Решение колебательного УЧП с логгированием событий

Входные параметры

`m` - Константа симметрии
`0` | `1` | `2`

`pdefun` - функция PDE для уравнений
указатель на функцию

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

Пример

`icfun` - Функция начальных условий
указатель на функцию

Кодирование начального условия

Пример

`bcfun` - Функция граничных условий
указатель на функцию

Кодирование граничных условий

Пример

`xmesh` - Пространственный mesh
вектор

`tspan` - Временной промежуток интегрирования
вектор

`options` - Структура опций
массив структур

Выходные аргументы

`sol` - Массив решений
трехмерный массив

`tsol` - Время решения
Вектор-столбец

`sole` - Решение во время событий
массив

`te` - Время событий
Вектор-столбец

`ie` - Индекс функции инициируемого события
Вектор-столбец

Совет

Алгоритмы

Ссылки

См. также

Темы

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

Поддержка

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

pdepe

Синтаксис

Описание

Примеры

Решение уравнения тепла в цилиндрических координатах

Решение колебательного УЧП с логгированием событий

Входные параметры

m - Константа симметрии 0 | 1 | 2

pdefun - функция PDE для уравнений указатель на функцию

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

Пример

icfun - Функция начальных условий указатель на функцию

Кодирование начального условия

Пример

bcfun - Функция граничных условий указатель на функцию

Кодирование граничных условий

Пример

xmesh - Пространственный mesh вектор

tspan - Временной промежуток интегрирования вектор

options - Структура опций массив структур

Выходные аргументы

sol - Массив решений трехмерный массив

tsol - Время решения Вектор-столбец

sole - Решение во время событий массив

te - Время событий Вектор-столбец

ie - Индекс функции инициируемого события Вектор-столбец

Совет

Алгоритмы

Ссылки

См. также

Темы

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

Поддержка

`m` - Константа симметрии
`0` | `1` | `2`

`pdefun` - функция PDE для уравнений
указатель на функцию

`icfun` - Функция начальных условий
указатель на функцию

`bcfun` - Функция граничных условий
указатель на функцию

`xmesh` - Пространственный mesh
вектор

`tspan` - Временной промежуток интегрирования
вектор

`options` - Структура опций
массив структур

`sol` - Массив решений
трехмерный массив

`tsol` - Время решения
Вектор-столбец

`sole` - Решение во время событий
массив

`te` - Время событий
Вектор-столбец

`ie` - Индекс функции инициируемого события
Вектор-столбец