pdepe

Решите 1D параболические и эллиптические УЧП

Синтаксис

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) решает систему параболических и эллиптических УЧП с одной пространственной переменной 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, и создание соответствующей функции: Значение, isterminal, direction] = myEventFcnMT, xmesh, umesh). xmesh введите содержит пространственную mesh и umesh решение в точках mesh.

Примеры

свернуть все

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

Скрипт Open 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 = 1pdepe автоматически осуществляет условие симметрии в $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$ .
Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения PDE, ожидаемой pdepe. Эти коэффициенты закодированы в терминах входных переменных xTU, и 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 = 1pdepe автоматически осуществляет условие симметрии в $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

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

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

pdepe возвращает решение в трехмерном массиве sol, где sol(i,j,k) аппроксимирует kкомпонент th решения $u_{k}$ оцененный в t(i) и x(j). Размер sol length(t)- length(x)- length(u0), начиная с u0 задает начальное условие для каждого компонента решения. Для этой проблемы, u имеет только один компонент, таким образом, sol 25 25 матрица, но в целом можно извлечь kкомпонент решения th с командой 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 object. The axes object 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 object. The axes object with title Temperature change at center of disc contains an object of type line.

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

Перечисленный здесь локальные функции помощника что решатель УЧП 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
%----------------------------------------------

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

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

Уравнение кода

Уравнение PDE уже находится в этой форме:

$\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$ .
Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения PDE, ожидаемой pdepe. Эти коэффициенты закодированы в терминах входных переменных xTU, и 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

Выберите Solution Mesh

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

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

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

Наконец, решите уравнение с помощью симметрии m, уравнение PDE, начальные условия, граничные условия, функция события и сетки для 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 object. The axes object 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 object. The axes object 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	1D Декартовы координаты без симметрии	$Δ f = \frac{\partial^{2} f}{\partial x^{2}}$
1	2D Цилиндрические координаты с азимутальной угловой симметрией	$Δ 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$ (азимутальный и зенит угловые симметрии).

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

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

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

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

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

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` — Функция УЧП для уравнений
указатель на функцию

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

pdefun кодирует коэффициенты для УЧП формы

Термины в уравнении:

$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 являются точками mesh.
Разрывы в c и s из-за материальных интерфейсов разрешены при условии, что точка mesh помещается в каждый интерфейс.

УЧП содержат для t ₀ ≤ t ≤ _tf и a ≤ x ≤ b. Значения tspan(1) и tspan(end) соответствуйте t ₀ и _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. Выходные аргументы cF, и 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 = ∂u / ∂ x и s = 0. Функция для этого уравнения:

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

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

`icfun` — Начальная функция условия
указатель на функцию

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

Для t = t ₀ и весь 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, а не ∂u / ∂ x. Кроме того, этих двух коэффициентов только 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, не необходимо использовать мелкую сетку около x = 0 с учетом координатной сингулярности.

Пример: xmesh = linspace(0,5,25) использует пространственную сетку 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компонент th вектора решения u. Элемент uiJK) = solJKi) аппроксимирует _ui в (t, x) = (tspanJ), xmeshK)).

ui = solJ,:i) аппроксимирует i компонента из решения во время tspanJ) и mesh указывает xmesh(:)Использование pdeval вычислить приближение и его частную производную ∂_ui / ∂ x в точках, не включенных в 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. и М. Берзиньш, "Метод для Пространственной Дискретизации Параболических уравнений в Одной Пространственной переменной", SIAM Journal на Научном и Статистическом Вычислении, Издании 11, 1990, pp.1-32.

Темы

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

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

pdepe

Синтаксис

Описание

Примеры

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

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

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

`m` — Постоянная симметрия
0 | 1 | 2

`pdefun` — Функция УЧП для уравнений
указатель на функцию

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

Пример

`icfun` — Начальная функция условия
указатель на функцию

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

Пример

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

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

Пример

`xmesh` — Пространственная mesh
вектор

`tspan` — Отрезок времени интегрирования
вектор

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

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

`sol` — Массив решения
Трехмерный массив

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

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

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

`ie` — Индекс инициированной функции события
вектор-столбец

Советы

Алгоритмы

Ссылки

Смотрите также

Темы

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

Поддержка

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

pdepe

Синтаксис

Описание

Примеры

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

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

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

m — Постоянная симметрия0 | 1 | 2

pdefun — Функция УЧП для уравнений указатель на функцию

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

Пример

icfun — Начальная функция условия указатель на функцию

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

Пример

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

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

Пример

xmesh — Пространственная mesh вектор

tspan — Отрезок времени интегрирования вектор

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

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

sol — Массив решения Трехмерный массив

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

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

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

ie — Индекс инициированной функции события вектор-столбец

Советы

Алгоритмы

Ссылки

Смотрите также

Темы

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

Поддержка

`m` — Постоянная симметрия
0 | 1 | 2

`pdefun` — Функция УЧП для уравнений
указатель на функцию

`icfun` — Начальная функция условия
указатель на функцию

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

`xmesh` — Пространственная mesh
вектор

`tspan` — Отрезок времени интегрирования
вектор

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

`sol` — Массив решения
Трехмерный массив

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

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

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

`ie` — Индекс инициированной функции события
вектор-столбец