Благовещенск, Амурская область, Россия
УДК 53 Физика
УДК 61 Медицина. Охрана здоровья
Основной задачей данной работы была разработка и исследование алгоритма численного решения системы нелинейных уравнений, описывающих тепловлагообмен в дыхательных путях человека. Разработанный алгоритм проанализирован на сходимость и погрешность. Получены удовлетворительные результаты.
дыхательные пути, слизистая оболочка, тепловлагообмен, математическая модель, система уравнений
УДК 612.221.4:001.891.573(.001.5)
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТЕПЛОВЛАГООБМЕНА В ДЫХАТЕЛЬНЫХ ПУТЯХ ЧЕЛОВЕКА (исследование)
© 2018. Н.В. Ульянычев2, канд. физ.-мат. наук, В.Ф. Ульянычева1, канд. физ.-мат. наук
1ФГБОУ ВО «Амурский государственный университет»,
2 Дальневосточный научный центр физиологии и патологии дыхания, Благовещенск
Основной задачей данной работы была разработка и исследование алгоритма численного решения системы нелинейных уравнений, описывающих тепловлагообмен в дыхательных путях человека. Разработанный алгоритм проанализирован на сходимость и погрешность. Получены удовлетворительные результаты.
Ключевые слова: дыхательные пути, слизистая оболочка, тепловлагообмен, математическая модель, система уравнений.
A MATHEMATICAL MODEL OF HEAT AND MOISTURE EXCHANGE IN THE RESPIRATORY TRACT OF A HUMAN (RESEARCH)
N.V. Ulyanychev, PhD; V.F. Ulyanycheva, PhD
Amur State University, Blagoveshchensk
Far-Eastern Scientific Center of Physiology and Pathology of Respiration, Blagoveshchensk
The main objective of this work was to develop and study an algorithm for the numerical solution of the system of linear equations describing the heat and moisture exchange in the human respiratory tract. The developed algorithm is analyzed for convergence and error. Satisfactory results were obtained.
Key words: respiratory tract, mucous membrane, heat and moisture exchange, mathematical model, system of equations.
В предыдущей статье нами получена система нелинейных дифференциальных уравнений, описывающих тепловлагообмен в дыхательных путях. Данная работа посвящена численному решению этой системы уравнений.
Методы численного решения
Для получения распределения влажности и температуры в однородном участке следует решить систему уравнений (1). Эта система является нелинейной и не имеет элементарного решения, чтобы ее решить численно, необходимо применить один из численных методов.
Систему (1) можно упростить. Во-первых, подставить второе уравнение системы в первое и получить соотношение:
kT*(TB – TM(x))/dY = hC*(TM(x) – TA(x)) + kC*hfg*(CM(x) – CA(x)). (2)
Это уравнение устанавливает зависимость температуры слизистой TM(x) от температуры и влажности воздуха а также влажности слизистой.
Также можно подставить четвертое уравнение системы (1) в третье и получить:
СM(x) = 486,081*10-11,04*(100 + TM(x))7,51/(1+0,00366*TM(x))/1000. (3)
Из этого уравнения видно, что влажность слизистой СM(x) зависит только от температуры слизистой TM(x).
Уравнение (3) можно было бы подставить в (2), тогда температура слизистой в получившейся системе зависела бы только от температуры и влажности воздуха, но делать этого не будем, так как влажность слизистой СM(x) – одна из величин, значение которой необходимо найти как одно из составляющих решение системы. Получаем систему их четырех уравнений, численное решение которой даст искомые распределения температуры и влажности воздуха и температуры и влажности слизистой:
В этой системе первое уравнение устанавливает зависимость для температуры слизистой, второе уравнение – зависимость влажности слизистой от температуры слизистой. По третьему уравнению можно рассчитать влажность воздуха, по четвертому – температуру воздуха.
Третье и четвертое уравнения системы (4) (для температуры и влажности воздуха) имеют вид:
y' = f(x,y). (5)
Также известны начальные условия:
y(x0) = y0. (6)
Для третьего и четвертого уравнения системы начальными условиями будут соответственно температура и влажность поступающего в однородный участок воздуха. Задача (5), (6) имеет вид задачи Коши [34, 35].
Для численного решения уравнений такого вида существует много различных методов. Среди них выделяют одноступенчатые и многоступенчатые методы. Все методы, входящие в оба этих класса методов представляют собой итерационный процесс, на каждом шаге которого вычисляется с определенной погрешностью yi для значения xi из некоторого промежутка.
Одноступенчатые методы используют информацию только об очередной точке решения и не используют информацию о ранее найденных точках. Однако при использовании этих методов приходится многократно вычислять функцию f(x,y) и затрачивать на это много машинного времени. Кроме того, при использовании этих методов весьма трудно получить оценку ошибки ограничения.
Среди многоступенчатых методов [34-36] есть различные семейства, обладающие многими достоинствами и недостатками. Здесь будет использован метод "прогноза и коррекции" [35, 36]. Этот метод в отличие от одноступенчатых методов (и как все многоступенчатые) использует информацию о предыдущих найденных точках и в отличие от многих других многоступенчатых [35] позволяет получить ошибку ограничения в качестве побочного продукта вычислений. Однако с помощью этого метода нельзя начать решение уравнения, так как в нем информация о ранее найденных точках используется обязательно.
Как указано в источниках [35, 36] придется использовать комбинацию из одноступенчатого метода и метода прогноза и коррекции.
Общая стратегия решения выглядит так:
1. Начать решение уравнения с помощью одноступенчатого метода.
2. Для дальнейшего решения использовать метод прогноза и коррекции.
3. В методе прогноза и коррекции использовать возможность получения ошибки ограничения в качестве побочного продукта вычислений.
В качестве одноступенчатого метода будем использовать классический метод Рунге – Кутта [34-36] четвертого порядка, который описывается системой следующих пяти соотношений:
где y – вычисляемая функция;
k1, k2, k3, k4 – коэффициенты.
Функция f(x,y) соответствует правым частям третьего и четвертого уравнения системы (4).
Следует отметить, что четвертое уравнение системы (4) зависит от третьего. Поэтому процесс численного решения этих уравнений должен происходить совместно. Сначала должно быть найдено значение CA(x) для некоторого xi, из третьего уравнения, а затем уже TA(x) для xi из четвертого уравнения.
Кроме того, третье и четвертое уравнение системы (4) связаны с первым и вторым уравнениями. А это означает, что перед тем, как вычислить температуру и влажность воздуха в i – й точке следует предварительно вычислить влажность и температуру слизистой.
После того, как будет с помощью системы (7) найдено необходимое количество точек, вычисления будут продолжены с помощью методов прогноза и коррекции.
Будет использован метод прогноза и коррекции второго порядка.
Вначале находится прогнозное значение по формуле [36]:
y(0)m+1 = ym-1 + 2h*f(xm,ym); (8)
где m – номер шага;
y(0)m+1 – прогнозное значение искомой функции в точке m+1;
h – шаг по x;
ym-1 – предварительно найденное значение (точка m-1).
Затем производится коррекция, то есть происходит приближение, которое вычисляется по формуле [36]:
y(i)m+1 = ym + h*[f(xm,ym) + f(xm+1,y(i-1)m+1)]/2; (9)
где y(i)m+1 – скорректированное значение искомой функции в точке m+1;
y(i-1)m+1 – значение искомой функции, вычисленное на предыдущем шаге коррекции.
Процесс коррекции – итерационный, он прекращается, когда следующее скорректированное значение отличается от предыдущего на пренебрежимо малую величину.
После того, как найдено значение искомой функции можно найти ошибку ограничения [36]:
E = 0.2*[y(0)m – y(i)m]. (10)
где E – ошибка ограничения;
y(0)m – прогнозное значение искомой функции на текущем шаге, вычисленное по формуле (8)
y(i)m – найденное значение искомой функции с помощью коррекции прогнозного значения на i-й итерации коррекции.
Далее, согласно рекомендациям [36] окончательное и более точное приближение можно найти по формуле:
ym = y(i)m + 0.2*(y(0)m – y(i)m). (11)
Проблемой остается выбор шага h. Очевидно, что чем меньше шаг h, тем быстрее сойдется итерационный процесс коррекции. В то же время, если взять малое значение h, то количество точек будет велико. Если же, напротив, взять большое h, то потребуется много итераций для каждой точки, но само количество точек будет меньше.
Кроме этого, известно, что уменьшение шага h уменьшает ошибку ограничения E.
В источнике [36] указывается, что если для какой то точки ошибка ограничения превысила допустимое значение, то величину шага следует уменьшить, и продолжать вычисления дальше. Для этого необходимо снова прибегнуть к одноступенчатому методу, а затем перейти к прогнозу и коррекции. Однако для данной системы важно, чтобы шаг оставался постоянным на всем промежутке. Поэтому в нашем случае если возникнет необходимость уменьшить шаг, то придется начинать расчет с начала всего интервала с меньшим шагом.
Также существуют соображения [37], согласно которым оптимальное число итераций равно двум. То есть, если для расчета понадобилось более двух итераций, то величину шага можно уменьшить.
Итак, численные методы решения третьего и четвертого уравнения системы (4) состоят в следующем.
1. О д н о с т у п е н ч а т ы й м е т о д р е ш е н и я
Изначально известны начальные значения температуры и влажности воздуха: TA0, CA0. Используя эти значения, находятся значения TA1 и CA1 для первого шага. Для этого используется метод Рунге – Кутта, формулы (7) которого для влажности воздуха будут иметь вид:
где CA1 – влажность воздуха в первой точке;
CA0 – начальная влажность воздуха (в нулевой точке);
CM0 – начальная влажность слизистой.
где TA1 – температура воздуха в первой точке;
TA0 – начальная температура воздуха (в нулевой точке);
TM0 температура слизистой в нулевой точке;
CM0 – начальная влажность слизистой;
CA0 – начальная влажность воздуха.
Сложностью метода Рунге – Кутта является вычисление ошибки ограничения. В отличие от методов прогноза и коррекции, где ошибка ограничения получается как побочный продукт вычислений, для вычисления ее в методе Рунге – Кутта необходимо провести еще более трудоемкие вычисления, чем для нахождения приближенного значения.
2. М е т о д п р о г н о з а и к о р р е к ц и и
Когда будут найдены температура и влажность воздуха для первой точки, можно будет воспользоваться методом прогноза и коррекции: (8), (9).
Прогноз для влажности воздуха будет иметь вид:
CA(0)m = CA m-2 + 2h*P*kC*(CM m-1 – CA m-1)/(A*v). (14)
где m – номер текущей точки;
CA(0)m – прогнозное значение влажности воздуха в текущей точке;
CA m-2 – влажность воздуха в точке m-2;
CM m-1 – влажность слизистой в точке m-1;
CA m-1 – влажность воздуха в точке m-1.
Прогноз для температуры воздуха имеет следующий вид:
TA(0)m = TAm-2 + 2h* P*[hC*(TM m-1 – TA m-1) + CPW*kC*(CM m-1 – CA m-1)* (TM m-1 – TA m-1)]/(A*p*CPA*v); (15)
где TA(0)m – прогнозное значение температуры воздуха в текущей точке;
TA m-2 – температура воздуха в точке m-2;
Далее необходимо выполнить коррекцию полученных прогнозных значений. Для этого используем формулу коррекции (9).
Если принять, что
CA' = fC(CA, CM), (16)
то формула коррекции для влажности воздуха примет вид:
CA(i)m+1 = CA m + h*[fC(CA m, CM m) + fC(CA(i-1)m+1, CM m+1)]/2; (17)
где i – номер итерации коррекции,
m – номер текущей точки,
fC – функция производной влажности воздуха,
Аналогично, если принять для температуры воздуха
TA' = fT(TA, TM, CA, CM), (18)
то формула коррекции для температуры воздуха будет:
TA(i)m+1=TA m+h*[fT(TA m, TM m, CA m, CM m)+fT(TA(i-1)m+1, TM m+1, CA m+1, CM m+1)]/2; (19)
где fT – функция производной температуры воздуха.
Итерации коррекции продолжаются до тех пор, пока предыдущее значение не будет отличаться от полученного текущего на пренебрежимо малую величину, которую должен задать пользователь. При этом эти итерации могут продолжаться бесконечно долго, если итерационный процесс коррекции расходится. Для этого следует предусмотреть счетчик количества итераций.
Далее для влажности воздуха и температуры вычисляется ошибка ограничения по формуле (10). Затем полученные значения уточняются с учетом ошибки по формуле (11).
И так, поочередно, производится прогноз и коррекция значений для каждой выбранной точки.
Когда процесс расчета значений влажности и температуры воздуха для всего протяжения однородного участка будет закончен, то полученные значения в последней точке участка станут начальными значениями для следующих непосредственно за ним участков, если есть разветвление, или участка, если разветвления нет.
И так последовательно произойдет расчет распределения температуры и влажности воздуха для полной модели тепловлагообмена легких человека.
Из формул (12), (13), (14), (15), (16), (17) ясно, что на каждом этапе нахождения значений температуры (TA) и влажности (CA) воздуха необходимо, чтобы были известны температура (TM) и влажность (CM) слизистой. Нахождение этих значений заключается в решении первых двух уравнений системы (4).
3. О ц е н к а п о г р е ш н о с т и м е т о д а Р у н г е – К у т т а
Если погрешность для метода прогноза и коррекции вычисляется довольно легко, то для метода Рунге – Кутта нет простых способов оценить его ошибку [36]. Как отмечено в источнике [36], ошибка ограничения при использовании метода Рунге – Кутта p-го порядка равна Khp+1, где K – некоторая постоянная. Верхние пределы для K при p = 2, 3 и 4 приведены в статье Ральстона [38]. Вывод формул для этих верхних пределов является вовсе не простой задачей; более того, этот вывод требует параметров, которые отсутствуют в (7).
Чтобы получить оценку ошибки ограничения, воспользуемся способом, изложенным в [36], предварительно внеся в него поправки и уточнения:
Допустим, что Ym есть точное решение при x = x0 + mh. Тогда для метода четвертого порядка (7) будет верно:
Ym = ym(h) + Kh5, (20)
где ym(h) – приближение к точному решению в m-й точке с шагом h.
Теперь, если вычислять то же самое значение Ym с вдвое меньшим шагом h/2, то, а это источник [36] упускает, количество точек на отрезке будет вдвое большим, то есть Ym будет точным решением при x = x0 + n*(h/2), где n – номер точки при расчете с вдвое меньшим шагом (n = 2m). Отсюда
Ym = yn(h/2) + K(h/2)5, (21)
где yn(h/2) приближение к точному решению в n-й точке с шагом h/2.
Вычитая (20) из (21), получаем:
ym(h) – y n(h/2) = – (31*Kh5)/32, (22)
и ошибка ограничения будет равна:
ET = Kh5 = 32*(ym(h) – y n(h/2))/31. (23)
Отметим, что для того, чтобы найти ошибку ограничения метода Рунге – Кутта, необходимо найти решение дважды для каждой точки. Это значительно сложнее формулы (10). Таким образом, общая стратегия расчета значений температуры и влажности воздуха в первой точке следующая:
1) Вычислить TA1, CA1 по формулам (12), (13) с шагом h;
2) Вычислить TA1*, CA1* с шагом h* = h/2 (формулы для этого те же: (12) и (13));
3) Вычислить TA2*, CA2* с шагом h*;
4) Вычислить ошибку ограничения ETA1 температуры воздуха:
ETA1 = 32*(TA1 – TA1*)/31, (24)
5) Вычислить ошибку ограничения ECA1 влажности воздуха:
ECA1 = 32*(CA1 – CA1*)/31. (25)
Заключение
Как было отмечено, численный расчет основной системы уравнений (4) происходит как "одновременное" решение всех четырех уравнений, входящих в ее состав. Под "одновременным" решением понимается нахождение решения системы (температуры воздуха, влажности воздуха, температуры слизистой, влажности слизистой) для каждой точки лишь после того, как найдено решение для предыдущей точки. То есть при последовательном нахождении решения для точек необходимо получить решение всех уравнений, входящих в состав системы, а затем переходить к нахождению решения для следующей точки.
Несмотря на то, что расчет системы должен происходить "одновременно", можно выделить две составляющие процесса решения системы. Это разбиение на составляющие основано на том, что для расчета температуры и влажности воздуха применяется один численный метод, а для расчета температуры и влажности слизистой другой.
ЛИТЕРАТУРА
- Вериго В.В. Системные методы в космической биологии и медицине. М.: Наука, 1987. 212 с.
- Прогресс биологической и медицинской кибернетики / под ред. акад. Л.И.Берга. М.: Медицина, 1974. 482 с.
- Судаков К.В. Функциональные системы организма как объект физиологического анализа // Вестник АМН СССР. 1985. N 3. С. 3-11.
- Лищук В.А., Потемкина И.С. Организация медицинских знаний и обеспечение решений современными алгоритмическими методами // Вестник АМН СССР. 1988. N 8. С. 47-52.
- West J.M., Skytte J. Anatomical modelling with computer aided design // Computers and biomedical research. 1986. Vol. 19. P. 535-542.
- Hanna L.M., Scherer P.W. Measurement of Local Mass Transfer Coefficients in a Cast Model of the Human Upper Respiratory // Biomechanical Engineering. 1986. Vol. 108. P. 12-18.
- Massad E., Engel A.B., and Nicolelis A.L. A Mathematical Model for Spirometry // Computers and biomedical research. 1987. Vol. 20. P. 105-112.
- Moslehi F., Ligas J.R., Pisani M.A. and Epstein M.A.F. The unsteady form of the Bernoulli eguation for estimating pressure drop in the airways // Respiration Physiology. 1989. Vol. 76. P. 319-326.
- Elad D., Kamm R.D., and Shapiro A.H. Mathematical simulation of forced expiration // J. Appl. Physiol. 1988. Vol. 65 (1). P. 14-25.
- Wilson T.A., Fredberg J.J., Rodarte J.R. and Hyatt R.E. Interdependence of regional expiratory flow // J. Appl. Physiol. 1985 . Vol. 59 P. 1924-1928.
- Humphrey J.D. A possible role of the pleura in lung mechanics // J. Biomechanics. 1987. Vol. 20. N. 8. P. 773-777.
- Sprigings E.J. Simulation of the force enhancement phenomenon in muscle // Comput. Biol. Med. 1986. Vol. 16. N. 6. P. 423-430.
- Bean J.C., Chaffin D.B. and Schultz A.B. Biomechanical model calculation of muscle contraction forces: a double linear programming method // J. Biomechanics. 1988. Vol. 21. N. 1. P. 59-66.
- Kelman G.R. Digital computer subroutine for the conversion oxygen tension into saturation // J. Appl. Physiol. 1966. Vol.21. P. 1375-1376.
- Kelman G.R. Digital computer procedure for the conversion of PCO2 into blood CO2 content // Resp. Physiol. 1967. Vol.3. P. 111-115.
- Gronlund J., Garby L., Lorenzen A.G. and Carter A.M. An improved algorithm and a computer program for the analysis of capillary gas exchange // Acta Physiol Scand. 1986. Vol. 126. P. 259-270.
- Peterman B.F., Longtin A. Multicompartment model of lung dynamics // Computers and Biomedical Research. 1984. Vol. 17. P. 580-589.
- Paiva M., Engel L.A. Model analysis of intra-acinar gas exchange // Resp. Physiol. 1985. Vol. 62. P. 257-272.
- Hanna L.M., Scherer P.W. A theoretical model of localized heat and water vapor transport in human respiratory tract // Asme Journal of Biomeechanical Engineering. 1986. V. 108. P. 19-27.
- Hanna L.M., Scherer P.W. Regional control of local airway heat and water vapor losses // J. Appl. Phisiol. 1986. V. 61. P. 624-632.
- Любимов Г.А., Скобелева И.М. Влияние физических параметров легких на форму кривой поток-обьем форсированного выдоха // Физиология человека. 1992. Т. 18. С. 32-42.
- Гродинз Ф. Теория регулирования и биологические системы. М.: Мир, 1966. 234 с.
- Computer and control in clinical medicine / Ed. by E.R.Carson and D.G.Cramp. New York: Plenum Press, 1985. 263 p.
- Rummel J.A. Integrated simulation of multiple physiological systems in the space life sciences: Report of Johnson Space Center NASA. Houston. 1979. 70p.
- Шабельников В.Г., Дьяченко А.И., Вериго В.В. Возможности оценки силы тяжести на легочное кровообращение по показателям газообмена // Моделирование систем в биологии и медицине. Прага: ЧСНТО, 1980. Т. 8. С.240.
- Ульянычев Н.В. Модель внешнего дыхания человека. Благовещенск, 1990. 57с.
- Дьяченко А.И., Шабельников В.Г. Математические модели действия гравитации на функции легких. М.: Наука, 1985. 274с.
- Шик Л.Л., Токарева Е.М. Количественный анализ зависимости артериализации крови от неравномерности вентиляционно-перфузионных отношений // Бюлл. экспер. биол. и мед. 1980. Т.89. С. 267-263.
- Macklem P.T. Action of the diaphragm on the rib cage // Proc. 28-th Int. Congr. Phisiol. Sci. Budapest. 1981. V. 10. P. 47-61.
- Rochester D.F. The diaphragm: Contractile properties and fatique // J. Clin. Invest. 1985. V. 75. P. 1397-1402.
- Любимов Г.А. Механика органов дыхания // Биомеханика кровообращения, дыхания и биологических тканей. Рига: Зинатне. 1981. С. 109-120.
- Drazen J.M. Physiological basis and interpretation of indices of pulmonary mechanics // Environ. Health Persp. 1984. V. 56. P. 3-9.
- Вейбель Э. Р. Морфометрия легких человека /Пер. с англ. Н. Н. Вольберг. М.:1970. С.174.
- Турчак Л.И. Основы численных методов: Учеб. пособие. М.: Наука. Гл. ред. физ-мат. лит., 1987. 320с.
- Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Лаборатория Базовых Знаний, 2000. 624 с.: ил.
- Мак-Кракен Д., Дорн У. Численные методы и программирование на фортране. М.: Мир, 1977 584с.
- Hull T. E., Greemer A. L. Efficiency of predictor – corrector procedures. Houston. 1963. 412.
- Ralston and Wilf, Mathematical methods for digital computers. Wiley, 1960. ch. 9.
- Вержбицкий В.М. Численные методы (линейная алгебра и нелинейные уравнения): учеб. пособие для вузов. М.: Высш. шк., 2000. 266с.: ил.
- Методы программирования. Учебное пособие. 2-ое издание / под ред. др. ф-м. наук Г.А.Угольницкого. М.: Вузовская книга, 2000. 280 с.
- Прата Стивен Язык программирования С. Лекции и упражнения: Пер. с англ. / Стивен Прата. К.: Издательство "ДиаСофт", 2000. 432 с.
Формулы для температуры воздуха TA1 будут выглядеть аналогично предыдущим:
1. Вериго В.В. Системные методы в космической биологии и медицине. М.: Наука, 1987. 212 с.
2. Прогресс биологической и медицинской кибернетики / под ред. акад. Л.И.Берга. М.: Медицина, 1974. 482 с.
3. Судаков К.В. Функциональные системы организма как объект физиологического анализа // Вестник АМН СССР. 1985. N 3. С. 3-11.
4. Лищук В.А., Потемкина И.С. Организация медицинских знаний и обеспечение решений современными алгоритмическими методами // Вестник АМН СССР. 1988. N 8. С. 47-52.
5. West J.M., Skytte J. Anatomical modelling with computer aided design // Computers and biomedical research. 1986. Vol. 19. P. 535-542.
6. Hanna L.M., Scherer P.W. Measurement of Local Mass Transfer Coefficients in a Cast Model of the Human Upper Respiratory // Biomechanical Engineering. 1986. Vol. 108. P. 12-18.
7. Massad E., Engel A.B., and Nicolelis A.L. A Mathematical Model for Spirometry // Computers and biomedical research. 1987. Vol. 20. P. 105-112.
8. Moslehi F., Ligas J.R., Pisani M.A. and Epstein M.A.F. The unsteady form of the Bernoulli eguation for estimating pressure drop in the airways // Respiration Physiology. 1989. Vol. 76. P. 319-326.
9. Elad D., Kamm R.D., and Shapiro A.H. Mathematical simulation of forced expiration // J. Appl. Physiol. 1988. Vol. 65 (1). P. 14-25.
10. Wilson T.A., Fredberg J.J., Rodarte J.R. and Hyatt R.E. Interdependence of regional expiratory flow // J. Appl. Physiol. 1985 . Vol. 59 P. 1924-1928.
11. Humphrey J.D. A possible role of the pleura in lung mechanics // J. Biomechanics. 1987. Vol. 20. N. 8. P. 773-777.
12. Sprigings E.J. Simulation of the force enhancement phenomenon in muscle // Comput. Biol. Med. 1986. Vol. 16. N. 6. P. 423-430.
13. Bean J.C., Chaffin D.B. and Schultz A.B. Biomechanical model calculation of muscle contraction forces: a double linear programming method // J. Biomechanics. 1988. Vol. 21. N. 1. P. 59-66.
14. Kelman G.R. Digital computer subroutine for the conversion oxygen tension into saturation // J. Appl. Physiol. 1966. Vol.21. P. 1375-1376.
15. Kelman G.R. Digital computer procedure for the conversion of PCO2 into blood CO2 content // Resp. Physiol. 1967. Vol.3. P. 111-115.
16. Gronlund J., Garby L., Lorenzen A.G. and Carter A.M. An improved algorithm and a computer program for the analysis of capillary gas exchange // Acta Physiol Scand. 1986. Vol. 126. P. 259-270.
17. Peterman B.F., Longtin A. Multicompartment model of lung dynamics // Computers and Biomedical Research. 1984. Vol. 17. P. 580-589.
18. Paiva M., Engel L.A. Model analysis of intra-acinar gas exchange // Resp. Physiol. 1985. Vol. 62. P. 257-272.
19. Hanna L.M., Scherer P.W. A theoretical model of localized heat and water vapor transport in human respiratory tract // Asme Journal of Biomeechanical Engineering. 1986. V. 108. P. 19-27.
20. Hanna L.M., Scherer P.W. Regional control of local airway heat and water vapor losses // J. Appl. Phisiol. 1986. V. 61. P. 624-632.
21. Любимов Г.А., Скобелева И.М. Влияние физических параметров легких на форму кривой поток-обьем форсированного выдоха // Физиология человека. 1992. Т. 18. С. 32-42.
22. Гродинз Ф. Теория регулирования и биологические системы. М.: Мир, 1966. 234 с.
23. Computer and control in clinical medicine / Ed. by E.R.Carson and D.G.Cramp. New York: Plenum Press, 1985. 263 p.
24. Rummel J.A. Integrated simulation of multiple physiological systems in the space life sciences: Report of Johnson Space Center NASA. Houston. 1979. 70p.
25. Шабельников В.Г., Дьяченко А.И., Вериго В.В. Возможности оценки силы тяжести на легочное кровообращение по показателям газообмена // Моделирование систем в биологии и медицине. Прага: ЧСНТО, 1980. Т. 8. С.240.
26. Ульянычев Н.В. Модель внешнего дыхания человека. Благовещенск, 1990. 57с.
27. Дьяченко А.И., Шабельников В.Г. Математические модели действия гравитации на функции легких. М.: Наука, 1985. 274с.
28. Шик Л.Л., Токарева Е.М. Количественный анализ зависимости артериализации крови от неравномерности вентиляционно-перфузионных отношений // Бюлл. экспер. биол. и мед. 1980. Т.89. С. 267-263.
29. Macklem P.T. Action of the diaphragm on the rib cage // Proc. 28-th Int. Congr. Phisiol. Sci. Budapest. 1981. V. 10. P. 47-61.
30. Rochester D.F. The diaphragm: Contractile properties and fatique // J. Clin. Invest. 1985. V. 75. P. 1397-1402.
31. Любимов Г.А. Механика органов дыхания // Биомеханика кровообращения, дыхания и биологических тканей. Рига: Зинатне. 1981. С. 109-120.
32. Drazen J.M. Physiological basis and interpretation of indices of pulmonary mechanics // Environ. Health Persp. 1984. V. 56. P. 3-9.
33. Вейбель Э. Р. Морфометрия легких человека /Пер. с англ. Н. Н. Вольберг. М.:1970. С.174.
34. Турчак Л.И. Основы численных методов: Учеб. пособие. М.: Наука. Гл. ред. физ-мат. лит., 1987. 320с.
35. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Лаборатория Базовых Знаний, 2000. 624 с.: ил.
36. Мак-Кракен Д., Дорн У. Численные методы и программирование на фортране. М.: Мир, 1977 584с.
37. Hull T. E., Greemer A. L. Efficiency of predictor – corrector procedures. Houston. 1963. 412.
38. Ralston and Wilf, Mathematical methods for digital computers. Wiley, 1960. ch. 9.
39. Вержбицкий В.М. Численные методы (линейная алгебра и нелинейные уравнения): учеб. пособие для вузов. М.: Высш. шк., 2000. 266с.: ил.
40. Методы программирования. Учебное пособие. 2-ое издание / под ред. др. ф-м. наук Г.А.Угольницкого. М.: Вузовская книга, 2000. 280 с.
41. Прата Стивен Язык программирования С. Лекции и упражнения: Пер. с англ. / Стивен Прата. К.: Издательство "ДиаСофт", 2000. 432 с.