#
|
Название артерии
|
Длина (см)
|
Площадь сечения (см2)
|
Beta (β) (кг/с2)
|
Коэффициент отражения (Rt)
|
1
|
Ascending Aorta
|
4.0
|
5.983
|
388
|
-
|
2
|
Aortic Arch I
|
2.0
|
5.147
|
348
|
-
|
3
|
Brachiocephalic
|
3.4
|
1.219
|
932
|
-
|
4
|
R. Subclavian I
|
3.4
|
0.562
|
1692
|
-
|
5
|
R. Carotid
|
17.7
|
0.432
|
2064
|
-
|
6
|
R. Vertebral
|
14.8
|
0.123
|
10360
|
0.302
|
7
|
R. Subclavian II
|
42.2
|
0.510
|
1864
|
-
|
8
|
R. radial
|
23.5
|
0.106
|
11464
|
0.273
|
9
|
R. Ulnar I
|
6.7
|
0.145
|
8984
|
-
|
10
|
R. Interosseous
|
7.9
|
0.031
|
51576
|
0.319
|
11
|
R. Ulnar II
|
17.1
|
0.133
|
9784
|
0.298
|
12
|
R. internal Carotid
|
17.6
|
0.121
|
10576
|
0.261
|
13
|
R. external Carotid
|
17.7
|
0.121
|
9868
|
0.26
|
14
|
Aortic Arch II
|
3.9
|
3.142
|
520
|
-
|
15
|
L. Carotid
|
20.8
|
0.430
|
2076
|
-
|
16
|
L. internal Carotid
|
17.6
|
0.121
|
10576
|
0.261
|
17
|
L. external Carotid
|
17.7
|
0.121
|
9868
|
0.264
|
18
|
Thoracic Aorta I
|
5.2
|
3.142
|
496
|
-
|
19
|
L. Subclavian I
|
3.4
|
0.562
|
1664
|
-
|
20
|
Vertebral
|
14.8
|
0.123
|
10360
|
0.302
|
21
|
L. Subclavian II
|
42.2
|
0.510
|
1864
|
-
|
22
|
L. Radial
|
23.5
|
0.106
|
11464
|
0.274
|
23
|
L. Ulnar I
|
6.7
|
0.145
|
8984
|
-
|
24
|
L. Interosseous
|
7.9
|
0.031
|
51576
|
0.319
|
25
|
L. Ulnar II
|
17.1
|
0.133
|
9784
|
0.298
|
26
|
Intercostals
|
8.0
|
0.196
|
3540
|
0.209
|
27
|
Thoracic Aorta II
|
10.4
|
3.017
|
468
|
-
|
28
|
Abdominal I
|
5.3
|
1.911
|
668
|
-
|
29
|
Celiac I
|
2.0
|
0.478
|
1900
|
-
|
30
|
Celiac II
|
1.0
|
0.126
|
7220
|
-
|
31
|
Hepatic
|
6.6
|
0.152
|
4568
|
0.308
|
32
|
Gastric
|
7.1
|
0.102
|
6268
|
0.307
|
33
|
Splenic
|
6.3
|
0.238
|
3224
|
0.31
|
Superior Mesenteric
|
5.9
|
0.430
|
2276
|
0.311
|
35
|
Abdominal II
|
1.0
|
1.247
|
908
|
-
|
36
|
L. Renal
|
3.2
|
0.332
|
2264
|
0.287
|
37
|
Abdominal III
|
1.0
|
1.021
|
1112
|
-
|
38
|
R. Renal
|
3.2
|
0.159
|
4724
|
0.287
|
39
|
Abdominal IV
|
10.6
|
0.697
|
1524
|
-
|
40
|
Inferior Mesenteric
|
5.0
|
0.080
|
7580
|
0.306
|
41
|
Abdominal V
|
1.0
|
0.578
|
1596
|
-
|
42
|
R.common Iliac
|
5.9
|
0.328
|
2596
|
-
|
43
|
L.common Iliac
|
5.8
|
0.328
|
2596
|
-
|
44
|
L. external iliac
|
14.4
|
0.252
|
5972
|
-
|
45
|
L. internal Iliac
|
5.0
|
0.181
|
12536
|
0.308
|
46
|
L. Femoral
|
44.3
|
0.139
|
10236
|
-
|
47
|
L. deep Femoral
|
12.6
|
0.126
|
10608
|
0.295
|
48
|
L. posterior Tibial
|
32.1
|
0.110
|
23232
|
0.241
|
49
|
L. anterior Tibial
|
34.3
|
0.060
|
36972
|
0.239
|
50
|
R. external Iliac
|
14.5
|
0.252
|
5972
|
-
|
51
|
R. internal Iliac
|
5.1
|
0.181
|
12536
|
0.308
|
52
|
R. Femoral
|
44.4
|
0.139
|
10236
|
-
|
53
|
R. deep Femoral
|
12.7
|
0.126
|
10608
|
0.296
|
54
|
L. posterior Tibial
|
32.3
|
0.110
|
23232
|
0.241
|
55
|
R. anterior Tibial
|
34.4
|
0.060
|
36972
|
0.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 (#"588128.files/image089.gif"> (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.
%******************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;
}
}