Математическое моделирование процессов аутостабилизации температуры в клеточной ткани для одномерного случая на основе неявных разностных схем

УДК 519.6

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ АУТОСТАБИЛИЗАЦИИ

ТЕМПЕРАТУРЫ В КЛЕТОЧНОЙ ТКАНИ ДЛЯ ОДНОМЕРНОГО СЛУЧАЯ НА ОСНОВЕ НЕЯВНЫХ РАЗНОСТНЫХ СХЕМ

© А.А. Арзамасцев, Д.В. Залевский

Ключевые слова: математическая модель; аутостабилизация; система дифференциальных уравнений в частных производных; неявная разностная схема; схема Кранка-Николсона; схема Ричардсона; схема Дюфорта-Франкела.

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

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

Биохимический реактор представляет собой емкость, в которой происходит превращение питательного вещества (субстрата) в целевой продукт - биомассу. Реакция протекает в аэробных условиях, и кислород поступает в реакционный объем из газовой фазы в ходе барботажа. Растворенный кислород расходуется микроорганизмами на окисление субстрата (экзогенные) и внутренние (эндогенные) нужды.

Биохимический реактор сообщается с окружающей средой по тепловому каналу через стенку. Стенка биореактора имеет тепловую рубашку, которая используется для поддержания температуры с помощью систем автоматики. В режиме аутостабилизации никакие устройства для поддержания температуры не используются. Таким образом, математическая модель аутостаби-лизации температуры клеточных тканей представляет собой систему дифференциальных уравнений в частных производных, которые связывают между собой значения температуры внутри объекта, концентрации субстрата и кислорода. Наиболее подробно данная математическая модель рассматривалась в работах А.А. Арзамасцева и Е.Н. Альбицкой [1-2]. В работе [3] была представлена математическая модель процесса аутостабилизации температуры в клеточных тканях для случая взаимодействия нескольких биологических объектов.

Рассмотрим применение неявных разностных схем к упрощенному случаю математической модели из [3], сделав следующие допущения:

- значения концентрации субстрата и кислорода постоянны в любой момент времени;

- основным объектом исследования будет являться одномерный стержень.

Основное уравнение математической модели при указанных выше условиях будет выглядеть следующим образом:

dT(x,t) X д T(x, t) дТ(x, t) _ - - + W-:--+ QT .

dt cp dx Начальные условия: T(x,0) = T0.

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

dx

dT

dx dT

ddx

= тцТ(0,t) — a(t)],

-^[T (Li,t) - gi(t)].

Дополнительные условия:

Qt = V max [T (x, t )]H cP

Vmax [T (x, t )]= a1 exp

E,

RT

— a2 exp

E,

RT (x, t).

В этих уравнениях: Т(х, Г) - значение температуры в клетке номер х в момент времени Г; X - коэффициент теплопроводности; с, р - удельная теплоемкость и

плотность клеточной ткани; - скорость конвекционного переноса; QT - распределенная функция источников для температуры; Е, Е - энергии активации; Н - интегральный тепловой эффект биохимической реакции; а1, а2 - предэкспоненциальные множители;

т^шах

V - максимальная удельная скорость биохимических реакций.

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

x=L

1805

Рис. 1. А - явная разностная схема; В - простая неявная разностная схема; С - явная разностная схема Кранка-Николсона

из МхМ точек, что определит шаги по времени и температуре Г и Т, соответственно. Тем самым определяются узлы, в которых будет осуществляться поиск решения. Затем необходимо заменить дифференциальные уравнения в частных производных аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого (г, п)-го узла сетки. Поскольку уравнения в частных производных по определению зависят от производных неизвестных функций по нескольким переменным, то вариантов дискретизации этих уравнений может быть довольно много. Конфигурацию узлов, используемую для разностной записи уравнений в частных производных на сетке, называют шаблоном. Существуют два основных наиболее часто применяемых типа шаблонов: явный и неявный.

Явной называют схему, в которой неизвестное значение искомой функции на (п+ 1)-м шаге стоит только в левой части, а в правой части стоят уже вычисленные ранее значения функций. Если в правой части стоят неизвестные значения функций, то схема - неявная. Таким образом, основные различия между неявными и явными разностными схемами заключаются в выборе точек, на основе которых будет составляться алгебраическая система уравнений для решения. Примеры таких схем представлены на рис. 1.

Неявные разностные схемы обладают определенными преимуществами.

1. Получение разностных уравнений, которые связаны системой алгебраических уравнений.

2. Использование метода Гаусса для решения трехдиагональной матрицы. Система линейных уравнений является хорошо обусловленной, а значит, решение системы уравнений с использованием метода Гаусса является оптимальным.

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

4. Неявные разностные схемы обладают большей сходимостью по сравнению с использованием явной разностной схемы.

Общая схема решения системы дифференциальных уравнений с использованием явных (неявных) разностных схем выглядит следующим образом.

1. Покрываем расчетную область сеткой МхМ узлов. Для решения уравнения следует определить искомую функцию в каждом узле, т. е. отыскать ММ неизвестных.

2. Для каждого внутреннего узла записываем дискретное представление уравнения в частных производных.

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

4. Реализуем разностную схему, т. е. решаем систему из ММ алгебраических уравнений каким-либо подходящим численным методом.

5. Полученное решение системы разностных уравнений считают решением уравнения в частных производных.

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

1) прямая неявная разностная схема;

2) неявная разностная схема с преобразованием через обратную матрицу;

3) неявная схема Кранка-Николсона;

4) схема Ричардсона;

5) схема Дюфорта-Франкела;

6) передаточная схема.

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

1. Прямая неявная разностная схема. Преобразования будем приводить к виду рис 1. А [4]. В уравнении (1) выразим все значения производных через их разностные составляющие, после этого подставим полученные значения в уравнение (1) и осуществим группировку элементов:

дТ(х, I) Т(1,, xJ) - Т(1,xJ) ,

А М '

1806

a2T(X,t) т(t,,Xj+1) + T(t,,Xj- 2T(t,,Xj) ,

dx2 Ax2

aT

dx 2Ax

Подставляем полученные разности в уравнение (1).

T (t,., Xj) - T (tM, Xj) = At

+ W

^Г T (t,, Xj+i) + T (t,, Xj-i) + 2T (t,, Xj) cp^ Ax2

Г T

2Ax

+ Q,.

Произведем группировку элементов в зависимости от значения Т, в результате преобразований получим следующее уравнение:

X W

T (t'-i, xj) = -At| CPAX2 + —X |-T (t,, ^ +

+ At| + — !• T(t,,Xj)-

cpAx2 At

X W - At|-—T (t,, x,.+i) +Qt.

cpAx2 2Ax

Введем дополнительные переменные, которыми обозначим соответствующие коэффициенты перед значениями функции Т.

в =-a,_^+

| cpAx2 2 Ax

A 2X i

B, , = At|-;- + —

j,j | cpAx2 At

л i X Wx

Bj,j+i = -At| зтг—-

cpAx2 2Ax

Обозначим через А следующую матрицу коэффициентов:

А =

Обозначим через С матрицу параметров, которые мы будем находить при решении системы алгебраических уравнений:

' Bi,i Bi,2 0

B2,i B2,2 B2,3

0 B3,i Вз,з

v ...

С =

Г Tj V

Tj-i,2 T

V j-i,n У

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

ратуры и значения температуры, рассчитанные на предыдущем шаге. Обозначим эту матрицу как X.

(т Л

T left

X =

T j,2

V ™ght у

В результате получаем систему алгебраических уравнений вида [4]:

АС = X.

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

2. Неявная разностная схема c преобразованием через обратную матрицу.

Осуществляем преобразования, аналогичные преобразованиям в предыдущей схеме, за исключением выбранных точек [5].

5T

dt

At

a2T

dx

Ax2

aT

dx

2Ax

Осуществляем группировку таким образом, чтобы значение Т(¡м, ^) осталось в левой части уравнения, а

все остальные показатели были перенесены в правую часть. В результате получим следующее:

X W

T(t,+i,Xj) = At| + —L. !• T(t,,Xj+i) +

j | cpAx2 2 Ax у j

+ A|-^r- !• T (t,, x,._i) +

vcpAx2 2 Ax

Г 2X i

+| ^-A К (t', ) +q, ,

( f ^ wx ^

2 + 2 Ax

T (t,+i, Xj) = B • T (t,, Xj) +

At

cpAx 0

2X W

cpAx2

X

2Ax

+Qt .

При этом в качестве матрицы В будут выступать следующие преобразования:

b = м

j,j- | cpAx2 2Ax

j,j | cpAx2 At

bjj+i =Atl

X Wr

cpAx 2Ax

+

1807

После получения выражения (2) осуществим преобразование его таким образом, чтобы выразить значение Т(^, х ) . Для этого осуществим обратное преобразование матрицы В и получим следующее уравнение:

Т &, х.) = Б- • Т&+1, х.) + Б-

А

М

Х

срАх

2Х

_ Е*-

2 + 2 Ах

срАх 2Ах

+я,'

В результате всех преобразований получаем уравнение (3), которое является основным уравнением для неявной разностной схемы с преобразованием через обратную матрицу [5].

3. Неявная схема Кранка-Николсона. Первоначально необходимо осуществить преобразования, которые характерны для явной разностной схемы [6].

дТ

Л

2А

д>ТхО_ Т

Лх

Ах2

дТ

Лх

Ах

Подставим указанные выше разностные аппроксимации в уравнение (1) и осуществим группировку элементов в зависимости от значения функции Т:

Таким образом, решая систему алгебраических уравнений (5), мы будем получать значение температуры для каждой клетки в заданный момент времени [6]. 4. Схема Ричардсона.

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

дТ

Л

2А

д2Т

Лх2 Ах2

дТ

Лх

Ах

В дальнейшем все преобразования осуществляются аналогично методу 1 (простая неявная разностная схема), составляется аналогичная матрица.

5. Схема Дюфорта-Франкела.

Схема Дюфорта-Франкела является преобразованием предыдущей схемы (схема Ричардсона) [8]. Основное отличие данной схемы заключается в способе выражения частной производной по х. Данная разностная схема относится к многоуровневым схемам [9]. Преобразование данной схемы осуществляется точно так же, как и преобразование схемы Ричардсона. Также, в отличие от схемы Ричардсона, схема Дюфорта-Франкела является устойчивой.

х.+1) + В.Т (1М, х.) + С.Т (I,-1, х.+1) = О.

где

ХА ШХА

А. =----

. 2срАх2 4Ах

Б = 1 -

ХА срАх2 ХА

- + -

С =-. 2срАх2 4Ах

Б = (

ХА ЖхА . гг/

-!• т (¡1, х.-1) +

2срАх2 4Ах 1 Iй

+ (_%А

V2срАх2 4Ах

^Ут(I,,х.+{)+1-Та,,х.)+а.

срАх

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

(В0 0 0

А Б С

А,

А В„

Т (11хп ) у

(0 Л

оо

о,

V Бп у

дТ

Л

2А

д2т

Лх

Ах2

дТ

Лх

Ах

6. Передаточная схема.

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

дТ

ЛГ

2А

д2Т(х,I) т(I,+!,х.+!) + Т(I,+!,х.-!) - 2Т(1М,х.),

Лх2 Ах2

дТ

Лх

Ах

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

0

1808

выше неявных разностных схем с явной разностной схемой и сделаем вывод о наиболее точной и производительной схеме.

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

Зададим начальное значение шага по времени в 0,000000001 для всех типов схем и сравним результаты, полученные в ходе работы каждой из схем. Они представлены в табл. 1.

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

Согласно полученным данным, точности неявных разностных схем аналогичны явной разностной схеме, что подтверждается на рис. 2. После этого проведем аналогичные исследования, но будем постепенно уменьшать значение шага по времени для построения разностной сетки. И сравним получаемые значения со значениями явной разностной схемы при точности 0,000000001. В дальнейшем будем приводить только графики сравнения точностей (рис. 3-4).

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

Таблица 1

Сравнение точности явной и неявных разностных схем при dt = 0,000000001, Г = 0,016 ч

Номер клетки Явная разностная схема Прямая неявная разностная схема Неявная разностная схема с преобразованием через обратную матрицу Неявная схема Кранка-Николсона Схема Ричардсона Схема Дюфорта-Франкеля Передаточная схема

1 36,997 36,897 36,787 36,897 36,769 36,820 36,903

2 36,542 36,442 36,654 36,215 36,542 36,050 36,519

3 35,998 35,788 36,651 36,348 36,544 36,051 36,165

4 35,712 35,445 35,263 35,310 35,210 35,321 35,948

5 35,588 35,458 35,521 35,060 35,847 35,320 35,545

6 35,478 35,378 35,321 35,301 35,652 35,542 35,210

7 35,558 35,722 35,695 35,063 35,524 35,541 35,651

8 35,712 36,128 35,012 35,023 35,120 35,321 35,456

9 35,988 36,128 35,560 35,021 35,231 35,321 35,365

10 36,542 36,642 36,514 36,321 36,210 36,654 36,564

11 36,997 36,897 36,561 36,654 36,654 36,623 36,812

Рис. 2. Сравнение точности вычислений неявных разностных схем: 1 - явная схема; 2 - неявная схема; 3 - схема через обратную матрицу; 4 - схема Кренка-Никольсона; 5 - схема Ричардсона; 6 - схема Дюфорта-Франкела

1809

37.5

37 36.5

36 35.5

35 34.5 34

123 4567 8 9 10 11

Рис. 3. Сравнение точности вычислений неявных разностных схем с шагом dt = 0,000001: 1 - явная схема; 2 - неявная схема; 3 - схема через обратную матрицу; 4 - схема Кренка-Никольсона; 5 - схема Ричардсона; 6 - схема Дюфорта-Франкела

38.5

38 37.5

37 36

35.5

35 34.5 34

123456789 10 11

Рис. 4. Сравнение точности вычислений неявных разностных схем с шагом dt = 0,0001: 1 - явная схема; 2 - неявная схема; 3 - схема через обратную матрицу; 4 - схема Кренка-Никольсона; 5 - схема Ричардсона; 6 - схема Дюфорта-Франкела

Таблица 2

Время вычисления каждой схемы

Название схемы / время выполнения (с) Л = 0,00000001; Л = 0,001 Л = 0,0000001; Л = 0,001 Л = 0,000001; Л = 0,001 Л = 0,00001; Л = 0,001 Л = 0,0001; Л = 0,001 & = 0,01; Л = 0,001 & = 0,1; Л = 0,001

Явная разностная схема 75,8 70,5 68,4 62,65 60,64 59,65 50,47

Прямая неявная разностная схема 32,5 30,59 28,65 27,89 25,98 20,32 17,48

Неявная схема с преобразованием через обратную матрицу 35,6 32,51 30,54 28,56 25,56 24,54 22,35

Неявная схема Кран-ка-Николсона 40,5 39,78 33,57 31,47 30,25 28,54 25,12

Схема Ричардсона 46,8 41,25 40,81 39,54 35,84 34,54 33,24

Схема Дюфорта-Франкела 55,19 55,1 52,32 50,54 48,98 45,68 44,46

Передаточная схема 55,89 55,02 52,11 50,85 48,65 45,98 42,54

1810

Рис. 5. Сравнение производительности вычислений при указанных точностях: 1 - Л = 0,1, Л = 0,001; 2 - Л = 0,01, Л = 0,001; 3 - Л = 0,001, Л = 0,001; 4 - Л = 0,00001, Л = 0,001; 5 - Л = 0,000001, Л = 0,001; 6 - Л = 0,0000001, Л = 0,001; 7 - Л = = 0,00000001, Л = 0,001

Рис. 6. Сравнение производительности вычислений при указанных точностях: 1 - Л = 0,1, Л = 0,01; 2 - Л = 0,01, Л = 0,001; 3 - Л = 0,01, Л = 0,001; 4 - Л = 0,00001, Л = 0,01; 5 - Л = 0,000001, Л = 0,01; 6 - Л = 0,0000001, Л = 0,01; 7 - Л = 0,00000001, Л = 0,01

Таким образом, исходя из полученных результатов вычислений, можно сделать вывод о том, что наиболее точными оказались схемы 1-4 (прямая неявная разностная схема; неявная разностная схема с преобразованием через обратную матрицу; схема Кранка-Никол-сона; схема Ричардсона).

Теперь необходимо сравнить производительность работы каждой из схем. Для этого оцениваем количество шагов (время) (табл. 2), которое потребовалось алгоритму для того, чтобы дойти до заданной точки с заданной погрешностью измерений (рис. 5). Аналогичным образом произведем сравнение при значении Л = = 0,01 (рис. 6).

В результате исследования рассмотрено 6 разностных схем и произведена оценка точности и производительности каждой из таких схем. Таким образом, можно сделать вывод о том, что наиболее производительными и точными являются следующие схемы:

1) прямая неявная разностная схема;

2) неявная разностная схема с преобразованием через обратную матрицу;

3) неявная схема Кранка-Николсона.

ЛИТЕРАТУРА

1. Арзамасцев А.А., Альбицкая Е.Н. Математическое моделирование и исследование саморегулирования температуры в популяциях микроорганизмов: два биообъекта // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2009. Т. 14. Вып. 2. С. 370-374.

2. Арзамасцев А.А., Альбицкая Е.Н. Вычислительные эксперименты по моделированию саморегулирования температуры в популяциях микроорганизмов // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2008. Т. 13. Вып. 1. С. 80-83.

3. Арзамасцев А.А., Альбицкая Е.Н. Математическое моделирование саморегулирования температуры в популяциях микроорганизмов: непрерывный процесс // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2007. Т. 12. Вып. 6. С. 709-714.

4. Kaus B.J.P. Numerische Methoden Explicit vs. implicit finite différentes // Johannes Gutenberg University Mainz. Germany, 2011. P. 1-4

5. Privault N. An Elementary Introduction to Stochastic Interest Rate Modeling. 2nd Edition. France: World Scientific, 2012. 242 p.

6. Numerically Solving PDE's: Crank-Nicholson Algorithm. Bologna: Università di Bologna, 2012.

7. Bulirsch R., Stoer J. Fehlerabschatzungen und Extrapolation mit rationalen Funktionen bei verfahren vom Richardson Typus // Numer Math. 1964. V. 6. P. 413-427.

8. Rezzolla L. Numerical Methods for the Solution of Partial Differential Equations // Lecture Notes for the COMPSTAR School on Computa-

1811

tional Astrophysics, 8-13/02/10, Caen, France; Albert Einstein Institute, Max-Planck-Institute for Gravitational Physics, Potsdam, Germany, 2011.

Gottlieb D., Gustafsson B. Generalized Du Fort-Frankel Methods for Parabolic Initial-Boundary Value Problems // SIAM Journal on Numerical Analysis. 1976. V. 13. № 1. P. 129-144.

Mastro M. Financial Derivative and Energy Market Valuation: Theory and Implementation in MATLAB. N. Y., 2013.

Поступила в редакцию 27 октября 2014 г.

Arzamastsev A.A., Zalevskiy D.V. MATHEMATICAL MODELING OF AUTO-STABILIZATION OF TEMPERATURE IN CELL TISSUE FOR ONE-DIMENSIONAL CASE ON THE BASIS OF IMPLICIT DIFFERENCE SCHEMES

A mathematical model with distributed parameters for the description of one-dimensional process of auto-stabilization of the temperature in the cellular tissue, the algorithm for calculating the system of equations based on the implicit difference schemes is developed. The results of model validation in computational experiment allow you to select the most productive implicit difference scheme.

Key words: mathematical model; auto-stabilization; system of differential equations; implicit difference scheme; Crank-Nicolson scheme; Richardson scheme; Dyufort-Frankel scheme.

Арзамасцев Александр Анатольевич, Тамбовский государственный университет им. Г.Р. Державина, г. Тамбов, Российская Федерация, доктор технических наук, профессор, зав. кафедрой компьютерного и математического моделирования, e-mail: arz_sci@mail.ru

Arzamastsev Alexander Anatolyevich, Tambov State University named after G.R. Derzhavin, Tambov, Russian Federation, Doctor of Technics, Professor, Head of Computer and Mathematical Simulation Department, e-mail: arz_sci@mail.ru

Залевский Дмитрий Вадимович, Тамбовский государственный университет им. Г.Р. Державина, г. Тамбов, Российская Федерация, аспирант, кафедра компьютерного и математического моделирования, e-mail: arz_sci@mail.ru

Zalevskiy Dmitriy Vadimovich, Tambov State University named after G.R. Derzhavin, Tambov, Russian Federation, Post-graduate Student, Computer and Mathematical Simulation Department, e-mail: arz_sci@mail.ru

1812