pdepe

Решите начальные краевые задачи для параболически-эллиптических УЧП в 1D

Синтаксис

sol = pdepe (m, pdefun, icfun, bcfun, xmesh, tspan)
sol = pdepe (m, pdefun, icfun, bcfun, xmesh, tspan, опции)
[sol, tsol, единственный, te, т.е.] = pdepe (m, pdefun, icfun, bcfun, xmesh, tspan, опции)

Аргументы

m

Параметр, соответствующий симметрии проблемы. m может быть плитой = 0, цилиндрический = 1, или сферический = 2.

pdefun

Указатель на функцию, которая задает компоненты УЧП.

icfun

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

bcfun

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

xmesh

Вектор [x0, x1..., xn] определение точек, в которых числовое решение требуют на каждое значение в tspan. Элементы xmesh должны удовлетворить       x0 < x1 < ... < xn. Длина xmesh должна быть >  = 3.

tspan

Вектор [t0, t1..., tf] определение точек, в которых решение требуют на каждое значение в xmesh. Элементы tspan должны удовлетворить       t0 < t1 < ... < tf. Длиной tspan должен быть >= 3.

опции

Некоторые опции базового решателя ОДУ доступны в pdepe: RelTol, AbsTol, NormControl, InitialStep, MaxStep и Events. В большинстве случаев значения по умолчанию для этих опций предоставляют удовлетворительные решения. Смотрите odeset для деталей.

Описание

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

Параметризация Функций объясняет, как предоставить дополнительные параметры функциям pdefun, icfun или bcfun, при необходимости.

pdepe решает УЧП формы:

c (x, t, u, ∂u∂x) ∂u∂t=x−m ∂∂ x (xmf (x, t, u, ∂u∂x)) +s (x, t, u, ∂u∂x)(1)

УЧП содержат для t0 ≤ t ≤ tf и xb. Интервал [a, b] должен быть конечным. m может быть 0, 1, или 2, соответствуя плите, цилиндрической, или сферической симметрии, соответственно. Если m> 0, то - ≥ 0.

В Уравнении 1, f (x, t, u, ∂u / ∂ x) термин потока, и s (x, t, u, ∂u / ∂ x) является характеристиками выброса. Связь частных производных относительно времени ограничивается умножением диагональной матрицей c (x, t, u, ∂u / ∂ x). Диагональные элементы этой матрицы являются или тождественно нулем или положительный. Элемент, который является тождественно нулевым, соответствует эллиптическому уравнению и в противном случае к параболическому уравнению. Должно быть по крайней мере одно параболическое уравнение. Элемент c, который соответствует параболическому уравнению, может исчезнуть в изолированных значениях x, если те значения x являются точками mesh. Разрывы в c и/или s из-за материальных интерфейсов разрешены при условии, что точка mesh помещается в каждый интерфейс.

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

u (x, t0) =u0 (x)(2)

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

p (x, t, u) +q (x, t) f (x, t, u, ∂u∂x) =0(3)

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

В вызове sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan):

  • m соответствует m.

  • xmesh(1) и xmesh(end) соответствуют a и b.

  • tspan(1) и tspan(end) соответствуют t0 и tf.

  • pdefun вычисляет условия c, f, и s (Уравнение 1). Это имеет форму

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

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

  • icfun оценивает начальные условия. Это имеет форму

    u = icfun(x)
    

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

  • bcfun оценивает условия p и q граничных условий (Уравнение 3). Это имеет форму

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

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

pdepe возвращает решение как многомерный массив sol. ui = ui = sol (:, :, i) является приближением к i th компонент вектора решения 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, чтобы вычислить приближение и его частную производную ∂ui / ∂ x в точках, не включенных в xmesh. Смотрите pdeval для деталей.

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) решает как выше с параметрами интегрирования по умолчанию, замененными значениями в options, аргумент, созданный с функцией odeset. Только некоторые опции базового решателя ОДУ доступны в pdepe: RelTol, AbsTol, NormControl, InitialStep и MaxStep. Значения по умолчанию, полученные, бросая входной параметр options, обычно будут удовлетворительными. Смотрите odeset для деталей.

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) со свойством 'Events' в наборе options к указателю на функцию Events, решает как выше, также находя, где функции события g(t,u(x,t)) являются нулем. Для каждой функции вы задаете, должно ли интегрирование отключить в нуле и имеет ли направление нулевого пересечения значение. Три вектора - столбца возвращены events: [value,isterminal,direction] = events(m,t,xmesh,umesh). xmesh содержит пространственную mesh, и umesh является решением в точках mesh. Используйте pdeval, чтобы оценить решение между точками mesh. Для функции события I-th value(i) является значением функции, ISTERMINAL(I) = 1, если интегрирование должно отключить в нуле этой функции события и 0 в противном случае. direction(i) = 0, если все нули должны быть вычислены (значение по умолчанию), +1, если только нули, где функция события увеличивается, и-1, если только нули, где функция события уменьшается. Вывод tsol является вектором - столбцом времен, заданных в tspan до первого терминального события. SOL(j,:,:) является решением в T(j). TE является вектором времен, в которые события имеют место. SOLE(j,:,:) является решением в TE(j), и индексы в векторном IE задают, который имело место событие.

Если UI = SOL(j,:,i) аппроксимирует компонент i из решения во время, TSPAN(j) и mesh указывают XMESH, pdeval оценивает приближение и его частную производную ∂ui / ∂ x в массиве точек XOUT и возвращает их в UOUT и DUOUTDX: [UOUT,DUOUTDX] = PDEVAL(M,XMESH,UI,XOUT)

Примечание

Частная производная ∂ui / ∂ x оценена здесь, а не поток. Поток непрерывен, но в материале взаимодействуют через интерфейс, частная производная может иметь скачок.

Примеры

Пример 1. Этот пример иллюстрирует прямую формулировку, вычисление и графическое изображение решения единственного УЧП.

π2∂u∂t = ∂∂ x (∂u∂x)

Это уравнение содержит через определенный интервал 0 ≤ x ≤ 1 в течение многих времен t ≥ 0.

УЧП удовлетворяет начальное условие

u (x, 0) =sinπx

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

u (0, t) ≡0πe−t + u∂x (1, t) =0

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

function pdex11

m = 0;
x = linspace(0,1,20);
t = linspace(0,2,5);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);

% A surface plot is often a good way to study a solution.
surf(x,t,u) 
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = pi^2;
f = DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = sin(pi*x);
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

В этом примере УЧП, начальное условие и граничные условия закодированы в локальных функциях pdex1pde, pdex1ic и pdex1bc.

Объемная поверхностная диаграмма показывает поведение решения.

Следующий график показывает профиль решений в окончательном значении t (i. e.,   t = 2).

Пример 2. Этот пример иллюстрирует решение системы УЧП. Проблема имеет пограничные слои в обоих концах интервала. Решение изменяется быстро для маленького t.

УЧП

∂u1∂t=0.024∂2u1∂x2−F(u1−u2) ∂u2∂t=0.170∂2u2∂x2+F (u1−u2)

где F (y) = exp (5,73 лет) – exp (–11.46y).

Это уравнение содержит через определенный интервал 0 ≤ x ≤ 1 в течение многих времен t ≥ 0.

УЧП удовлетворяет начальные условия

u1 (x, 0) ≡1u2 (x, 0) ≡0

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

∂u1∂x (0, t) ≡0u2 (0, t) ≡0u1 (1, t) ≡1∂u2∂x (1, t) ≡0

В форме, ожидаемой pdepe, уравнения

[11] * ∂∂ t [u1u2] = ∂∂ x [0.024 (∂u1 / x) 0.170 (∂u2 / x)] + [−F (u1−u2) F (u1−u2)]

Граничные условия на частных производных вас должны быть записаны с точки зрения потока. В форме, ожидаемой pdepe, левое граничное условие

[0u2] + [10] * [0.024 (∂u1 / x) 0.170 (∂u2 / x)] = [00]

и правильное граничное условие

[u1−10] + [01] * [0.024 (∂u1 / x) 0.170 (∂u2 / x)] = [00]

Решение изменяется быстро для маленького t. Программа выбирает размер шага вовремя, чтобы разрешить это резкое изменение, но видеть это поведение в графиках, пример должен выбрать выходные времена соответственно. Существуют пограничные слои в решении в обоих концах [0,1], таким образом, точки mesh мест в качестве примера около 0 и 1, чтобы разрешить эти резкие изменения. Часто некоторое экспериментирование необходимо, чтобы выбрать mesh, которая показывает поведение решения.

function pdex4
m = 0;
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [1; 1]; 
f = [0.024; 0.17] .* DuDx; 
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F]; 
% --------------------------------------------------------------
function u0 = pdex4ic(x);
u0 = [1; 0]; 
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [0; ul(2)]; 
ql = [1; 0]; 
pr = [ur(1)-1; 0]; 
qr = [0; 1]; 

В этом примере УЧП, начальные условия и граничные условия закодированы в локальных функциях pdex4pde, pdex4ic и pdex4bc.

Объемные поверхностные диаграммы показывают поведение компонентов решения.

Советы

  • Массивы xmesh и tspan играют различные роли в pdepe.

    tspan Функция pdepe выполняет интегрирование времени с решателем ОДУ, который выбирает и временной шаг и формулу динамически. Элементы tspan просто задают, где вы хотите ответы, и стоимость зависит слабо от длины tspan.

    xmesh – Приближения второго порядка к решению сделаны на mesh, заданной в xmesh. Обычно лучше использовать близко расположенные точки mesh, где решение изменяется быстро. pdepe не выбирает mesh в x автоматически. Необходимо обеспечить соответствующую фиксированную mesh в xmesh. Стоимость зависит строго от длины xmesh. Когда m> 0, не необходимо использовать мелкую сетку рядом x = 0, чтобы составлять координатную особенность.

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

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

    Никакая корректировка не необходима для элементов вектора начальных условий, которые соответствуют параболическим уравнениям.

Ссылки

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

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

Была ли эта тема полезной?