Решение PDE с прерыванием

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

В этом примере показано, как решить PDE, который взаимодействует с материалом. Интерфейс материала создает разрыв в задаче при $x =$ 0,5, а исходное условие имеет разрыв в правой $границе x$ = 1.

Рассмотрим кусочно PDE

${\begin{matrix} \frac{}{}^{} \frac{}{∂u∂t=x-2∂∂x} (^{} x2 \frac{}{} 5\partialu\partialx) -^{} & 1000eu ( \\ \frac{0≤x≤0.5}{)}^{} \frac{}{}^{} ∂u∂t=x-2∂∂x \frac{(}{x2}^{} & \partialu\partialx) - eu ( \end{matrix}$ 0.5≤x≤1)

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

$\begin{array}{l} u (x, 0) = 0 ( \\ 0≤x<1), u (1,0) = 1 \end{array}$ (x = 1).

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

$\begin{array}{l} \frac{}{} ∂u∂x=0 (x = 0) \\ , u (1, t) = 1 ( \end{array}$ x = 1).

Чтобы решить это уравнение в 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 находится в правильной форме, поэтому можно считывать значения коэффициентов.

Значения для члена потока $f (\frac{}{} x,t,u,\partialu\partialx$ ) и члена источника $s (\frac{}{}$ x.t,u,∂u∂x) изменяются в зависимости от значения x. Коэффициенты:

$m =$ 2

$c (\frac{}{} x,t,u,\partialu\partialx)$ = 1

${\begin{matrix} f (\frac{}{} x,t,u,\partialu\partialx) \frac{}{} & =5∂u∂x ( \\ 0\leqx\leq0.5) f (\frac{}{} \frac{}{x,t,u,\partialu\partialx)} & =∂u∂x ( \end{matrix}$ 0.5≤x≤1)

${\begin{matrix} s (\frac{}{} x,t,u,\partialu\partialx) = -^{} & 1000eu ( \\ 0≤x≤0.5) s (\frac{}{}^{} & x,t,u,\partialu\partialx) = - eu \end{matrix}$ (0.5≤x≤1)

Теперь можно создать функцию для кодирования уравнения. Функция должна иметь подпись [c,f,s] = pdex2pde(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] = pdex2pde(x,t,u,dudx)
c = 1;
if x <= 0.5
    f = 5*dudx;
    s = -1000*exp(u);
else
    f = dudx;
    s = -exp(u);
end
end

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

Исходные условия кода

Затем запишите функцию, которая возвращает начальные условия. Начальное условие применяется при первом значении и предоставляет значение $u (x,_{}$ t0) для любого значения x. Используйте сигнатуру функцииu0 = pdex2ic(x) для записи функции.

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

$\begin{array}{l} u (x, 0) = 0 (0\leqx<1), u \\ (1,0) = 1 (x = 1). \end{array}$

Соответствующая функция:

function u0 = pdex2ic(x)
if x < 1
    u0 = 0;
else
    u0 = 1;
end
end

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

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

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

Поскольку этот пример имеет сферическую симметрию ( $m =$ 2), pdepe решатель автоматически применяет левое граничное условие для привязки решения в начале координат и игнорирует любые условия, указанные для левой границы в граничной функции. Поэтому для условия левой границы можно задать $_{pL} =_{} qL$ = 0. Для правого граничного условия можно переписать граничное условие в стандартной форме и считать значения коэффициентов $_{}$ для pR $_{}$ и qR.

Для $x =$ 1 уравнение равно $u (1, t) =1→ (\frac{u-1)}{}$ +0⋅∂u∂x=0. Коэффициенты:

$_{pR} (1, t, u)$ = u-1
$_{qR} (1, t)$ = 0

Граничная функция должна использовать сигнатуру функции [pl,ql,pr,qr] = pdex2bc(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] = pdex2bc(xl,ul,xr,ur,t)
pl = 0; 
ql = 0; 
pr = ur - 1;
qr = 0;
end

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

Пространственная сетка должна включать несколько значений около $x =$ 0,5 для учета прерывистого интерфейса, а также точки $около x$ = 1 из-за противоречивости начального $значения (u ($ 1,0) = 1) и граничного $значения (u$ (1, t) = 0) в этой точке. Решение быстро изменяется для малых t, поэтому используйте временной шаг, который может разрешить это резкое изменение.

x = [0 0.1 0.2 0.3 0.4 0.45 0.475 0.5 0.525 0.55 0.6 0.7 0.8 0.9 0.95 0.975 0.99 1];
t = [0 0.001 0.005 0.01 0.05 0.1 0.5 1];

Уравнение решения

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

m = 2;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,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 является матрицей 8 на 18, но в целом можно извлечь kКомпонент решения th с командой u = sol(:,:,k).

Извлечь первый компонент раствора из sol.

u = sol(:,:,1);

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

Создайте график поверхности решения $u$ , построенного в выбранных точках сетки для $x$ и t. Поскольку $m =$ 2, задача ставится в сферической геометрии со сферической симметрией, поэтому решение изменяется только в $радиальном$ направлении x.

surf(x,t,u)
title('Numerical solution with nonuniform mesh')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution u')

Figure contains an axes. The axes with title Numerical solution with nonuniform mesh contains an object of type surface.

Теперь постройте график только $x$ и $u$ , чтобы получить вид сбоку контуров на графике поверхности. Добавьте строку при $x =$ 0,5 для выделения эффекта интерфейса материала.

plot(x,u,x,u,'*')
line([0.5 0.5], [-3 1], 'Color', 'k')
xlabel('Distance x')
ylabel('Solution u')
title('Solution profiles at several times')

Figure contains an axes. The axes with title Solution profiles at several times contains 17 objects of type line.

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

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

function [c,f,s] = pdex2pde(x,t,u,dudx) % Equation to solve
c = 1;
if x <= 0.5
    f = 5*dudx;
    s = -1000*exp(u);
else
    f = dudx;
    s = -exp(u);
end
end
%----------------------------------------------
function u0 = pdex2ic(x) %Initial conditions
if x < 1
    u0 = 0;
else
    u0 = 1;
end
end
%----------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) % Boundary conditions
pl = 0;
ql = 0;
pr = ur - 1;
qr = 0;
end
%----------------------------------------------

См. также

pdepe

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