Кодовое уравнение

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

Поскольку в системе ПДЭ существует два уравнения, система ПДЭ может быть переписана как

Значения коэффициентов в уравнении затем

$$\mathit{c}\mathit{(}\mathit{}\mathit{}\frac{\mathit{}}{\mathit{}}x,t,u,\partial u\partial x)\begin{array}{c}\mathrm{=}\\ \mathrm{}\end{array}[$$11] (только диагональные значения)

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

В результате уравнения в этом примере могут быть представлены функцией:

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

(Примечание: Все функции включены в качестве локальных функций в конце примера.)

Исходные условия кода

Затем запишите функцию, которая возвращает начальное условие. Начальное условие применяется при первом значении и предоставляет значение $\mathit{n}\mathit{(}x{\mathit{,}}_{\mathrm{}}$t0) $\mathit{и}\mathit{c}{\mathit{(}}_{\mathrm{x}},$ t0) для любого значения x. Используйте сигнатуру функцииu0 = angioic(x) для записи функции.

Эта проблема имеет постоянное, устойчивое состояние, когда

$\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\end{array}[n0c0]=\begin{array}{c}\mathrm{}\\ \mathrm{}\mathrm{}\end{array}[10.5]$.

Однако анализ стабильности предсказывает эволюцию системы в негомогенное решение [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

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

Теперь напишите функцию, которая вычисляет граничные условия

Для проблем, возникающих в интервале $\mathit{}\mathit{}\mathit{a\le x\le b}$, граничные условия применяются для всех $\mathit{t}$ и $\mathit{либо\; x}\mathit{=}$ a, $\mathit{либо}\mathit{x}$= b. Стандартная форма для граничных условий, ожидаемых решателем:

Для $\mathit{x}\mathrm{=}$0 уравнение граничного условия

Таким образом, коэффициенты:

Для $\mathit{x}\mathrm{=}$ 1 граничные условия одинаковы, ${\mathit{так\; что}}_{\mathit{}}\mathit{pR}\mathit{(}x\mathit{,}t,\begin{array}{c}\mathrm{}\\ \mathrm{u}\end{array})$ = [${\mathit{00}}_{\mathit{]}}\mathit{}\mathit{и}qR\begin{array}{c}\mathrm{}\\ \mathrm{(}\end{array}x$, t) = [11].

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

Граничные условия в этом примере представлены функцией:

function [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t)
pl = [0; 0];
ql = [1; 1];
pr = pl;
qr = ql;
end

Выбрать сетку решения

Для просмотра ограничивающего поведения уравнений требуется длительный временной интервал, поэтому используйте 10 точек в $\mathrm{}\mathit{}\mathrm{0\le t\le 200\; интервала}$. Кроме того, предельное распределение $\mathit{c}\mathit{(}x\mathit{,}$t) изменяется только примерно на 0,1% по $\mathrm{}\mathit{}\mathrm{}$0≤x≤1 интервалу, поэтому подходящая относительно тонкая пространственная сетка с 50 точками.

x = linspace(0,1,50);
t = linspace(0,200,10);

Уравнение решения

Наконец, решите уравнение, используя симметрию $\mathit{m}$, уравнение PDE, начальное условие, граничные условия и сетки для $\mathit{x}$ и $\mathit{t}$.

m = 0;
sol = pdepe(m,@angiopde,@angioic,@angiobc,x,t);

pdepe возвращает решение в массиве 3-D sol, где sol(i,j,k) аппроксимирует kтретий компонент раствора ${\mathit{}}_{\mathit{uk}}$ оценен на t(i) и x(j). Извлеките компоненты решения в отдельные переменные.

n = sol(:,:,1);
c = sol(:,:,2);

Решение для построения графика

Создайте график поверхности компонентов решения $\mathit{n}$ и $\mathit{c}$, построенный в выбранных точках сетки для $\mathit{x}$ и $\mathit{t}$.

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')

Теперь постройте график только конечных распределений решений при ${\mathit{}}_{\mathit{tf}}\mathrm{=}$200. Эти графики соответствуют рисункам 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] Хамфрис, М.Е. и М.А.Дж. Капеллан. «Математическая модель первых шагов ангиогенеза, связанного с опухолью: образование капиллярных ростков и вторичное разветвление». IMA Journal of Mathematics Applied in Medicine & Biology, 13 (1996), стр. 73-98.

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

Здесь перечислены локальные вспомогательные функции, которые вызывает pdepe решателя PDE для вычисления решения. Кроме того, эти функции можно сохранить в виде собственных файлов в каталоге по пути 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
% ---------------------------------------------