Решите систему УЧП с начальными ступенчатыми функциями условия

Скрипт Open Live Script

В этом примере показано, как решить систему дифференциальных уравнений с частными производными, которая использует ступенчатые функции в начальных условиях.

Рассмотрите УЧП

$\begin{array}{l} \frac{\partial n}{\partial t} = \frac{\partial}{\partial x} [d \frac{\partial n}{\partial x} - a n \frac{\partial c}{\partial x}] + S r n (N - n), \\ \frac{\partial c}{\partial t} = \frac{\partial^{2} c}{\partial x^{2}} + S (\frac{n}{n + 1} - c) . \end{array}$

Уравнения включают постоянные параметры $d$ , $a$ , $S$ , $r$ , и $N$ , и заданы для $0 \leq x \leq 1$ и $t \geq 0$ . Эти уравнения возникают в математической модели первых шагов связанного с опухолью ангиогенеза [1]. $n (x, t)$ представляет плотность ячейки эндотелиальных клеток, и $c (x, t)$ концентрация белка они выпускают в ответ на опухоль.

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

$[\begin{matrix} n_{0} \\ c_{0} \end{matrix}] = [\begin{matrix} 1 \\ 0.5 \end{matrix}]$ .

Однако анализ устойчивости предсказывает эволюцию системы к неоднородному решению [1]. Таким образом, ступенчатые функции используются в качестве начальных условий, чтобы встревожить устойчивое состояние и стимулировать эволюцию системы.

Граничные условия требуют, чтобы оба компонента решения имели нулевой поток в $x = 0$ и $x = 1$ .

$\begin{array}{l} \frac{\partial}{\partial x} n (0, t) = \frac{\partial}{\partial x} n (1, t) = 0, \\ \frac{\partial}{\partial x} c (0, t) = \frac{\partial}{\partial x} c (1, t) = 0 . \end{array}$

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

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

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

$c (x, t, u, \frac{\partial u}{\partial x}) \frac{\partial u}{\partial t} = x^{- m} \frac{\partial}{\partial x} (x^{m} f (x, t, u, \frac{\partial u}{\partial x})) + s (x, t, u, \frac{\partial u}{\partial x}) .$

С тех пор существует два уравнения в системе УЧП, система УЧП может быть переписана как

$[\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}] \frac{\partial}{\partial t} [\begin{matrix} n \\ c \end{matrix}] = \frac{\partial}{\partial x} [\begin{matrix} d \frac{\partial n}{\partial x} - a n \frac{\partial c}{\partial x} \\ \frac{\partial c}{\partial x} \end{matrix}] + [\begin{matrix} S r n (N - n) \\ S (\frac{n}{n + 1} - c) \end{matrix}] .$

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

$m = 0$

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

$f (x, t, u, \frac{\partial u}{\partial x}) = [\begin{matrix} d \frac{\partial n}{\partial x} - a n \frac{\partial c}{\partial x} \\ \frac{\partial c}{\partial x} \end{matrix}]$

$s (x, t, u, \frac{\partial u}{\partial x}) = [\begin{matrix} S r n (N - n) \\ S (\frac{n}{n + 1} - c) \end{matrix}]$

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

x независимая пространственная переменная.
t независимая переменная времени.
u зависимая переменная, дифференцируемая относительно x и t. Это - двухэлементный вектор где u(1) $n (x, t)$ и u(2) $c (x, t)$ .
dudx частичная пространственная производная $\partial u / \partial x$ . Это - двухэлементный вектор где dudx(1) $\partial n / \partial x$ и dudx(2) $\partial c / \partial x$ .
Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения PDE, ожидаемой 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

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

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

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

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

$[\begin{matrix} n_{0} \\ c_{0} \end{matrix}] = [\begin{matrix} 1 \\ 0.5 \end{matrix}]$ .

Однако анализ устойчивости предсказывает эволюцию системы к негомогенному решению [1]. Так, ступенчатые функции используются в качестве начальных условий, чтобы встревожить устойчивое состояние и стимулировать эволюцию системы.

$\begin{array}{l} u (x, 0) = [\begin{matrix} n_{0} \\ c_{0} \end{matrix}], \\ u (x, 0) = {\begin{matrix} 1.05 u_{1} & 0.3 \leq x \leq 0.6 \\ 1.0005 u_{2} & 0.3 \leq x \leq 0.6 \end{matrix} \end{array}$

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

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

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

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

$\begin{array}{l} \frac{\partial}{\partial x} n (0, t) = \frac{\partial}{\partial x} n (1, t) = 0, \\ \frac{\partial}{\partial x} c (0, t) = \frac{\partial}{\partial x} c (1, t) = 0 . \end{array}$

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

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

Для $x = 0$ , уравнение граничного условия

$[\begin{array}{c} 0 \\ 0 \end{array}] + [\begin{array}{c} 1 \\ 1 \end{array}] \cdot [\begin{array}{c} d \frac{\partial n}{\partial x} - a n \frac{\partial c}{\partial x} \\ \frac{\partial c}{\partial x} \end{array}] = 0 .$

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

$p_{L} (x, t, u) = [\begin{matrix} 0 \\ 0 \end{matrix}],$
$q_{L} (x, t) = [\begin{matrix} 1 \\ 1 \end{matrix}]$ .

Для $x = 1$ граничные условия являются тем же самым, таким образом, $p_{R} (x, t, u) = [\begin{matrix} 0 \\ 0 \end{matrix}]$ и $q_{R} (x, t) = [\begin{matrix} 1 \\ 1 \end{matrix}]$ .

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

Входные параметры xl и ul соответствовать $u$ и $x$ для левого контура.
Входные параметры xr и ur соответствовать $u$ и $x$ для правильного контура.
t независимая переменная времени.
Выходные параметры pl и ql соответствовать $p_{L} (x, t, u)$ и $q_{L} (x, t)$ для левого контура ( $x = 0$ для этой проблемы).
Выходные параметры pr и qr соответствовать $p_{R} (x, t, u)$ и $q_{R} (x, t)$ для правильного контура ( $x = 1$ для этой проблемы).

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

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

Выберите Solution Mesh

Долговременный интервал требуется, чтобы видеть ограничивающее поведение уравнений, так используйте 10 точек в интервале $0 \leq t \leq 200$ . Кроме того, предельное распределение $c (x, t)$ варьируется только на приблизительно 0,1% на интервале $0 \leq x \leq 1$ , таким образом, относительно прекрасная пространственная mesh с 50 точками является соответствующей.

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

Решите уравнение

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

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

pdepe возвращает решение в трехмерном массиве sol, где sol(i,j,k) аппроксимирует kкомпонент th решения $u_{k}$ оцененный в t(i) и x(j). Извлеките компоненты решения в отдельные переменные.

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

Постройте решение

Создайте объемную поверхностную диаграмму компонентов решения $n$ и $c$ построенный в выбранной mesh указывает для $x$ и $t$ .

surf(x,t,c)
title('c(x,t): Concentration of Fibronectin')
xlabel('Distance x')
ylabel('Time t')

Figure contains an axes. The axes with title c(x,t): Concentration of Fibronectin contains an object of type surface.

surf(x,t,n)
title('n(x,t): Density of Endothelial Cells')
xlabel('Distance x')
ylabel('Time t')

Figure contains an axes. The axes with title n(x,t): Density of Endothelial Cells contains an object of type surface.

Теперь постройте только итоговые распределения решений в $t_{f} = 200$ . Эти графики соответствуют рисункам 3 и 4 в [1].

plot(x,n(end,:))
title('Final distribution of n(x,t_f)')

Figure contains an axes. The axes with title Final distribution of n(x,t_f) contains an object of type line.

plot(x,c(end,:))
title('Final distribution of c(x,t_f)')

Figure contains an axes. The axes with title Final distribution of c(x,t_f) contains an object of type line.

Ссылки

[1] Humphreys, M.E. и M.A.J. Капеллан. "Математическая модель первых шагов связанного с опухолью ангиогенеза: Капиллярное формирование ростка и вторичное ответвление". Журнал IMA Математики, Прикладной в Medicine & Biology, 13 (1996), стр 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
% ---------------------------------------------

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

pdepe

Документация