Этот пример показывает, как решить систему дифференциальных уравнений с частными производными, которая использует функции шага в начальных условиях.
Рассмотрите PDE
Уравнения включают постоянные параметры , , , , и , и определены для и . Эти уравнения возникают в математической модели первых этапов ангиогенеза, связанного с опухолью [1]. представляет плотность камер эндотелиальных камер, и концентрация белка, который они выделяют в ответ на опухоль.
Эта задача имеет постоянное, устойчивое состояние, когда
.
Однако анализ устойчивости предсказывает эволюцию системы в неоднородное решение [1]. Таким образом, шаговые функции используются в качестве начальных условий для возмущения устойчивого состояния и стимулирования эволюции системы.
Граничные условия требуют, чтобы оба компонента решения имели нулевой поток и .
Чтобы решить эту систему уравнений в MATLAB, необходимо кодировать уравнения, начальные условия и граничные условия, затем выбрать подходящий mesh решения перед вызовом решателя pdepe
. Можно либо включить необходимые функции в качестве локальных функций в конце файла (как это сделано здесь), либо сохранить их как отдельные, именованные файлы в директории по пути MATLAB.
Прежде чем вы сможете кодировать уравнение, необходимо убедиться, что это в том виде, в котором pdepe
решатель ожидает:
Поскольку в системе PDE существует два уравнения, система PDE может быть переписана как
Значения коэффициентов в уравнении тогда
(только по диагонали)
Теперь можно создать функцию, чтобы кодировать уравнение. Функция должна иметь подпись [c,f,s] = angiopde(x,t,u,dudx)
:
x
- независимая пространственная переменная.
t
- независимая временная переменная.
u
- зависимая переменная, дифференцируемая по отношению к x
и t
. Это двухэлементный вектор, где u(1)
является и u(2)
является .
dudx
- частная пространственная производная . Это двухэлементный вектор, где dudx(1)
является и dudx(2)
является .
Выходные выходы c
, f
, и s
соответствуют коэффициентам в стандартной форме уравнения УЧП, ожидаемым по pdepe
.
В результате уравнения в этом примере могут быть представлены функцией:
function [c,f,s] = angiopde(x,t,u,dudx) d = 1e-3; a = 3.8; S = 3; r = 0.88; N = 1; c = [1; 1]; f = [d*dudx(1) - a*u(1)*dudx(2) dudx(2)]; s = [S*r*u(1)*(N - u(1)); S*(u(1)/(u(1) + 1) - u(2))]; end
(Примечание. Все функции включены в качестве локальных функций в конце примера.)
Затем напишите функцию, которая возвращает начальное условие. Начальное условие применяется в первый раз и предоставляет значение и для любого значения x. Используйте сигнатуру функции u0 = angioic(x)
для записи функции.
Эта задача имеет постоянное, устойчивое состояние, когда
.
Однако анализ устойчивости предсказывает эволюцию системы в неоднородное решение [1]. Таким образом, шаговые функции используются в качестве начальных условий для возмущения устойчивого состояния и стимулирования эволюции системы.
Соответствующая функция является
function u0 = angioic(x) u0 = [1; 0.5]; if x >= 0.3 && x <= 0.6 u0(1) = 1.05 * u0(1); u0(2) = 1.0005 * u0(2); end end
Теперь напишите функцию, которая оценивает граничные условия
Для задач, поставленных на интервале , граничные условия применяются для всех и либо или . Стандартная форма для граничных условий, ожидаемых решателем,
Для , граничное уравнение условия
Таким образом, коэффициенты:
.
Для граничные условия совпадают, поэтому и .
Граничная функция должна использовать сигнатуру функции [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t)
, где:
Входы xl
и ul
соответствуют и для левого контура.
Входы xr
и ur
соответствуют и для правого контура.
t
- независимая временная переменная.
Выходные выходы pl
и ql
соответствуют и для левого контура ( для этой задачи).
Выходные выходы pr
и qr
соответствуют и для правого контура ( для этой задачи).
Граничные условия в этом примере представлены функцией:
function [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t) pl = [0; 0]; ql = [1; 1]; pr = pl; qr = ql; end
Чтобы увидеть предельное поведение уравнений, требуется длительный временной интервал, поэтому используйте 10 точек в интервале . Кроме того, предельное распределение изменяется только примерно на 0,1% в течение интервала таким образом, подходящей является относительно тонкий пространственный mesh с 50 точками.
x = linspace(0,1,50); t = linspace(0,200,10);
Наконец, решите уравнение, используя симметрию , уравнение УЧП, начальное условие, граничные условия и сетки для и .
m = 0; sol = pdepe(m,@angiopde,@angioic,@angiobc,x,t);
pdepe
возвращает решение в трехмерный массив sol
, где sol(i,j,k)
аппроксимирует k
первый компонент решения оценивается в t(i)
и x(j)
. Извлеките компоненты решения в отдельные переменные.
n = sol(:,:,1); c = sol(:,:,2);
Создать объемную поверхностную диаграмму компонентов решения и нанесенный на график в выбранных точках mesh для и .
surf(x,t,c) title('c(x,t): Concentration of Fibronectin') xlabel('Distance x') ylabel('Time t')
surf(x,t,n) title('n(x,t): Density of Endothelial Cells') xlabel('Distance x') ylabel('Time t')
Теперь постройте график только окончательных распределений решений в . Эти графики соответствуют Фигурам 3 и 4 в [1].
plot(x,n(end,:))
title('Final distribution of n(x,t_f)')
plot(x,c(end,:))
title('Final distribution of c(x,t_f)')
[1] Хамфрис, M.E. и M.A.J. Капеллан. «Математическая модель первых шагов ангиогенеза, связанного с опухолью: Образование капиллярного ростка и вторичное ответвление». IMA Journal of Mathematics Applied in Medicine & Biology, 13 (1996), pp. 73-98.
Здесь перечислены локальные вспомогательные функции, которые УЧП решателя pdepe вызывает для вычисления решения. Также можно сохранить эти функции как собственные файлы в директории по пути MATLAB.
function [c,f,s] = angiopde(x,t,u,dudx) % Equation to solve d = 1e-3; a = 3.8; S = 3; r = 0.88; N = 1; c = [1; 1]; f = [d*dudx(1) - a*u(1)*dudx(2) dudx(2)]; s = [S*r*u(1)*(N - u(1)); S*(u(1)/(u(1) + 1) - u(2))]; end % --------------------------------------------- function u0 = angioic(x) % Initial Conditions u0 = [1; 0.5]; if x >= 0.3 && x <= 0.6 u0(1) = 1.05 * u0(1); u0(2) = 1.0005 * u0(2); end end % --------------------------------------------- function [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t) % Boundary Conditions pl = [0; 0]; ql = [1; 1]; pr = pl; qr = ql; end % ---------------------------------------------