В этом примере показано, как решить УЧП, который взаимодействует через интерфейс с материалом. Материальный интерфейс создает разрыв в проблеме в , и начальное условие имеет разрыв на правильном контуре .
Рассмотрите кусочный УЧП
Начальные условия
Граничные условия
Чтобы решить это уравнение в MATLAB, необходимо закодировать уравнение, начальные условия и граничные условия, затем выбрать подходящую mesh решения прежде, чем вызвать решатель pdepe
. Вы любой может включать необходимые функции как локальные функции в конце файла (как сделано здесь) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.
Прежде чем можно будет закодировать уравнение, необходимо убедиться что именно в форме pdepe
решатель ожидает. Стандартная форма, что pdepe
ожидает
.
В этом случае УЧП находится в соответствующей форме, таким образом, можно прочитать значения коэффициентов.
Значения для термина потока и характеристики выброса изменитесь в зависимости от значения . Коэффициенты:
Теперь можно создать функцию, чтобы закодировать уравнение. Функция должна иметь подпись [c,f,s] = pdex2pde(x,t,u,dudx)
:
x
независимая пространственная переменная.
t
независимая переменная времени.
u
зависимая переменная, дифференцируемая относительно x
и t
.
dudx
частичная пространственная производная .
Выходные параметры c
F
, и s
соответствуйте коэффициентам в стандартной форме уравнения УЧП, ожидаемой pdepe
. Эти коэффициенты закодированы в терминах входных переменных 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
(Примечание: Все функции включены как локальные функции в конце примера.)
Затем запишите функцию, которая возвращает начальные условия. Начальное условие применяется в первой временной стоимости и вводит значение для любого значения . Используйте функциональную подпись u0 = pdex2ic(x)
записать функцию.
Начальные условия
Соответствующая функция
function u0 = pdex2ic(x) if x < 1 u0 = 0; else u0 = 1; end end
Теперь запишите функцию, которая оценивает граничные условия. Для проблем, созданных на интервале , граничные условия применяются ко всему t и также или . Стандартная форма для граничных условий, ожидаемых решателем,
Поскольку этот пример имеет сферическую симметрию (), pdepe
решатель автоматически осуществляет левое граничное условие к связанному решение в начале координат и игнорирует любые условия, заданные для левого контура в граничной функции. Таким образом для левого граничного условия, можно задать . Для правильного граничного условия можно переписать граничное условие в стандартной форме и прочитать содействующие значения для и .
Для , уравнение . Коэффициенты:
Граничная функция должна использовать функциональную подпись [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
:
Входные параметры xl
и ul
соответствуйте вам и x для левого контура.
Входные параметры xr
и ur
соответствуйте вам и x для правильного контура.
t
независимая переменная времени.
Выходные параметры pl
и ql
соответствовать и для левого контура ( для этой проблемы).
Выходные параметры pr
и qr
соответствовать и для правильного контура ( для этой проблемы).
Граничные условия в этом примере представлены функцией:
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) pl = 0; ql = 0; pr = ur - 1; qr = 0; end
Пространственная mesh должна включать несколько значений рядом составлять прерывистый интерфейс, а также точки рядом из-за противоречивого начального значения () и граничное значение () в той точке. Решение изменяется быстро для маленького , так используйте временной шаг, который может разрешить это резкое изменение.
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 решения оцененный в 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);
Создайте объемную поверхностную диаграмму решения построенный в выбранной mesh указывает для и . С тех пор проблема создана в сферической геометрии со сферической симметрией, таким образом, решение только изменяется в шине с радиальным кордом направление.
surf(x,t,u) title('Numerical solution with nonuniform mesh') xlabel('Distance x') ylabel('Time t') zlabel('Solution u')
Теперь постройте только и получить вид сбоку контуров в объемной поверхностной диаграмме. Добавьте линию в подсвечивать эффект материального интерфейса.
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 %----------------------------------------------