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

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

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

{ut=x-2x(x25ux)-1000eu(0x0.5)ut=x-2x(x2ux)-eu(0.5x1)

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

u(x,0)=0  (0x<1),u(1,0)=1  (x=1).

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

       ux=0  (x=0),u(1,t)=1  (x=1).

Чтобы решить это уравнение в MATLAB, необходимо закодировать уравнение, начальные условия и граничные условия, затем выбрать подходящую mesh решения прежде, чем вызвать решатель pdepe. Вы любой может включать необходимые функции как локальные функции в конце файла (как сделано здесь) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

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

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

c(x,t,u,ux)ut=x-mx(xmf(x,t,u,ux))+s(x,t,u,ux).

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

{ut=x-2x(x25ux)-1000eu(0x0.5)ut=x-2x(x2ux)-eu(0.5x1)

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

m=2

c(x,t,u,ux)=1

{f(x,t,u,ux)=5ux(0x0.5)f(x,t,u,ux)=ux(0.5x1)

{s(x,t,u,ux)=-1000eu(0x0.5)s(x,t,u,ux)=-eu(0.5x1)

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

  • x независимая пространственная переменная.

  • t независимая переменная времени.

  • u зависимая переменная, дифференцируемая относительно x и t.

  • dudx частичная пространственная производная u/x.

  • Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения PDE, ожидаемой 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,t0) для любого значения x. Используйте функциональную подпись u0 = pdex2ic(x) записать функцию.

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

u(x,0)=0        (0x<1),u(1,0)=1        (x=1).

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

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

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

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

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

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

Для x=1, уравнение u(1,t)=1(u-1)+0ux=0. Коэффициенты:

  • pR(1,t,u)=u-1

  • qR(1,t)=0

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

  • Входные параметры xl и ul соответствуйте вам и x для левого контура.

  • Входные параметры xr и ur соответствуйте вам и 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

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

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

pdepe возвращает решение в трехмерном массиве sol, где sol(i,j,k) аппроксимирует kкомпонент th решения 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 построенный в выбранной 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')

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.

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

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

Смотрите также

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте