Построение, численное моделирование и анализ одномерной модели гемодинамики

 

Содержание


Введение

1. Основная часть

1.1 Вывод модели

1.2 Анализ математических свойств модели

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

1.4 Реализация модели

1.4.1 Сведение к обыкновенному дифференциальному уравнению

1.4.2 Линеаризация граничных условий

1.4.3 Метод ортогональной прогонки

1.4.4 Теория возмущений (замечание о точности решения)

1.4.5 Генерация матриц граничных условий для метода ортогональной прогонки для системы из 55 артерий

1.4.6 Описание работы программы

1.4.7 Описание структуры программы

1.5 Тестовые расчёты

1.5.1 Моделирование гемодинамики в норме

1.5.2 Моделирование физической нагрузки и тахикардии

1.5.3 Применение метода итераций для оптимизации тестовых расчётов

Заключение

Список литературы

Приложение


Введение


Сердечно-сосудистая система (ССС) человека состоит из сердца и сосудов - артерий, капилляров и вен. Основная функция ССС состоит в обеспечении продвижения крови по замкнутой цепи сосудов. Назначение постоянной циркуляции крови в организме заключается в доставке к различным органам и тканям веществ, необходимых для их нормального функционирования, и удалении продуктов их жизнедеятельности. Артерии доставляют обогащённую кислородом кровь из левого желудочка сердца к различным органам, а по венам насыщенная углекислым газом кровь возвращается к правому предсердию и направляется в малый круг кровообращения для газообмена.

Артериальное давление (АД) является важным параметром, характеризующим работу кровеносной системы. Нормальные численные значения АД лежат в пределах 120-135 мм рт ст для систолического и 75-85 мм рт. Ст. для диастолического в режиме бодрствования. Стойкое превышение этих значений характеризуется как артериальная гипертония (АГ), которая является одним из самых распространённых заболеваний ССС, и без должного и вовремя начатого лечения может привести к поражению различных органов человека. АГ занимает первое место по смертности в большинстве развитых стран мира, поэтому разработка мер и средств по ее предупреждению, лечению и диагностике является одной из важных задач современной науки.

Моделирование движения крови по организму человека является важной задачей. Созданием моделей процессов, происходящих в ССС, люди начали заниматься достаточно давно. Началом формального описания физиологии кровообращения можно считать работу Уильяма Гарвея "Анатомическое исследование о движении сердца и крови у животных (Exercitatio anatomica de motu cordis et sanguinis in animalibus), выпущенной в 1628 году. Он заключил, что кровь в организме человека движется по кругу. В этой книге Гарвей точно описал работу сердца и выделил малый и большой круги кровообращения. [6] Наряду с этим, Гарвей доказал, что сердце ритмически бьется до тех пор, пока в организме теплится жизнь, причем после каждого сокращения сердца наступает короткий перерыв в его работе, во время которого этот важный орган отдыхает (таким образом он описал систолу и диастолу).

Математические модели стали появляться гораздо позже. Одномерное моделирование артериальной системы человека было предложено Эйлером в 1775 году, который вывел частные дифференциальные уравнения, выражающие закон сохранения массы и закон сохранения импульса для идеальной (невязкой) жидкости. Для того, чтобы замкнуть систему уравнений, он предложил два возможных, но экспериментально нереалистичных, уравнения состояния, описывающих поведение эластичной стенки в зависимости от изменений давления в сосуде. Вероятно, он не определил природу потока, как волновую, и не мог найти решения своих уравнений, испытывая, как он сам сказал, "непреодолимые трудности". Волновая природа артериального течения была впервые научно описана Юнгом, который вывел волновую скорость, используя аргументы, базирующиеся на интуиции и аналогии с ньютоновской теорией скорости звука в воздухе [3]. Уравнения, выведенные Эйлером в 1775 году - это система нелинейных дифференциальных уравнений в частных производных, аналогичная волновому уравнению для мелкой воды в гидродинамике или уравнению, описывающему одномерное невязкое течение в газовой динамике.

Следует отметить, что физиологические параметры артериальной системы человека таковы, что большинство характеристик потока сохраняются при линеаризации системы уравнений. Это, главным образом, было использовано в методе Вумерслея (Womersley, 1957), который линеаризовал двумерные уравнения потока в прямых, круглых эластичных трубках и получил волновое решение с помощью преобразования Фурье.

Развивая математический аппарат, предложенный Эйлером, Навье (1822) и Стокс (1845) вывели дифференциальные уравнения, описывающие движения вязкой несжимаемой жидкости (уравнения Навье-Стокса). В 1842 году физик и врач Пуазейль связал расход жидкости с давлением. В начале XX в Франк предложил идею, что система кровообращения аналогична электрической цепи, что положило начало новому способу представления ССС, используемому до настоящего времени.

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

Используя предположение, что плотность крови постоянна, физические процессы, происходящие в артериальной системе человека, можно описать трёхмерными нестационарными уравнениями Навье-Стокса, представляющими собой систему дифференциальных уравнений второго порядка в частных производных, совместно с уравнениями динамики эластичных оболочек сосудов, предложенными Пуазейлем [4]. Аналитическое решение этих уравнений найдено лишь для некоторых частных случаев. Применение такой сложной многомерной модели требует нахождения численных решений задач со свободной границей для уравнения Навье-Стокса в сложных областях. Это приводит к большим вычислительным затратам, поэтому многомерные модели на практике для глобального описания артериальной системы не применяются, а используются только для детального описания кровотока в интересующих исследователя локальных областях. Более того, часто нет необходимости в детальном представлении процесса гемодинамики для дальнейшего анализа медиками-специалистами, поэтому даже упрощённые модели могут обеспечить достаточным объемом информации при разумных вычислительных затратах. Таким образом, в настоящее время для глобального описания кровотока используют различные приближения исходной трёхмерной модели. Модели кровообращения, описывающие гемодинамику системы "в целом", являются важнейшим инструментом для исследования взаимодействия гемодинамических параметров, изучения и диагностики повреждений ССС, оценки последствий тех или иных методов вмешательства в её функционирование и т.д.

В настоящее время в моделировании гемодинамики используются несколько моделей, отличающихся степенью детальности описания системы. Наиболее простыми (и характеризующимися "грубым" приближением реальной трехмерной модели ССС) являются модели с сосредоточенными параметрами (или "нульмерные" (0D) модели). Главной их особенностью является сопоставление кровеносной системе электрических цепей. При таком подходе аналогами давления и объемного кровотока являются напряжение и электрический ток, а аналогами сосудов - сопротивления отдельных сегментов электрической цепи.0D модели являются упрощением одномерных (1D) моделей гемодинамики. Одномерные модели являются гидравлическим приближением исходной математической 3D модели, и получаются её усреднением по "поперечному" направлению. Третья группа упрощённых моделей - это гибридные модели. Они представляют собой стыковку, например, 3D и 0D моделей, 3D и 1D моделей и т.д.

В настоящее время наиболее интересной для анализа представляется 1D модель в силу предоставляемой ею возможности проследить за изменением параметров на всей протяженности сосуда по сравнению с 0D моделью, а также ввиду большей простоты моделирования при достаточной детальности описания системы по сравнению с 2D и 3D моделями. Для решения уравнений 1D модели могут быть применены различные вычислительные методы, характеризующиеся различной точностью, устойчивостью к погрешностям, скоростью сходимости итераций и простотой использования.

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

1. Основная часть


1.1 Вывод модели


Одномерная модель описывает течение крови в артериях и ее взаимодействие с подвижными стенками. Пусть t - время, а (x,y,z) - декартовы координаты. Будем считать, что идеализацией артерии является эластичная трубка в виде прямого кругового цилиндра. Рассмотрим кусок этой трубки длины L=const от z=0 до z=L (Рис.1.1).


Рис.1.1 Упрощённая геометрия кровеносного сосуда.


Через D (t) мы будем обозначать пространственную область, которая является указанной трубкой, заполненной кровью. Эта область меняется со временем под воздействием пульсирующей жидкости (крови). Осевое сечение {z=const} цилиндра D (t) будем обозначать через S=S (t,z). Описанная упрощенная геометрия артерии естественным образом приводит нас к использованию цилиндрических координат (r, ?, z).

Основные уравнения выводятся при следующих упрощающих предположениях:

) Осевая симметрия. Все величины не зависят от угловой координаты ?. Как следствие, каждое осевое сечение S остается круговым в течении всего времени движения стенок сосуда. Радиус трубки R является функцией времени t и осевой координаты z.

2) Радиальные смещения. Стенка сосуда перемещается только вдоль радиального направления, т.е. в каждой точке поверхности трубки вектор смещения


? = ?,


где ? = R ? смещение относительно исходного радиуса .

3) Фиксированная ось цилиндра. Сосуд растягивается и сжимается вдоль оси цилиндра, которая фиксирована во времени. (Эта гипотеза согласуется с гипотезой об осевой симметрии. Однако, она исключает возможность учета эффектов смещения оси артерии, имеющих, например, место в коронарных артериях и вызванных движением сердца.)

) Постоянство давления на каждом осевом сечении. Предполагается, что давление p постоянно на каждом сечении, т.е. оно зависит только от z и t.

) Отсутствие массовых сил. Мы пренебрегаем массовыми силами. (Xотя, например, учет силы гравитации не создает больших трудностей и связан только лишь с добавкой к давлению слагаемого вида gh. Немного сложнее учесть изменение гравитационных сил, возникающего при вставании человека из горизонтального положения.)

) Преобладание осевой компоненты скорости. Считается, что компоненты скорости, ортогональные оси z, являются пренебрежимо малыми по сравнению с осевой компонентой скорости . Предполагается, что


(1.1)


где - средняя скорость на каждом осевом сечении, а v: [0, 1] ? R - профиль скорости. (Тот факт, что профиль скорости явно не зависит от времени и пространственных переменных, вообще говоря, противоречит экспериментальным наблюдениям, а также численным результатам для полной трехмерной модели. Однако, мы можем считать, что v является профилем скорости для усредненной конфигурации течения.)

Через A=A (t,z) мы будем обозначать площадь осевого сечения S (t,z). Она находится по формуле


. (1.2)


Средняя скорость


(1.3)


где


(1.4)


средний объёмный кровоток. Из (1.1) следует, что


(1.5)


Введём в рассмотрение коэффициент


(1.6)


Который называют коэффициентом Кориолиса. Используя неравенство Коши-Буняковского, из (1.1) мы получаем, что ? ? 1. В общем случае этот коэффициент является величиной, зависящей от времени и пространственных переменных, но в рамках нашей модели ? = const, что является следствием упрощающего предположения (1.1).

В качестве возможного профиля скорости можно взять параболический профиль что соответствует известному решению Пуазейля, характерному для стационарных течений в круговых трубках. В этом случае ? = 4/3. Однако, известно, что для кровотока в артериях профиль скорости в среднем достаточно "плоский. Как правило, он описывается функцией где обычно ? = 9 (величина ? = 2 соответствует параболическому профилю). Тогда ? = (? + 2) / (? + 1) = 1.1 Более того, выбор ? = 1, соответствующий "абсолютно плоскому профилю скорости (полагая формально ? = ?, имеем v (y) = 1 для и v (1) = 0), значительно упрощает последующий анализ модели. Заметим, что в большинстве работ (см., например, [7]) по умолчанию рассматриваются модели с ? = 1.

Одномерная модель выводится с помощью интегрирования уравнений Навье-Стокса по произвольному осевому сечению S, т.е. усреднением этих уравнений по "поперечному" направлению.

Одномерная модель представляется системой из двух дифференциальных уравнений с частными производными:


(1.7)


где A, Q - неизвестные,

? - коэффициент Кориолиса,

? - плотность,

- коэффициент трения, который зависит от выбранного профиля скорости, т.е. от выбора функции v в (1.1).

Система, рассматривается при для интервала времени (0,T).

Для того чтобы замкнуть систему надо задать давление. Для задания закона давления необходимо учесть движение стенки кровеносного сосуда (в нашем случае артерии). Для описания упругих напряжений стенки сосуда, на практике чаще всего используется линейный закон (закон Гука)


,


где , E - модуль Юнга, - толщина стенки сосуда.

Тогда с учётом (1.2) получаем


где , .


Мы рассматриваем одномерную модель гемодинамики (1.7) с учетом обычно используемой на практике упрощающей гипотезы, которая состоит в предположении о пренебрежимой малости инерциальных членов по сравнению с упругими напряжениями. При таком предположении закон давления для системы (1.7) имеет вид:


(1.8)


И система (1.7), (1.8) является системой уравнений в частных производных первого порядка.

Рассматривается случай ? = 1, (для параболического профиля)

? - постоянный коэффициент вязкости.

Случай ? = 1 соответствует профилю скорости с ? = ?. Для этого профиля скорости коэффициент трения . Заметим, что ? ? при ? ? ?. Поэтому, уравнения (1.7), (1.8) с = 8?? являются некоторым огрублением общей одномерной модели гемодинамики. С другой стороны, модель (1.7), (1.8) с = 8?? часто вполне адекватно отражает реальные физические процессы, протекающие в артериальной системе человека. Поэтому и мы численные эксперименты проводим именно для этой модели, т.е. для уравнений (1.7) с ? = 1 и = 8?? (и используем закон давления (1.8)).


1.2 Анализ математических свойств модели


Для простоты будем предполагать, что и (т.е. , и модуль Юнга E = const). С учетом (1.8) система уравнений (1.7) записывается в квазилинейном виде


(2.1)


для вектора неизвестных, где


, .


При естественном требовании A>0 система (2.1) является одномерной квазилинейной гиперболической системой уравнений. Матрица В имеет два вещественных и различных собственных значения


, (2.2)


где



(). Причём при стандартных физиологических условиях типичные значения скорости крови u = Q/A и механических характеристик стенок сосуда таковы, что , т.е. и . (2.3)

Таким образом, гиперболическая система (2.1) на интервале (0,L) имеет по одной уходящей характеристике на границах z=0 и z=L и поэтому она требует по одному граничному условию на каждой границе.

Систему (2.1) на гладких решениях можно переписать также в консервативном виде, точнее в виде гиперболической системы балансовых законов


, (2.4)


где


.


Обсудим теперь кратко вопрос о приведении системы (2.1) к каноническому виду. Заметим, что гиперболическая система двух уравнений всегда может быть приведена к каноническому виду локально, т.е. в достаточно малой окрестности произвольной точки U. Однако, для частного случая ? = 1 существуют также и глобальные инварианты Римана. Пусть (,) и (,) пары левых и правых собственных векторов матрицы B. Введем в рассмотрение матрицы


, , ,


где, не нарушая общности, LR = I. Тогда матрица B может быть представлена в виде

= L?R,


а система (2.1) эквивалентно переписана так:


.


Инвариантами Римана (для нашей системы) являются, как известно, величины и , удовлетворяющие соотношениям


,


Полагая , приводим систему (2.1) к диагональному виду



Нахождение инвариантов Римана и подробно обсуждается в [2]. Мы же здесь приведем только результат (для случая ? = 1):


(2.5)


Заметим, что знание инвариантов Римана важно, например, для правильной постановки граничных условий при z = 0 и z = L. Так помимо того, что в силу (2.3) на каждой границе для системы (2.1) должно быть поставлено по одному граничному условию, необходимо, чтобы эти граничные условия разрешались для "уходящих инвариантов Римана, т. е для инвариантов Римана, соответствующих уходящим характеристикам. Вопрос о граничных условиях обсуждается в следующем параграфе.


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


Рассмотрим вначале случай одной артерии. Как уже отмечалось, для артерии длины L на левой (при z = 0) и на правой (при z = L) границах необходимо поставить по одному граничному условию. С математической точки зрения абсолютно неважно какими будут эти условия, а важно лишь то, чтобы они могли быть приведены к виду, когда "уходящие инварианты Римана выражаются через "приходящие. Так как в нашем случае в силу (2.3)"уходящим инвариантом Римана на левой границе является , а на правой границе - , необходимо, чтобы граничные условия приводились к виду


, , (3.1)


где и - некоторые функции.

С другой стороны, граничные условия должны быть разумными с физической точки зрения.

На левой границе таким разумным граничным условием является задание давления:


, (3.2)


где - некоторый входной профиль давления, (в нормальном случае мы рассматривали кусок синусоиды) моделирующий давление на выходе из сердца (Рис.3.1.). С учётом (1.8), (2.5) можно проверить, что граничное условие (3.2) является хорошим и с математической точки зрения, поскольку оно может быть приведено к виду первого условия из (3.1).


Рис.3.1 График давление на выходе из сердца.


Что касается правой границы, то мы поставили там условие постоянства давлений,


=9.32*104Дин/см2=70 мм рт ст, (3.3)


которое также является хорошим с математической точки зрения.

С физической точки зрения, является вполне логичным задавать

в качестве правого граничного условия - условие постоянства площади осевого сечения A, на концах последних артерий

.


с физической точки зрения, является вполне логичным задавать

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

Кроме этого, рассматривался вариант задания, в качестве левого граничного условия, в виде зависимости потока от времени


,


но оказалось, что это некорректно ввиду условия 3.1.

Планируется в дальнейшем для правых граничных условий записать закон Дарси (закон фильтрации крови в тканях) и ввести его в модель.

Граничные условия в точках ветвления артерий.


Рис.3.2 Одномерная модель бифуркации артерии (метод разбиения области).


Артериальная система характеризуется наличием ветвлений. Течение крови в местах ветвления артерий является существенно трехмерным. Однако, используя метод разбиения области, можно описать его в рамках одномерной модели. На Рис.3.2 изображена схема одномерной бифуркации. Мы упрощаем реальную геометрическую структуру бифуркации, полагая, что ветвление происходит в одной точке, и пренебрегая эффектами, связанными с углами ветвления.

Для того, чтобы решить три системы уравнений в областях (основная ветвь), и соответственно, нам снова необходимо найти интерфейсные условия в точке z = ? (= ?±). Так как системы являются гиперболическими и имеют по одной уходящей характеристике на границе {z = ?, }, требуется поставить три интерфейсных условия. Для уравнения неразрывности (первое в системе (1.7)) получаем непрерывности объемного кровотока в точке z = ? (условие сохранения массы крови при переходе через точку бифуркации), т.е.


при z = ?. (3.4).


Что касается интерфейсного условия для уравнения движения (второго в системе (1.7)), то в литературе часто используется условие непрерывности давления p.

В работах Formaggia и Quarteroni [8] было показано, что для случая ? = 1 интерфейсным условием, которое для задачи в подобластях , и гарантирует те же условия "численной" устойчивости, что и для исходной задачи в области D, является условие непрерывности на разрыве полного давления


. (3.5)


Таким образом, второе условие - непрерывность полного давления, в точке ветвления:


при z=Г. (3.6)


В итоге мы получаем три требуемых интерфейсных условия (3.4) и (3.6).

Аналогичным образом можно получить одномерную модель гемодинамики для всей артериальной системы человека, состоящей из 55 крупных артерий тела человека (Рис 3.3), список которых был использован в работе Lamponi (Табл.3.1) [1].


Таблица 3.1 Список и характеристики 55 крупных артерий тела человека.

#Название артерииДлина (см) Площадь сечения (см2) Beta (?) (кг/с2) Коэффициент отражения (Rt) 1Ascending Aorta4.05.983388-2Aortic Arch I2.05.147348-3Brachiocephalic3.41.219932-4R. Subclavian I3.40.5621692-5R. Carotid17.70.4322064-6R. Vertebral14.80.123103600.3027R. Subclavian II42.20.5101864-8R. radial23.50.106114640.2739R. Ulnar I6.70.1458984-10R. Interosseous7.90.031515760.31911R. Ulnar II17.10.13397840.29812R. internal Carotid17.60.121105760.26113R. external Carotid17.70.12198680.2614Aortic Arch II3.93.142520-15L. Carotid20.80.4302076-16L. internal Carotid17.60.121105760.26117L. external Carotid17.70.12198680.26418Thoracic Aorta I5.23.142496-19L. Subclavian I3.40.5621664-20Vertebral14.80.123103600.30221L. Subclavian II42.20.5101864-22L. Radial23.50.106114640.27423L. Ulnar I6.70.1458984-24L. Interosseous7.90.031515760.31925L. Ulnar II17.10.13397840.29826Intercostals8.00.19635400.20927Thoracic Aorta II10.43.017468-28Abdominal I5.31.911668-29Celiac I2.00.4781900-30Celiac II1.00.1267220-31Hepatic6.60.15245680.30832Gastric7.10.10262680.30733Splenic6.30.23832240.3134Superior Mesenteric5.90.43022760.31135Abdominal II1.01.247908-36L. Renal3.20.33222640.28737Abdominal III1.01.0211112-38R. Renal3.20.15947240.28739Abdominal IV10.60.6971524-40Inferior Mesenteric5.00.08075800.30641Abdominal V1.00.5781596-42R.common Iliac5.90.3282596-43L.common Iliac5.80.3282596-44L. external iliac14.40.2525972-45L. internal Iliac5.00.181125360.30846L. Femoral44.30.13910236-47L. deep Femoral12.60.126106080.29548L. posterior Tibial32.10.110232320.24149L. anterior Tibial34.30.060369720.23950R. external Iliac14.50.2525972-51R. internal Iliac5.10.181125360.30852R. Femoral44.40.13910236-53R. deep Femoral12.70.126106080.29654L. posterior Tibial32.30.110232320.24155R. anterior Tibial34.40.060369720.239

Рис.3.3 Схема дерева из 55 крупных артерий человека.


Для этого на входе первой артерии (аорты) ставится граничное условие вида (3.2), моделирующее работу сердца, на выходе конечных артерий граничные условия вида (3.3), а в точках ветвления артерий мы требуем выполнения интерфейсных условий (3.4) и (3.6). При этом в каждой отдельной артерии течение крови описывается системой гемодинамики (1.7). Считая, что исходный радиус артерий не меняется с изменением их длины, мы ставим начальные данные для площади осевого сечения A для каждой артерии:


для (3.7)


где - площади сечения отдельных артерий взяты из таблицы параметров артериальной системы человека (Табл.3.1).

Что касается начального условия для объемного кровотока Q, то на практике в численных расчётах артериальной системы обычно используют условие .

Значения таких физических параметров как длина сосудов и коэффициент также были взяты из работы Lamponi [1].


1.4 Реализация модели


1.4.1 Сведение к обыкновенному дифференциальному уравнению

Рассмотрим уравнение (2.1), к которому свелась система гемодинамики


(4.1)


где вектор неизвестных, и при ? = 1


, .


Требуется найти распределение U=U (t,z) вдоль оси z в любой момент времени t, при заданных начальных и граничных условиях.

Дискретизация по времени. Рассмотрим пространственно-временную систему координат {t,z}. В полуполосе построим по оси t полосы с шагом по времени


k=0,1,… Пусть Тогда


И вместо (4.1) получаем


(4.2)


Перепишем последнее выражение в виде


(4.3)


Таким образом, у нас получилось дифференциальное уравнение


(4.4)


где матрица G имеет вид


(4.5)


а вектор f есть


. (4.6)


Здесь в элементы G и f входят известные величины , которые берутся из предыдущего временного слоя.


1.4.2 Линеаризация граничных условий

Граничные условия (3.1), (3.3), (3.6) не являются линейными. Для линеаризации произведём следующее преобразование


(4.7)


Для линеаризации (3.6) также используются ис предыдущего шага по времени:


(4.8)


Обозначим


, (4.9)

, (4.10)

, (4.11)


тогда


. (4.12)


После этого, граничные условия (3.2), (3.3) для обыкновенного дифференциального уравнения (4.4), с учётом (4.9) - (4.12) формулируются в виде


, (4.13)

. (4.14)


Граничное условие (3.6) (в общем случае ) будет выглядеть так


. (4.15)


Решение. Для численного решения дифференциального уравнения (4.4) с граничными условиями (4.13), (4.14), (3.3), (4.15) используем метод ортогональной прогонки.


.4.3 Метод ортогональной прогонки

Метод ортогональной прогонки предназначен для численного решения краевой задачи


, , , (4.16)


где А - матрица размера nЧn, f - вектор-функция размера n, L,R - прямоугольные матрицы размера kЧn и pЧn, их ранги соответственно равны k и p. Элементы матрицы A (x) и вектор-функции предполагаются непрерывными на отрезке. Метод основан на представлении решения задачи (4.16) через решения серий задач Коши. Метод ортогональной прогонки зарекомендовал себя на практике как эффективное и надежное средство численного решения краевой задачи. Он обладает повышенной устойчивостью к погрешностям округлений при реализации метода на компьютере.

Подробное описание алгоритма метода ортогональной прогонки представлено в работах [9] и [10].


1.4.4 Теория возмущений (замечание о точности решения)

Пусть краевая задача


, , , (4.17)


правильно поставлена. (Здесь A (x) - матрица размера n Ч n, L, R - прямоугольные матрицы размера k Ч n и p Ч n соответственно, причем k + p =n, вектор-функция f (x) - размера n). Это означает, что краевая задача (4.17) имеет единственное решение для любых , ? и любой непрерывной на отрезке [0, d] вектор-функции f (x). Тогда решение задачи представляется c помощью матриц Грина G (x, s), , (см. [11])


, (4.18)


где - матрица-функция размера n Ч k, решение задачи



матрица-функция размера n Ч p, решение задачи


где - нулевая матрица размера l Ч m.

Будем говорить, что задача (4.17) хорошо обусловлена, если для любых x и s из интервала верны оценки


(4.19)


Из (4.18) следует, что для решения хорошо обусловленной задачи (4.17) справедлива оценка


. (4.20)


Это позволяет утверждать, что решение задачи (4.17) устойчиво по отношению к возмущениям. Предположим, что мы одновременно рассматриваем две "близкие" краевые задачи. Первая задача состоит в отыскании на отрезке 0 ? x ? d решения u (x) системы


,


удовлетворяющего граничным условиям


, .


Во второй задаче требуется на том же отрезке найти решение системы


(4.21)


для которого


,

.


Матрицы , , входящие в граничные условия, предполагаются имеющими то же число строк и столбцов, что и L, R соответственно. Утверждение, что первая и вторая задачи "близкие, можно понимать как утверждение, что имеют место неравенства


(4.22)


где


- достаточно маленькое число.


Всюду f (x), - непрерывные вектор-функции. Оказывается, что из разрешимости первой задачи при не слишком большом ? вытекает разрешимость второй.

Теорема 1. Пусть u (x) - решение правильно поставленной задачи (4.17), для функций Грина, которой справедливы оценки (4.19).

Пусть наряду с краевой задачей (4.17) имеется близкая к ней краевая задача (4.21), причем для элементов, определяющих эту краевую задачу, справедливы оценки (4.22).

Тогда при


? < min{1/ [2K (2 + d)], 1/ (2 + d) } (4.23)

существует и единственно - решение задачи (4.21) и для него верна оценка


, (4.24)


причем µ - число обусловленности задачи равно


µ = K (2 + d) (1 + 2KF),


где


.


(Доказательство - в [10].)

Оценка означает, что решение u (x) краевой задачи непрерывно зависит от коэффициентов A, L, R и правых частей , ?, f (x). В оценку непрерывности входит длина d отрезка, на котором ищется решение, и оценка K норм матриц Грина.


1.4.5 Генерация матриц граничных условий для метода ортогональной прогонки для системы из 55 артерий

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

=l, Ru=r


В случае 55 артерий L, R - матрицы размера (55Ч110), вектор неизвестных



(110*1), вектора и соответственно (55*1).

Артериальное дерево укладывается определённым образом как показано на рисунке 4.1


Рис.4.1 Укладка артериального дерева.


Рассмотрим, как выписываются матрицы L, R и вектора и .

Рассмотрим k-ый шаг по времени.

Граничное условие для первой артерии выражается формулой (4.13), таким образом первая строчка матрицы L будет выглядеть так:



Первый элемент вектора


.


Граничные условия для всех конечных артерий выражаются формулой (4.14), и если m - конечная артерия, то соответствующая ей строчка в матрице будет



И соответствующий ей элемент вектора


.


Граничные условия для бифуркации артерий, выражаются формулами (3.3), (4.15) или в общем случае, когда i артерия разветвляется на j и l



Соответствующие этим условиям строчки матрицы L


,


соответствующие элементы вектора


.

Аналогичным образом составляется матрица R и вектор .


1.4.6 Описание работы программы

Для работы с данной моделью был создан программный код на языке Java для интеграции ее с базой данных BMOND (#"15" src="doc_zip185.jpg" /> (tDelta), число разбиений одного сосуда по длине (Number of segmentation) и число разбиений сосуда для интегрирования (Segmentation for integrating), задаваемые в закладке Parameters (Рис.4.2.). В закладке Table of vessels задаётся длина сосуда, площадь его сечения (в начальный момент времени) и коэффициент (Рис.4.1.).

Выходными данными программы являются два массива матриц - A [i], Q [i], где сохранено значение площади сечения сосудов и потока в каждой точке разбиения по времени и по пространству, которые используются BioUML для построения графиков зависимости давления в заданных точках заданных сосудов от времени.


Рис.4.1 Пользовательский интерфейс пакета программ BioUML для плагина "Hemodynamics" и типа диаграмм "Arterial tree". Активна закладка Table of vessels.


Рис.4.2 Пользовательский интерфейс пакета программ BioUML для плагина "Hemodynamics" и типа диаграмм "Arterial tree". Активна закладка Parameters.

1.4.7 Описание структуры программы

В системе BioUML любая биологическая система описывается в виде графа, где с каждой вершиной или ребром графа связан набор определенных свойств. В модели гемодинамики вершины и ребра графа могут быть следующих типов

вершины: сердце, точка разветвления сосудов, контрольная точка (в которой можно наблюдать изменение давления);

рёбра - сосуды.

Свойства сосудов: длина, радиус, коэффициент .

Свойства точки ветвления: номер сосуда, который разветвляется; номера двух дочерних сосудов.

Система BioUML обеспечивает пользовательский интерфейс для графического представления редактирования структуры графа. Для удобства редактирования параметров вершин графа (параметры сосудов) в виде таблицы создана специальная вкладка.


Таблица 4.1 UML диаграмма, описывающая структуру классов в пакете "Hemodynamics.


Сплошные стрелки - включение одного объекта в другой. Пунктирные стрелки - реализация интерфейса. Красные - использование одним объектом другого.


Таблица 4.1 UML диаграмма, описывающая структуру классов в пакете "Hemodynamics. (Продолжение.)


Описание работы алгоритма

При создании графа сосудистого древа в BioUML, в программе автоматически создаётся объект класса ArterialTreeDiagramType, затем с помощью метода convert класса DiagramToArterialTreeConvertor из диаграммы генерируется объект (atm) класса ArterialBinaryTreeModel.

После того, как мы запускаем моделирование (нажимаем на кнопку "Start Model"), запускается метод solve из класса HemodynamicsModelSolver с параметром atm.

В методе solve реализована ортогональная прогонка, которая в том числе использует метод QR-разложения матрицы (QRDecomposition).

Для генерации матриц граничных условий L и R используется класс ArterialBinaryTreeMatrixResolver, в котором есть функция resolveMatrix, использующая метод traverseTree, реализующий интерфейс TreeTraversal. С Помощью TreeTraversal реализуется обход дерева.

Чтобы сгенерировать программный код для Matlab, в тесте создаётся объект класса ArterialBinaryTreeModel, для него запускается метод getMatlabPresentation класса ArterialBinaryTreeMatlabGenerator, после запуска теста программа автоматически генерирует m-файл, который в свою очередь является готовой программой на Matlabе.

Вторая часть таблицы 4.1 демонстрирует реализацию пользовательского интерфейса. Класс HemodynamicsViewPart содержит ArterialBinaryTreePane - панель "Hemodynamics в BioUML. На этой панели есть две вкладки Parameters" и Table of Vessels, которые реализованы в классах ParameterListTab и TableOfVesselsTab соответственно.


1.5 Тестовые расчёты


После реализации модели в BioUML был проведен ряд экспериментов in silico. Для всех 55 артерий человека были получены графики зависимости давления в них от времени и подтверждение их соответствия диапазону реальных физиологических значений параметров in vivo в норме и при патологии.


.5.1 Моделирование гемодинамики в норме

Некоторые результаты расчёта приведены ниже в виде графиков давления, посчитанного для срединных точек каждого из сосудов (Рис.5.1 - 5.4). Используемый математический и программный аппарат фактически позволяет отслеживать давление в любом сосуде в любой его точке в любой момент времени. Полученные графики показывают, что осцилляции давления спадают по мере удаления от сердца, что также хорошо согласуется с экспериментальными данными (Рис.5.3, 5.4).


Рис.5.1 Графики давления крови в восходящей ветви аорты (красным), дуге аорты (синим) и плечеголовной артерии (зеленым).


Рис.5.2 Графики давления крови в восходящей ветви аорты (желтым), правой сонной артерии (красным) и ее внешней (зеленым) и внутренней (синим) ветвях.


Рис.5.3 Графики давления крови в восходящей ветви аорты (зеленым), в левой (красным) и правой (синим) почечных артериях.

Рис.5.4 Графики давления крови в восходящей ветви аорты (желтым), бедренной (красным) левой задней голенной (синим) и левой передней голенной (зеленым) артериях.


1.5.2 Моделирование физической нагрузки и тахикардии

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

Входной профиль артериального давления был взят из промежуточных результатов дипломной работы Семисалова Б. "Построение, численное моделирование и анализ комплексной модели регуляции артериального давления, включая биофизические и биохимические блоки" по исследованию модели Солодянникова (Рис.5.5).


Рис.5.5 Динамика артериального давления на выходе из сердца при моделировании физической нагрузки на сердечно-сосудистую систему.


Далее приведены графики давления в некоторых артериях в зависимости от времени (Рис 5.6-5.9), которые были получены с помощью одномерной модели, в которой на вход был подан профиль давления, приведённый на Рис.5.5.


Рис.5.6 Графики давления крови в восходящей ветви аорты (красным), дуге аорты (синим) и плечеголовной артерии (зеленым) при физической нагрузке (начинается с 30 сек.).


Рис.5.7 Графики давления крови в восходящей ветви аорты (желтым), правой сонной артерии (красным) и ее внешней (зеленым) и внутренней (синим) ветвях.


Рис.5.8 Графики давления крови в восходящей ветви аорты (зеленым), в левой (красным) и правой (синим) почечных артериях при физической нагрузке (начинается с 30 сек.).


Рис.5.9 Графики давления крови в восходящей ветви аорты (желтым), бедренной (красным) левой задней голенной (синим) и левой передней голенной (зеленым) артериях при физической нагрузке (начинается с 30 сек.).


Из представленных графиков хорошо видно влияние физической нагрузки не только на динамику давления крови в восходящей ветви аорты, но также и в удалённых от сердца артериях.

Затем были проведены эксперименты по моделированию тахикардии. Входной профиль давления был взят из результата модели Солодянникова, и подставлен на вход в одномерную модель гемодинамики. Результаты эксперимента приведены в дипломной работе Семисалова Б.

1.5.3 Применение метода итераций для оптимизации тестовых расчётов

В начальный момент времени мы задаём нереалистичные условия (Q (0) =0, A (z) =const, , d - длина сосуда). Поэтому выход на действительное решение происходит после определённого промежутка времени. Если применить метод добавления итераций, то выход происходит быстрее.

Рассмотрим уравнение (4.1)


(5.1)

,


где . Тогда уравнение (5.1) будет следующим


. (5.2)


В пункте 4.1 мы использовали уравнение


, (5.3)


которое является приближением уравнения (5.2).

В модель вместо уравнения (5.3) можно ввести метод итераций

1 шаг: вычисляем также, как и в пункте (4.1) по формуле (5.3)


. (5.4)


шаг: вычисляем по формуле:


. (5.5)

шаг: вычисляем по формуле:


. (5.6)


Для (5.5) и (5.6) получается дифференциальное уравнение, аналогичное (4.4)


, (5.7)


где матрица G имеет вид


(4.5)


а вектор f есть


. (4.6)


Рис.5.6 Графики давления в сосудах при шаге разбиения по времени =0.02, одна итерация.


Рис.5.7 Графики давления в сосудах при шаге разбиения по времени =0.02, две итерации.

Заключение


В ходе выполненной работы были получены результаты:

·Созданы компьютерные программы на Java и на MatLab, реализующие одномерную модель гемодинамики с помощью метода прямых и метода ортогональной прогонки, позволяющие проводить моделирование и исследование артериальной системы человека.

·На языке Java был написан удобный пользовательский интерфейс, позволяющий использовать программу как математикам, так и биологам, не являющимся специалистами в области вычислительных методов.

·Данная программа в виде плагина интегрирована в пакет программ Biouml, созданный в Институте Системной Биологии, предназначенный для формального описания биологических систем, создания моделей и их изучения и способный работать с различными базами данных.

·Проведены эксперименты с созданной моделью, в ходе которых рассматривались различные состояния сердечно-сосудистой системы:

oнормальное состояние в покое;

oпри физической нагрузке;

oв случае патологического состояния - тахикардии.

Используемый математический и программный аппарат фактически позволяет отслеживать давление в любом сосуде в любой его точке в любой момент времени.

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

гемодинамика одномерная модель дифференциальная

Список литературы


[1] Lamponi D. One dimensional and multiscale models for blood flow circulation. // Lausanne, EPFL, 2004.

[2] Formaggia L., Veneziani A. Reduced and multiscale models for the human cardiovascular system. // Von Karman, 2003. Lecture notes of the 7th VKI lecture series, "Biological Fluid Dynamics

[3] Sherwin S., Franke V., Peiro J., Parker K. One-dimensional modeling of a vascular network in space-time variables. // Journal of Engineering Mathematics, Volume 47, Numbers 3-4, December 2003, pp.217-250 (34).

[4] Педли Т. Гидродинамика крупных кровеносных сосудов. // М: Мир, 1983.

[5] Павловский Ю.Н., Регирер С.А., Скобелева И.М. Гидродинамика крови // Итоги науки. М.: ВИНИТИ, 1970. Сер. Гидромеханика, 1968.

[6] Quarteroni A. Modeling the Cardiovascular System - A Mathematical Adventure: Part I.

// Siam News, Volume 34, Number 5.2001

[7] Абакумов М.В., Гаврилюк К.В., Есикова Н.Б., Кошелев В.Б., Лукшин А.В., Мухин С.И., Соснин Н.В., Тишкин В.Ф., Фаворский А.П. Математическая модель гемодинамики сердечно-сосудистой системы. // Дифференциальные уравнения. 1997, 33 (7), с.892-898.

[8] Quarteroni A., Formaggia L. Mathematical modelling and numerical simulation of the cardiovascular system. // In: Ciarlet P. G. (ed.). Handbook of numerical analyis, V.12, special volume: Ayache N. (ed.).computational Models for the Human Body. Amsterdam, Elsevier Science & Technology, 2003,P.3-129.

[9] Бибердорф Э.А., Попова Н.И. Глава 3: Формальное описание и моделирование сердечно-сосудистой системы. // А.М. Блохин, Л.Н. Иванова и др. Система кровообращения и артериальная гипертония: биофизические и генетико-физиологические механизмы, математическое и компьютерное моделирование. (Монография). Издательство СОРАН, Новосибирск, 2008. (В печати).

[10] Кузнецов С.В. Развитие метода ортогональной прогонки.

Диссертация. // Институт математики СО РАН, Новосибирск: 1988.

[11] Неймарк М.А. Линейные дифференциальные операторы. // М: Наука, 1969.

Приложение


. Пример сгенерированного Matlab кода

function [Energy] =hemodynamicModelSimulation ()=2; %***********number of unknown functions*****************=5; %*************segmentation for integrating*******************

%******************physical parameters**************************=1;_=0.05;=0.035;=3*10^6;=10; %***********number of segmentation************************=0.05; %*********time step***********************************=0; %***********start time************************************=50; %***********end time************************************=55; %********number of vessels*****************************

%******************vessels length******************************(1) =4.0;(2) =2.0;

length (54) =32.3;(55) =34.4;

%******************vessels initial area***************************(1) =5.983;(2) =5.147;

A0 (54) =0.11;(55) =0.06;

%*************sign*******************************************= [1; - 1; 1; - 1; 1; 1; - 1; 1; 1; - 1; - 1; 1; 1; - 1; - 1; 1; - 1; 1; - 1; - 1; 1; - 1; 1; - 1; 1; - 1; - 1; 1; - 1; - 1; 1; - 1; - 1; 1; - 1; 1; 1; - 1; - 1; 1; 1; 1; - 1; - 1; - 1; 1; - 1; - 1; 1; 1; - 1; - 1; 1; - 1; - 1];

%******************derived physical parameters*******************(1) =388.0;(2) =348.0;

gamma (54) =23232.0;(55) =36972.0;=8*pi*nu;

%******************segmentation*******************************j=1: Nvet,(j) =length (j) /N;;

%******************dimension*********************************_a=Nfun*Nvet;

%******************solution***********************************=zeros ([Nfun, (N+1),Nvet]);i=1: N+1,for j=1: Nvet,(1, i,j) =A0 (j);;;

%*************time iteration***********************************=round (T/tau);=zeros (Kmax);pr=1: 4,for k=1: Kmax,=t+tau;(k) =t;

%*************marching **************************************

%*************left boundary condition***************************j=1: Nvet,(j) =V (1,1,j);(j) =V (2,1,j);(j) =gamma (j) / (A0 (j) * (sqrt (An (j)) +sqrt (A0 (j))));(j) = (rho*Qn (j)) / (2*An (j) ^2);(j) =gamma (j) / (sqrt (An (j)) +sqrt (A0 (j)));;=t-0.86*fix (t/0.86);pr==1sh<0.37 c (1) =c (1) +0.001* (max ( (-679.328*sh^2+502.703*sh+7),70));(1) =c (1) +0.001* (max ( (387.339*sh^2-666.222*sh+293.476),70));; end;pr==2sh<0.37 c (1) =c (1) +0.01* (max ( (-679.328*sh^2+502.703*sh+7),70));(1) =c (1) +0.01* (max ( (387.339*sh^2-666.222*sh+293.476),70));; end;pr==3sh<0.37 c (1) =c (1) +0.1* (max ( (-679.328*sh^2+502.703*sh+7),70));(1) =c (1) +0.1* (max ( (387.339*sh^2-666.222*sh+293.476),70));; end;pr==4sh<0.37 c (1) =c (1) + (max ( (-679.328*sh^2+502.703*sh+7),70));(1) =c (1) + (max ( (387.339*sh^2-666.222*sh+293.476),70));; end;= [a (1),b (1),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;

,0,a (2),b (2),0,0,-a (4),-b (4),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,a (9),b (9),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];= [c (1); c (2) - c (4); c (2) - c (5); 0; c (18) - c (26); c (18) - c (27); 0; c (28) - c (34); c (28) - c (35); 0; c (36) +9; c (37) - c (38); c (37) - c (39); 0; c (40) +9; c (41) - c (42); c (41) - c (43); 0; c (44) - c (46); c (44) - c (47); 0; c (48) +9; c (49) +9; c (45) +9; c (50) - c (52); c (50) - c (53); 0; c (54) +9; c (55) +9; c (51) +9; c (29) - c (30); c (29) - c (31); 0; c (32) +9; c (33) +9; c (19) - c (20); c (19) - c (21); 0; c (22) +9; c (23) - c (24); c (23) - c (25); 0; c (16) +9; c (17) +9; c (3) - c (6); c (3) - c (7); 0; c (10) +9; c (11) - c (12); c (11) - c (13); 0; c (14) +9; c (15) +9; c (8) +9; c (9) +9];

[n_l,n] =size (L);

%*************right boundary condition**************************j=1: Nvet,(j) =V (1,N+1,j);(j) =V (2,N+1,j);(j) =gamma (j) / (A0 (j) * (sqrt (An (j)) +sqrt (A0 (j))));(j) = (rho*Qn (j)) / (2*An (j) ^2);(j) =gamma (j) / (sqrt (An (j)) +sqrt (A0 (j)));;= [a (1),b (1),-a (2),-b (2),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;(1),b (1),0,0,-a (3),-b (3),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;

0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];= [c (1) - c (2); c (1) - c (3); 0; c (4) - c (18); c (4) - c (19); 0; c (26) +9; c (27) - c (28); c (27) - c (29); 0; c (34) +9; c (35) - c (36); c (35) - c (37); 0; c (38) +9; c (39) - c (40); c (39) - c (41); 0; c (42) - c (44); c (42) - c (45); 0; c (46) - c (48); c (46) - c (49); 0; c (47) +9; c (43) - c (50); c (43) - c (51); 0; c (52) - c (54); c (52) - c (55); 0; c (53) +9; c (30) - c (32); c (30) - c (33); 0; c (31) +9; c (20) +9; c (21) - c (22); c (21) - c (23); 0; c (24) +9; c (25) +9; c (5) - c (16); c (5) - c (17); 0; c (6) - c (10); c (6) - c (11); 0; c (12) +9; c (13) - c (14); c (13) - c (15); 0; c (7) - c (8); c (7) - c (9); 0];

%*************assignment ************************************_od=n-n_l; % number of linealy independent solutions homogeneous left boundary condition for one vessel_0=zeros ([n,n_od,N+1]);_f=zeros ([n,N+1]);_0=zeros ([n,n_od,N+1]);_f=zeros ([n,N+1]);=zeros ([n_od+1,n_od+1,N]); %orthogonalising matrices=zeros ([Nfun,N+1,Nvet]); %solution=zeros ([n_od+1,N+1]); %coefficients of return marching

%*************processing of left boundary condition****************

[Q,R_left] =qr (L');_trans=Q';= (R_left (1: n_l,1: n_l)) ';_f (1: n_l,1) =inv (RR) *l;_f (:,1) =Q*Y_f (:,1);_0 (:,:,1) =Q (:,n_l+1: n);_0 (:,:,1) =Y_0 (:,:,1);_f (:,1) =Y_f (:,1);

%**************straight marching********************************j=1: Nvet;=V (2,1,j);=V (1,1,j);= (gamma (j) *sqrt (a)) / (2*rho*A0 (j));(:,:,j) =-znak (j) /tau* [2*q/ (q^2/a-c*a), 1/ (- (q/a) ^2+c);

,0];(:,j) =znak (j) * [q* (1/tau+Kr/a) / ( (q/a) ^2-c);/tau];;

%**************dimensional segmentation loop*********************i=1: N,

% vessels loopj=1: Nvet=V (2, i+1,j);=V (1, i+1,j);= (gamma (j) *sqrt (a)) / (2*rho*A0 (j));(:,:,j) =-znak (j) /tau* [2*q/ (q^2/a-c*a), 1/ (- (q/a) ^2+c);

,0];(:,j) =znak (j) * [q* (1/tau+Kr/a) / ( (q/a) ^2-c);/tau];

% solution of Cauchy problem=h (j) /M;_mat= (G2 (:,:,j) - G1 (:,:,j)) /M;_fun= (f2 (:,j) - f1 (:,j)) /M;_0=Z_0 ( (j-1) *Nfun+1: j*Nfun,:, i);_f=Z_f ( (j-1) *Nfun+1: j*Nfun, i);

% integrationjj=1: M,=G1 (:,:,j) + (jj-1/2) *del_mat; % intermediate matrix=f1 (:,j) + (jj-1/2) *del_fun; % intermediate right part_0=A*z_0*del;_0=z_0+del_0; %solution of homogeneous equations_f= (A*z_f+ff) *del;_f=z_f+del_f; %solution of nonhomogeneous equations;

% new value_0 ( (j-1) *Nfun+1: j*Nfun,:, i+1) =z_0;_f ( (j-1) *Nfun+1: j*Nfun, i+1) =z_f;;

% orthogonalizaiton= [Y_0 (:,:, i+1),Y_f (:, i+1)];

[Qpr,Rpr] =qr (Y);_f (:, i+1) =Qpr (:,n_od+1) *Rpr (n_od+1,n_od+1);(n_od+1,n_od+1) =1;_0 (:,:, i+1) =Qpr (:,1: n_od);(:,:, i) =Rpr (1: (n_od+1),1: (n_od+1));;

%***************inversion of right boundary condition**************=inv (R*Z_0 (:,:,N+1)) * (r-R*Z_f (:,N+1));

%***************return marching********************************(:,N+1) = [alpha; 1];(:,N+1) = [Z_0 (:,:,N+1),Z_f (:,N+1)] *beta (:,N+1);N>1,for i=1: N,=N+1-i;(:,m) =inv (Ort (:,:,m)) *beta (:,m+1);(:,m) = [Z_0 (:,:,m),Z_f (:,m)] *beta (:,m);;(:,1) =inv (Ort (:,:,1)) *beta (:,2);;(:,1) = [Y_0 (:,:,1),Y_f (:,1)] *beta (:,1);

%***************vessels loop***********************************j=1: Nvet,(:,:,j) =Z ( (j-1) *Nfun+1: j*Nfun,:); %solution;

%***************reserved operators for output result****************= [U (:,1,1)];i=1: Nvet-2,UUL= [UUL (1: 2*i); U (:,1, i+1)];;= [UUL (1: 2*Nvet-2); U (:,1,Nvet)];= [U (:,N+1,1)];i=1: Nvet-2,UUR= [UUR (1: 2*i); U (:,N+1, i+1)];;= [UUR (1: 2*Nvet-2); U (:,N+1,Nvet)];=L*UUL-l;=R*UUR-r;(k) =norm (WL);(k) =norm (WR);

%***************move to next step******************************=U;

% saving results

% for j=1: Nvet

% OA (k,:,j) =U (1,:,j);

% OQ (k,:,j) =U (2,:,j);

%end;(k) =gamma (5) * ( (sqrt (U (1,round (N/2),5)) - sqrt (A0 (5))) /A0 (5));;

%**********graphing******************************************i=1: Kmax(i) =tau*i;;

%F=figure;

%F1=figure;

%F2=figure;=figure;

% graphing middle of 2 vessel area and flux

%figure (F); hold;

%plot (mas,p (:,2),'red'); hold; pause;

%figure (F1); hold;

%plot (mas,p (:,3),'blue'); hold; pause;

% graphing middle of 5 vessel area and flux

%figure (F2); hold;

%plot (mas,p (:,4),'red'); hold; pause;(F3); % hold;(mas,p (:),'blue'); %hold; pause;;(max (Lbound));(max (Rbound)); pause;=1;

2. Пример сгенерированного Java кода

public class HemodynamicsModelSolver

implements TreeTraversal

{

private int [] sign;

protected ArterialBinaryTreeModel atm;

private double [] vesselAreas;

private Matrix [] flux;

public double [] vesselBeta;

private double [] vesselLegths; // static

double [] An, Qn, a, b, c;[] OA, OQ;

final int NumberOfUnknownFunction = 2;

int SegmentationForIntegrating = 10, NumberOfSegmentation = 20, rho = 1;

double nu = 0.035, gamma = 265868.077635827404094

double tDelta = 0.02, tMin = 0, tMax = 40;

public void solve (ArterialBinaryTreeModel atm, Matrix [] flux, double [] heartPressure)

{

this. atm = atm;

int NumberOfVessels = atm. size ();= new double [atm. size ()];= new double [atm. size ()];= new double [atm. size ()];= new int [atm. size ()];. traverseTree (this);

for (int i = 0; i < atm. size (); i++)

{[i] = vesselBeta [i] * 1000;

}

double [] An = new double [NumberOfVessels];

double [] Qn = new double [NumberOfVessels];

double [] a = new double [NumberOfVessels];

double [] b = new double [NumberOfVessels];

double [] c = new double [NumberOfVessels];[] Y_0 = new Matrix [NumberOfSegmentation + 1];[] Z_0 = new Matrix [NumberOfSegmentation + 1];[] Ort = new Matrix [NumberOfSegmentation];[] G1 = new Matrix [NumberOfVessels];[] G2 = new Matrix [NumberOfVessels];[] U = new Matrix [NumberOfVessels];[] U0 = new Matrix [NumberOfVessels];[] V = new Matrix [NumberOfVessels];= new Matrix [NumberOfVessels];= new Matrix [NumberOfVessels];

int i, j, k, jj;

int Kmax = (int) Math. round (tMax / tDelta);

double sh = 0.0;

for (i = 0; i < NumberOfSegmentation + 1; i++)

{_0 [i] = new Matrix (NumberOfVessels * 2, NumberOfVessels);_0 [i] = new Matrix (NumberOfVessels * 2, NumberOfVessels);

if (i < NumberOfSegmentation)

{[i] = new Matrix (NumberOfVessels + 1, NumberOfVessels + 1);

}

}

for (i = 0; i < NumberOfVessels; i++)

{[i] = new Matrix (2, NumberOfSegmentation + 1);[i] = new Matrix (2, NumberOfSegmentation + 1);[i] = new Matrix (2, NumberOfSegmentation + 1);[i] = new Matrix (Kmax + 1, NumberOfSegmentation + 1);[i] = new Matrix (Kmax + 1, NumberOfSegmentation + 1);[i] = new Matrix (2,2);[i] = new Matrix (2,2);

}L = new Matrix (NumberOfVessels, NumberOfVessels * 2);L1 = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);Right = new Matrix (NumberOfVessels, NumberOfVessels * 2);Q = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);Q1 = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);R_left = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);Rr = new Matrix (NumberOfVessels, NumberOfVessels);Y_f = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);Z_f = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);z_f = new Matrix (2, 1);z_0 = new Matrix (2, NumberOfVessels);Y = new Matrix (NumberOfVessels * 2, NumberOfVessels + 1);Y1 = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);Qpr = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);Rpr = new Matrix (NumberOfVessels + 1, NumberOfVessels + 1);f1 = new Matrix (2, NumberOfVessels);f2 = new Matrix (2, NumberOfVessels);ff = new Matrix (2, 1);del_mat = new Matrix (2,2);del_fun = new Matrix (2, 1);A_ = new Matrix (2,2);alpha = new Matrix (NumberOfVessels, 1);Beta = new Matrix (NumberOfVessels + 1, NumberOfSegmentation + 1);Z = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);Z_ = new Matrix (NumberOfVessels * 2, NumberOfVessels + 1);r = new Matrix (NumberOfVessels, 1);l = new Matrix (NumberOfVessels, 1);Zero = new Matrix (NumberOfVessels * 2, NumberOfVessels + 1);ZeroY_f = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);

for (i = 0; i < NumberOfSegmentation + 1; i++)

{

for (j = 0; j < NumberOfVessels; j++)

{[j]. set (0, i, vesselAreas [j]);

}

}

for (j = 0; j < NumberOfVessels; j++)

{[j]. setMatrix (0, 0, 0, NumberOfSegmentation, V [j]. getMatrix (0, 0, 0, NumberOfSegmentation));[j]. setMatrix (0, 0, 0, NumberOfSegmentation, V [j]. getMatrix (1, 1, 0, NumberOfSegmentation));

}

for (k = 0; k < Kmax; k++)

{+= tDelta;

// experiment repeat marching

for (int k2 = 0; k2 < 1; k2++)

{

for (j = 0; j < NumberOfVessels; j++)

{

if (k2 == 0)

{[j] = V [j]. get (0, 0);[j] = V [j]. get (1, 0);

}

else

{[j] = U0 [j]. get (0, 0);[j] = U0 [j]. get (1, 0);

}[j] = vesselBeta [j] / (vesselAreas [j] * (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j])));[j] = (rho * Qn [j]) / (2 * An [j] * An [j]);[j] = vesselBeta [j] / (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j]));

}

if (k < Kmax / 2)

{. out. println ("tMin " + tMin);

}[0] +=heartPressure [k] * (1 - Math. exp ( - (tMin) * (tMin) / 100));

double [] flux1 = new double [NumberOfVessels];

int col = flux [0]. getColumnDimension ();

int fluxLength = flux. length;

if (k < Kmax / 2)

{. out. println ("a " + a [0] + " b " + b [0] + " c " + c [0]);

}

// resolve left matrix {7}leftResolver = new MatrixResolver (atm, a, b, c, true, vesselAreas, flux1, tMin);leftMatrix = leftResolver. resolveMatrix ();= leftResolver. getMatrix (leftMatrix);= leftResolver. getVector (leftMatrix);

for (j = 0; j < NumberOfVessels; j++)

{

if (k2 == 0)

{[j] = V [j]. get (0, NumberOfSegmentation);[j] = V [j]. get (1, NumberOfSegmentation);

}

else

{[j] = U0 [j]. get (0, NumberOfSegmentation);[j] = U0 [j]. get (1, NumberOfSegmentation);

}[j] = vesselBeta [j] / (vesselAreas [j] * (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j])));[j] = (rho * Qn [j]) / (2 * An [j] * An [j]);[j] = vesselBeta [j] / (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j]));

}

// resolve right matrix {8}rightResolver = new MatrixResolver (atm, a, b, c, false, vesselAreas, flux1, tMin);rightMatrix = rightResolver. resolveMatrix ();= rightResolver. getMatrix (rightMatrix);= rightResolver. getVector (rightMatrix);

// experiment repeat marching was here. out. println ("k2 " + k2);

// null matrixes Y_f,Z_f,beta, Z_0, Y_0, Ort, U_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfSegmentation, ZeroY_f);_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfSegmentation, ZeroY_f);. setMatrix (0, NumberOfVessels, 0, NumberOfSegmentation, ZeroY_f. getMatrix (0, NumberOfVessels, 0,NumberOfSegmentation));

for (int i1 = 0; i1 < NumberOfSegmentation + 1; i1++)

{_0 [i1]. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Zero. getMatrix (0, NumberOfVessels * 2 - 1, 0,NumberOfVessels - 1));_0 [i1]. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Zero. getMatrix (0, NumberOfVessels * 2 - 1, 0,NumberOfVessels - 1));

if (i1 < NumberOfSegmentation)

{[i1]. setMatrix (0, NumberOfVessels, 0, NumberOfVessels, Zero

. getMatrix (0, NumberOfVessels, 0, NumberOfVessels));

}

}

for (int i2 = 0; i2 < NumberOfVessels; i2++)

{[i2]. setMatrix (0, 1, 0, NumberOfSegmentation, ZeroY_f. getMatrix (0, 1, 0, NumberOfSegmentation));

}

// end of null matrixes. setMatrix (0, 2 * NumberOfVessels - 1, 0, NumberOfVessels - 1, L. transpose ());= new QRDecomposition (L1). getQFull ();_left = new QRDecomposition (L. transpose ()). getR ();= (R_left. getMatrix (0, NumberOfVessels - 1, 0, NumberOfVessels - 1)). transpose ();_f. setMatrix (0, NumberOfVessels - 1, 0, 0, (Rr. inverse ()). times (l));_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, 0, Q. times (Y_f. getMatrix (0, NumberOfVessels * 2 - 1, 0, 0)));_0 [0]. setMatrix (0, 2 * NumberOfVessels - 1, 0, NumberOfVessels - 1, Q. getMatrix (0, NumberOfVessels * 2 - 1,NumberOfVessels, NumberOfVessels * 2 - 1));_0 [0] = Y_0 [0];_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, 0, Y_f. getMatrix (0, NumberOfVessels * 2 - 1, 0, 0));

for (j = 0; j < NumberOfVessels; j++) // straight marching

{

double q, a_, q_k, a_k;

if (k2 == 0)

{= V [j]. get (1, 0);_ = V [j]. get (0, 0);_k = 0;_k = 0;

}

else

{= U0 [j]. get (1, 0);_ = U0 [j]. get (0, 0);_k = V [j]. get (1, 0);_k = V [j]. get (0, 0);

}

double c_ = (vesselBeta [j] * Math. sqrt (a_)) / (2 * rho * vesselAreas [j]);[j]. set (0, 0, ( - sign [j] / tDelta) * 2 * q / (q * q / a_ - c_ * a_));[j]. set (0, 1, ( - sign [j] / tDelta) / ( - (q / a_) * (q / a_) + c_));[j]. set (1, 0, - sign [j] / tDelta);[j]. set (1, 1, 0);

if (k2 == 0)

{. set (0, j, sign [j] * (q * (1/tDelta + (8 * Math. PI * nu) / a_) / ( (q / a_) * (q / a_) - c_)));. set (1, j, sign [j] * a_ / tDelta);

}

else

{. set (0, j, sign [j]

* ( (1/tDelta) * (q_k + (8 * Math. PI * nu) * tDelta * q / a_) / ( (q / a_) * (q / a_) - c_)));. set (1, j, sign [j] * a_k / tDelta);

}

}

for (i = 0; i < NumberOfSegmentation; i++) // dimensional

// segmentation loop

{

for (j = 0; j < NumberOfVessels; j++) // vessel loop

{

double q, a_, q_k, a_k;

if (k2 == 0)

{= V [j]. get (1, i + 1);_ = V [j]. get (0, i + 1);_k = 0;_k = 0;

}

else

{= U0 [j]. get (1, i + 1);_ = U0 [j]. get (0, i + 1);_k = V [j]. get (1, i + 1);_k = V [j]. get (0, i + 1);

}

double c_ = (vesselBeta [j] * Math. sqrt (a_)) / (2 * rho * vesselAreas [j]);[j]. set (0, 0, ( - sign [j] / tDelta) * 2 * q / (q * q / a_ - c_ * a_));[j]. set (0, 1, ( - sign [j] / tDelta) / ( - (q / a_) * (q / a_) + c_));[j]. set (1, 0, - sign [j] / tDelta);[j]. set (1, 1, 0);

if (k2 == 0)

{. set (0, j, sign [j] * (q * (1/tDelta + (8 * Math. PI * nu) / a_) / ( (q / a_) * (q / a_) - c_)));. set (1, j, sign [j] * a_ / tDelta);

}

else

{. set (0, j,[j]

* ( (1/tDelta) * (q_k + (8 * Math. PI * nu) * tDelta * q / a_) / ( (q / a_)

* (q / a_) - c_)));. set (1, j, sign [j] * a_k / tDelta);

}

// solution of Cauchy problem

double del = vesselLegths [j] / (NumberOfSegmentation * SegmentationForIntegrating);_mat. setMatrix (0, 1, 0, 1, (G2 [j]. minus (G1 [j])). times (1.0/SegmentationForIntegrating));_fun. set (0, 0, (f2. get (0, j) - f1. get (0, j)) / SegmentationForIntegrating);_fun. set (1, 0, (f2. get (1, j) - f1. get (1, j)) / SegmentationForIntegrating);_0 = Z_0 [i]. getMatrix (j * 2, j * 2 + 1, 0, NumberOfVessels - 1);_f = Z_f. getMatrix (j * 2, j * 2 + 1, i, i);

for (jj = 0; jj < SegmentationForIntegrating; jj++)

{_. setMatrix (0, 1, 0, 1, G1 [j]. plus (del_mat. times (jj + 0.5))); // intermediate

// matrix. setMatrix (0, 1, 0, 0, (f1. getMatrix (0, 1, j, j)). plus (del_fun. times (jj + 0.5))); // intermediate

// right

// part_0. plusEquals (A_. times (z_0. times (del)));_f. plusEquals ( ( (A_. times (z_f)). plus (ff)). times (del));

}

// new value_0 [i + 1]. setMatrix (j * 2, j * 2 + 1, 0, NumberOfVessels - 1, z_0);_f. setMatrix (j * 2, j * 2 + 1, i + 1, i + 1, z_f);

for (int i3 = 0; i3 < NumberOfVessels; i3++)

{[i3] = G2 [i3]. copy ();= f2. copy ();

}

}

// orthogonalizaiton. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Y_0 [i + 1]);. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Y_f. getMatrix (0, NumberOfVessels * 2 - 1,i + 1, i + 1));. setMatrix (0, 2 * NumberOfVessels - 1, 0, NumberOfVessels, Y);= new QRDecomposition (Y1). getQFull ();= new QRDecomposition (Y). getR (); // k==1 +_f. setMatrix (0, NumberOfVessels * 2 - 1, i + 1, i + 1, Qpr. getMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels,). times (Rpr. get (NumberOfVessels, NumberOfVessels)));. set (NumberOfVessels, NumberOfVessels, 1);_0 [i + 1]. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Qpr. getMatrix (0, NumberOfVessels * 2 - 1,0, NumberOfVessels - 1));[i]. setMatrix (0, NumberOfVessels, 0, NumberOfVessels, Rpr. getMatrix (0, NumberOfVessels, 0, NumberOfVessels));

}. setMatrix (0, NumberOfVessels - 1, 0, 0, ( (Right. times (Z_0 [NumberOfSegmentation])). inverse ())

. times (r. minus (Right. times (Z_f. getMatrix (0, NumberOfVessels * 2 - 1, NumberOfSegmentation,)))));

// return marching. setMatrix (0, NumberOfVessels - 1, NumberOfSegmentation, NumberOfSegmentation, alpha);. set (NumberOfVessels, NumberOfSegmentation, 1);_. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Z_0 [NumberOfSegmentation]. getMatrix (0,NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1));_. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Z_f. getMatrix (0, NumberOfVessels * 2 - 1,NumberOfSegmentation, NumberOfSegmentation));. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfSegmentation, NumberOfSegmentation, Z_. times (Beta. getMatrix (0,NumberOfVessels, NumberOfSegmentation, NumberOfSegmentation)));

if (NumberOfSegmentation > 1)

{

for (i = 0; i < NumberOfSegmentation; i++)

{

int m = NumberOfSegmentation - i;. setMatrix (0, NumberOfVessels, m - 1, m - 1, (Ort [m - 1]. inverse ()). times (Beta. getMatrix (0,NumberOfVessels, m, m)));_. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Z_0 [m - 1]. getMatrix (0,NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1));_. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Z_f. getMatrix (0,NumberOfVessels * 2 - 1, m - 1, m - 1));. setMatrix (0, NumberOfVessels * 2 - 1, m - 1, m - 1, Z_. times (Beta. getMatrix (0, NumberOfVessels, m - 1,m - 1)));

}

}

else

{

. setMatrix (0, NumberOfVessels, 0, 0, (Ort [0]. inverse ()). times (Beta

. getMatrix (0, NumberOfVessels, 1, 1)));

}_. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Y_0 [0]. getMatrix (0, NumberOfVessels * 2 - 1, 0,NumberOfVessels - 1));_. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Y_f. getMatrix (0, NumberOfVessels * 2 - 1, 0,0));. setMatrix (0, NumberOfVessels * 2 - 1, 0, 0, Z_. times (Beta. getMatrix (0, NumberOfVessels, 0, 0)));

// vessels loop

for (j = 0; j < NumberOfVessels; j++)

{[j]. setMatrix (0, 1, 0, NumberOfSegmentation, Z. getMatrix (j * 2, j * 2 + 1, 0, NumberOfSegmentation)); // solution

}

for (i = 0; i < NumberOfVessels; i++)

{[i] = U [i]. copy ();

}

if ( (k2 == 0) && (tMin < 2))

{

break;

}

} // end of experiment repeat marching

for (i = 0; i < NumberOfVessels; i++)

{[i] = U [i]. copy ();

}

for (j = 0; j < NumberOfVessels; j++)

{[j]. setMatrix (k + 1, k + 1, 0, NumberOfSegmentation, U [j]. getMatrix (0, 0, 0, NumberOfSegmentation));[j]. setMatrix (k + 1, k + 1, 0, NumberOfSegmentation, U [j]. getMatrix (1, 1, 0, NumberOfSegmentation));

}

} // end of general loop

}

public void setTDelta (double value)

{= value;

}

public double getTDelta ()

{

return tDelta;

}

public void setTMin (double value)

{= value;

}

public void setTMax (double value)

{= value;

}

public void setNumberOfSegmentation (int value)

{= value;

}

public void setSegmentationForIntegrating (int value)

{= value;

}

public Matrix [] getArea ()

{

return OA;

}

public Matrix [] getFlux ()

{

return OQ;

}

public void visit (Vessel v)

{

int i = v. index;. out. println (atm. size ()); // for debug[i] = v. length;[i] = v. area;[i] = v. getBeta ();

if (atm. isEvenNode (v))[i] = 1;

else[i] = - 1;

if (v. left! = null)(v. left);

if (v. right! = null)(v. right);

}

public static class MatrixResolver

implements TreeTraversal

{

public Matrix getMatrix (Matrix m)

{resultMatrix = new Matrix (atm. size (), atm. size () * 2);= (m. getMatrix (0, 2 * atm. size () - 1, 0, atm. size () - 1)). transpose ();

return resultMatrix;

}

public Matrix getVector (Matrix v)

{resultVector = new Matrix (atm. size (), 1);= (v. getMatrix (2 * atm. size (), 2 * atm. size (), 0, atm. size () - 1)). transpose ();

return resultVector;

}

private double [] a;

private double [] b;

private double [] c;

private Matrix result = null;

private int index = 0;

private boolean left = true;

private double [] vesselAreas;

private int k;

private ArterialBinaryTreeModel atm;

private double [] flux;

private double t;

public MatrixResolver (ArterialBinaryTreeModel atm, double [] a, double [] b, double [] c, boolean left, double [] vesselAreas,

double [] flux, double t)

{

this. atm = atm;

this. a = a;

this. b = b;

this. c = c;

this. left = left;

this. vesselAreas = vesselAreas;

this. flux = flux;

this. t = t;= new Matrix (2 * atm. size () + 1, atm. size ());

}

public Matrix resolveMatrix ()

{

if (left)

{

double flux0 = flux [0];

double fluxFin = flux [0];singlePressureFirstCondition = resolveSinglePressureEndCondition (0, a, b, c, flux0, fluxFin, t);. setMatrix (0, 2 * atm. size (), index, index, singlePressureFirstCondition);++;

}. traverseTree (this);

return result;

}

public void visit (Vessel v)

{[v. index] = v. area;

if (atm. isEvenNode (v) == left)

return;

if (v. left! = null && v. right! = null)

{singlePressureConditionLeft = resolveSinglePressureCondition (v, v. left. index, a, b, c);. setMatrix (0, 2 * atm. size (), index, index, singlePressureConditionLeft);++;singlePressureConditionRight = resolveSinglePressureCondition (v, v. right. index, a, b, c);. setMatrix (0, 2 * atm. size (), index, index, singlePressureConditionRight);++;singleFluxCondition = resolveSingleFluxCondition (v. index, v. left. index, v. right. index, a, b, c);. setMatrix (0, 2 * atm. size (), index, index, singleFluxCondition);++;

}

else

{

double flux0 = flux [0];

double fluxFin = flux [v. index];singlePressureEndCondition = resolveSinglePressureEndCondition (v. index, a, b, c, flux0, fluxFin, t);. setMatrix (0, 2 * atm. size (), index, index, singlePressureEndCondition);++;

}

}

private Matrix resolveSinglePressureCondition (Vessel vessel, int childIndex, double [] a, double [] b, double [] c)

{result = new Matrix (2 * atm. size () + 1, 1);

int parentIndex = vessel. index;. set (2 * parentIndex, 0, a [parentIndex]);. set (2 * parentIndex + 1, 0, b [parentIndex]);. set (2 * childIndex, 0, - a [childIndex]);. set (2 * childIndex + 1, 0, - b [childIndex]);. set (2 * atm. size (), 0, c [parentIndex] - c [childIndex]);

return result;

}

private Matrix resolveSingleFluxCondition (int i, int j, int k, double [] a, double [] b, double [] c)

{result = new Matrix (2 * atm. size () + 1, 1);. set (2 * i + 1, 0, 1);. set (2 * j + 1, 0, - 1);. set (2 * k + 1, 0, - 1);

return result;

}

private Matrix resolveSinglePressureEndCondition (int vesselIndex, double [] a, double [] b, double [] c, double flux0,double fluxFin, double t)

{result = new Matrix (2 * atm. size () + 1, 1);. set (2 * vesselIndex, 0, a [vesselIndex]);. set (2 * vesselIndex + 1, 0, b [vesselIndex]);

if (vesselIndex == 0)

{. set (2 * atm. size (), 0, c [vesselIndex]);

}

else

{. set (2 * atm. size (), 0, c [vesselIndex] + 70 * (1 - Math. exp ( - (t * t) / 100)));

}

return result;

}

}

}


Содержание Введение 1. Основная часть 1.1 Вывод модели 1.2 Анализ математических свойств модели 1.3 Граничные условия 1.4 Реализация модел

Больше работ по теме:

КОНТАКТНЫЙ EMAIL: [email protected]

Скачать реферат © 2017 | Пользовательское соглашение

Скачать      Реферат

ПРОФЕССИОНАЛЬНАЯ ПОМОЩЬ СТУДЕНТАМ