pdepe

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

Синтаксис

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

sol = pdepe (m, pdefun, icfun, bcfun, xmesh, tspan, опции)

[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 функция возвращает значения решения на сетке, предоставленной в 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). xmesh вход содержит пространственную сетку и umesh - решение в точках сетки.

Примеры

свернуть все

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

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

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

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

$\frac{}{} \frac{}{} \frac{}{∂u∂t=1x∂∂x} (\frac{}{} x\partialu\partialx$ ).

Уравнение определяется для $0≤x≤1$ в моменты времени $t≥0$ . Начальное условие определяется в терминах бессельной функции $_{J0} (x$ ) и её первого нуля $n =$ 2.404825557695773 как

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

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

$u (1, t)_{} = {J0}^{(^{n})}$ e-n2t.

Исходные и граничные условия выбираются в соответствии с аналитическим решением задачи.

$u (x, t)_{} = J0^{(^{nx})}$ e-n2t.

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

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

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

Написанный в этой форме PDE становится

$\frac{}{}^{} \frac{}{∂u∂t=x-1∂∂x} (^{} \frac{}{} x1\partialu\partialx$ )

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

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

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

x является независимой пространственной переменной.
t - независимая переменная времени.
u является зависимой переменной, дифференцируемой по отношению к x и t.
dudx - частная пространственная производная $∂u/∂x$ .
Продукция c, f, и s соответствуют коэффициентам в стандартной форме уравнения PDE, ожидаемым 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)_{} = J0$ (nx)

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

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

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

$u (1, t)_{} = {J0}^{(^{n})}$ e-n2t.

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

Для проблем, возникающих в интервале $a≤x≤b$ , граничные условия применяются для всех $t$ и $либо x =$ a, $либо x$ = b. Стандартная форма для граничных условий, ожидаемых решателем:

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

Записанные в этой форме граничные условия для частных производных $u$ должны быть выражены в терминах потока $f (\frac{}{} x,t,u,\partialu\partialx$ ). Таким образом, правое граничное условие для этой задачи

$u (1, t)_{} -J0^{(n^{)}} e-n2t+0\cdotf (\frac{}{}$ x,t,u,∂u∂x) = 0.

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

Исходные данные xl и ul соответствуют $u$ и $x$ для левой границы.
Исходные данные xr и ur соответствуют $u$ и $x$ для правой границы.
t - независимая переменная времени.
Продукция pl и ql соответствуют $_{pL} (x, t,$ u) $_{и} qL ($ x, t) для левой $границы$ (x = 0 для этой задачи).
Продукция pr и qr соответствуют $_{pR} (x, t,$ u) $_{и} qR ($ 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

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

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

Для этой задачи используйте сетку с 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 возвращает решение в массиве 3-D sol, где sol(i,j,k) аппроксимирует kтретий компонент раствора $_{uk}$ оценен на 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. 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
%----------------------------------------------

Решение колебательного PDE с протоколированием событий

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

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

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

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

Уравнение определяется для $0≤x≤1$ и $t≥0.1$ . Начальное условие:

$u (x, 0,1)$ = 1.

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

$u (0, t)$ = 1,

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

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

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

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

$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 уже имеет следующую форму:

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

Таким образом, вы можете прочитать соответствующие термины:

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

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

x является независимой пространственной переменной.
t - независимая переменная времени.
u является зависимой переменной, дифференцируемой по отношению к x и t.
dudx - частная пространственная производная $∂u/∂x$ .
Продукция c, f, и s соответствуют коэффициентам в стандартной форме уравнения PDE, ожидаемым 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) \frac{f}{(}$ x,t,u,∂u∂x) = 0.

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

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

$u (1, t) - cos (δ t) +0\cdotf (\frac{}{}$ x,t,u,∂u∂x) = 0.

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

Исходные данные xl и ul соответствуют $x$ и $u$ для левой границы.
Исходные данные xr и ur соответствуют $x$ и $u$ для правой границы.
t - независимая переменная времени.
Продукция pl и ql соответствуют $_{pL} (x, t,$ u) $_{и} qL ($ x, t) для левой $границы$ (x = 0 для этой задачи).
Продукция pr и qr соответствуют $_{pR} (x, t,$ u) $_{и} qR ($ 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 - пространственная сетка.
umesh содержит решение в точках сетки.
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

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

Перед решением уравнения необходимо указать нужные точки сетки $(t, x).$ pdepe для оценки решения. Для решения этой задачи используйте мелкую сетку из 50 точек в интервалах $0≤x≤1$ и $0.1≤t≤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. 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 решателя PDE для вычисления решения. Кроме того, эти функции можно сохранить в виде собственных файлов в каталоге по пути 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 декартовы координаты без симметрии	$\frac{^{}}{^{Δf=∂2f∂x2}}$
`1`	2-D Цилиндрические координаты с азимутальной угловой симметрией	$\frac{}{} \frac{}{Δf=1ρ∂∂ρ} (\frac{}{} ρ\partialf\partialρ) \frac{}{^{}} \frac{^{}}{^{}}$ +1ρ2∂2f∂φ2, где $\frac{}{} ∂f∂φ=0$ (азимутальная симметрия).
`2`	3-D Сферические координаты с азимутальной и зенитной угловой симметрией	$\frac{}{^{}} \frac{}{Δf=1r2∂∂r} (^{} \frac{}{} r2\partialf\partialr) \frac{}{^{}} \frac{}{} +1r2sinθ∂∂θ (\frac{}{} \frac{sinθ∂f∂θ}{)^{}^{}} \frac{^{}}{^{}}$ +1r2sin2θ∂2f∂φ2, где $\frac{}{} \frac{}{} ∂f∂θ=∂f∂φ=0$ (азимутальные и зенитные угловые симметрии).

0

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

\frac{^{}}{^{Δf=∂2f∂x2}}

1

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

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

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

2

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

$\frac{}{^{}} \frac{}{Δf=1r2∂∂r} (^{} \frac{}{} r2\partialf\partialr) \frac{}{^{}} \frac{}{} +1r2sinθ∂∂θ (\frac{}{} \frac{sinθ∂f∂θ}{)^{}^{}} \frac{^{}}{^{}}$ +1r2sin2θ∂2f∂φ2,

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

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

`pdefun` - функция PDE для уравнений
дескриптор функции

Функция PDE для уравнений, заданная как дескриптор функции, определяющий коэффициенты PDE.

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

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

Члены в уравнении:

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

PDE удерживаются для t0 ≤ t ≤ _tf и a ≤ x ≤ b. Ценности tspan(1) и tspan(end) соответствуют t0 и _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 сохраняет диагональные элементы матрицы с.

Пример

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

$1 · \frac{}{}^{} \frac{}{\partialu\partialt=x0\partial\partialx} (^{} \frac{}{x0∂u∂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 = t0 и всех x компоненты раствора удовлетворяют начальным условиям формы

$u (x,_{} t0)_{} = u0$ (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) \frac{f}{(}$ x,t,u,∂u∂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≤x≤1. В форме, ожидаемой решателем, граничными условиями являются

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

В этой форме ясно, что оба 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` - Пространственная сетка
вектор

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

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

В случае проблем с неоднородностями следует размещать точки сетки на разрывах, чтобы проблема была гладкой в пределах каждого субинтервала. Однако при 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` - Массив решений
3-D массив

Массив решения, возвращаемый как массив 3-D.

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) и точек сетки 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) и точки сетки xmesh, то pdeval вычисляет аппроксимацию и ее частную производную _∂ui/∂x в массиве точек xout и возвращает их в uout и duoutdx: [uout,duoutdx] = pdeval(m,xmesh,uji,xout). pdeval функция вычисляет частную производную _∂ui/∂x а не поток. Поток является непрерывным, но на поверхности раздела материала частная производная может иметь скачок.

Алгоритмы

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

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

Ссылки

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

См. также

ode15s | odeget | odeset | pdeval

Темы

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

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

pdepe

Синтаксис

Описание

Примеры

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

Решение колебательного PDE с протоколированием событий

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

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

`pdefun` - функция PDE для уравнений
дескриптор функции

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

Пример

`icfun` - Функция начального состояния
дескриптор функции

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

Пример

`bcfun` - Функция граничного условия
дескриптор функции

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

Пример

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

`tspan` - Временной интервал интеграции
вектор

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

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

`sol` - Массив решений
3-D массив

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

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

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

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

Совет

Алгоритмы

Ссылки

См. также

Темы

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

Поддержка

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

pdepe

Синтаксис

Описание

Примеры

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

Решение колебательного PDE с протоколированием событий

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

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

pdefun - функция PDE для уравнений дескриптор функции

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

Пример

icfun - Функция начального состояния дескриптор функции

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

Пример

bcfun - Функция граничного условия дескриптор функции

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

Пример

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

tspan - Временной интервал интеграции вектор

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

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

sol - Массив решений 3-D массив

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

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

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

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

Совет

Алгоритмы

Ссылки

См. также

Темы

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

Поддержка

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

`pdefun` - функция PDE для уравнений
дескриптор функции

`icfun` - Функция начального состояния
дескриптор функции

`bcfun` - Функция граничного условия
дескриптор функции

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

`tspan` - Временной интервал интеграции
вектор

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

`sol` - Массив решений
3-D массив

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

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

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

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