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

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

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

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

Похожие темы

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