Решите УЧП с разрывом

В этом примере показано, как решить УЧП, который взаимодействует через интерфейс с материалом. Материальный интерфейс создает разрыв в проблеме в $x = 0.5$ , и начальное условие имеет разрыв на правильном контуре $x = 1$ .

Рассмотрите кусочный УЧП

${\begin{matrix} \frac{\partial u}{\partial t} = x^{- 2} \frac{\partial}{\partial x} (x^{2} 5 \frac{\partial u}{\partial x}) - 1000 e^{u} & (0 \leq x \leq 0.5) \\ \frac{\partial u}{\partial t} = x^{- 2} \frac{\partial}{\partial x} (x^{2} \frac{\partial u}{\partial x}) - e^{u} & (0.5 \leq x \leq 1) \end{matrix}$

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

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

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

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

Чтобы решить это уравнение в 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})$ .

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

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

$m = 2$

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

${\begin{matrix} f (x, t, u, \frac{\partial u}{\partial x}) = 5 \frac{\partial u}{\partial x} & (0 \leq x \leq 0.5) \\ f (x, t, u, \frac{\partial u}{\partial x}) = \frac{\partial u}{\partial x} & (0.5 \leq x \leq 1) \end{matrix}$

${\begin{matrix} s (x, t, u, \frac{\partial u}{\partial x}) = - 1000 e^{u} & (0 \leq x \leq 0.5) \\ s (x, t, u, \frac{\partial u}{\partial x}) = - e^{u} & (0.5 \leq x \leq 1) \end{matrix}$

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

x независимая пространственная переменная.
t независимая переменная времени.
u зависимая переменная, дифференцируемая относительно x и t.
dudx частичная пространственная производная $\partial u / \partial x$ .
Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения УЧП, ожидаемой pdepe. Эти коэффициенты закодированы в терминах входных переменных xTU, и 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, t_{0})$ для любого значения $x$ . Используйте функциональную подпись u0 = pdex2ic(x) записать функцию.

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

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

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

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

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

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

Поскольку этот пример имеет сферическую симметрию ( $m = 2$ ), pdepe решатель автоматически осуществляет левое граничное условие к связанному решение в начале координат и игнорирует любые условия, заданные для левого контура в граничной функции. Таким образом для левого граничного условия, можно задать $p_{L} = q_{L} = 0$ . Для правильного граничного условия можно переписать граничное условие в стандартной форме и прочитать содействующие значения для $p_{R}$ и $q_{R}$ .

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

$p_{R} (1, t, u) = u - 1$
$q_{R} (1, t) = 0$

Граничная функция должна использовать функциональную подпись [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t):

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

Выберите Solution Mesh

Пространственная mesh должна включать несколько значений рядом $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, уравнение УЧП, начальные условия, граничные условия и сетки для x и t.

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

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

u = sol(:,:,1);

Постройте решение

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

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

Теперь постройте только $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')

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

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

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