Уравнение кода

Прежде чем можно будет закодировать уравнение, необходимо убедиться что именно в форме 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

Выберите Solution 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, уравнение PDE, начальные условия, граничные условия и сетки для x и t.

m = 2;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

pdepe возвращает решение в трехмерном массиве sol, где sol(i,j,k) аппроксимирует kкомпонент th решения ${\mathit{u}}_{\mathit{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);

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

Создайте объемную поверхностную диаграмму решения $\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')

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

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