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

Прежде чем вы сможете кодировать уравнение, вы должны убедиться, что это в форме, которую pdepe решатель ожидает. Стандартная форма, которая pdepe ожидает, что

$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{-\mathit{m}}\frac{\partial }{\partial \mathit{x}}\left({\mathit{x}}^{\mathit{m}}\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\right)+\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$.

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

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

Теперь можно создать функцию, чтобы кодировать уравнение. Функция должна иметь подпись [c,f,s] = pdex2pde(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

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

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

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

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

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

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

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

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

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

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

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

Граничные условия в этом примере представлены функцией:

function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = 0; 
ql = 0; 
pr = ur - 1;
qr = 0;
end

Выберите Mesh решения

Пространственный mesh должна включать несколько значений около $\mathit{x}=0.5$ для расчета прерывистого интерфейса, а также точек около $\mathit{x}=1$ из-за несогласованного начального значения ($\mathit{u}\left(1,0\right)=1$) и краевое значение ($\mathit{u}\left(1,\mathit{t}\right)=0$) в эту точку. Решение быстро изменяется для малых $\mathit{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первый компонент решения ${\mathit{u}}_{\mathit{k}}$ оценивается в t(i) и x(j). Размер sol является length(t)-by- length(x)-by- length(u0), с u0 задает начальное условие для каждого компонента решения. Для этой задачи u имеет только один компонент, так что sol является матрицей 8 на 18, но в целом можно извлечь kкомпонент решения с командой u = sol(:,:,k).

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

u = sol(:,:,1);

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

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

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

Теперь постройте график просто $\mathit{x}$ и $\mathit{u}$ для получения вида сбоку контуров на объемной поверхностной диаграмме. Добавьте линии в $\mathit{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')

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

Здесь перечислены локальные вспомогательные функции, которые решатель 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
%----------------------------------------------