Numerical solution of the direct problem of electroimpedance tomography in a complete electrode formulation.pdf Введение На сегодняшний день существует несколько томографических методов: рентгеновские, ультразвуковые, магнитно-резонансные. Наиболее многообещающим методом является электроимпедансная томография (ЭИТ). Этот метод имеет несколько преимуществ. Во-первых, в отличие от рентгеновского метода, пациент не подвергается облучению. Во-вторых, этот метод можно использовать для длительного наблюдения за деятельностью внутренних органов человека. В-третьих, он очень экономичен и подходит для использования в портативных устройствах [1-3]. Для повышения точности и достоверности медицинских исследований можно комбинировать ЭИТ и другие методы диагностики. Существует два типа задач ЭИТ [4]: - поиск значений электрических потенциалов на поверхности биологического объекта при заданном распределении удельной проводимости с и подаваемого тока I (прямая задача); - восстановление неизвестного распределения проводимости с при заданном токе и значениях потенциалов на поверхности биологического объекта (обратная задача). Математические задачи, возникающие в области ЭИТ, относятся к классу нелинейных некорректно поставленных обратных задач. Численное решение таких задач представляют большую сложность и требует разработки специальных методов. При численном решении дифференциальные задачи обычно приводят к дискретному аналогу с помощью сеточных методов. Выделяют следующие группы этих методов: конечно-разностные (МКР), конечно-элементные (МКЭ), гранично-элементные (МГЭ), конечно-объемные (МКО) и др. Большинство дифференциальных задач решается с помощью МКЭ, преимуществом которого является гибкость, когда нужно быстро улучшить качество аппроксимации без сложных изменений алгоритма. Однако в работе [5] подчеркивалось, что в некоторых вариантах МКЭ не выполняются условия непрерывности электрического тока, а значит, не соблюдается закон сохранения на элементе сетки и поверхности в целом. Результат сравнения МКЭ и МКО в работе [5] показывает, что численное решение, полученное методом МКО, близко к точным результатам даже на очень грубой сетке, и особенно в граничных узлах. Это очень полезно для проблемы ЭИТ. В работе [6] построены конечно-объемные аппроксимации для численного решения прямой задачи ЭИТ на неструктурированных сетках. В качестве конечных объемов были выбраны барицентрические ячейки, ячейки Дирихле-Вороного и 7 Математика / Mathematics треугольные конечные объемы. На основе сравнительного анализа полученных численных решений для построенных разностных схем с точным решением и с численным решением на основе МКЭ установлено, что разностные схемы обеспечивают высокую точность расчетов на барицентрических ячейках, ячейках Дирихле-Вороного. Результаты, полученные с использованием треугольных разностных схем конечного объема, оказались хуже, чем в указанных видах конечных объемов, но стоит отметить, что они позволяют использовать граничные условия Неймана без погрешности аппроксимации. В статье [7] авторы используют метод конечных объемов с барицентрическими ячейками для решения прямых и обратных задач ЭИТ. Результаты численных экспериментов показывают, что численные решения могут хорошо приближаться к точным решениям при уменьшении размера сетки, а точность численных решений обратных задач ЭИТ зависит от выбора начальных значений для итерационного процесса. По сравнению с МКЭ МКО требует меньше времени для достижения требуемой точности. Хотя обратная задача ЭИТ некорректна, полученный в [7] алгоритм сходится, и его устойчивость при граничных условиях Неймана лучше, чем при условиях Дирихле. Однако требуется дополнительный теоретический анализ сходимости и устойчивости. В публикации [8] предлагаются результаты анализа чувствительности внутреннего распределения электрического потенциала в зависимости от небольших изменений электрической проводимости в объекте и / или контактных сопротивлений поверхностных электродов. При дискретизации дифференциальных постановок задач ЭИТ используется МКЭ. Авторы статьи указывают, что полезная информация для решения обратных задач ЭИТ содержится в основном в небольшой приповерхностной области объекта вблизи электродов. Поэтому использование полной электродной модели способствует лучшему пониманию, как решение обратной задачи EIT обнаруживает изменения в неоднородных средах с аномалиями и без них. Целью данной работы является построение математической постановки и численного метода для решения прямой задачи ЭИТ с учетом сопротивления электродов. Математическая постановка представляет собой смешанную граничную задачу для эллиптического уравнения с кусочно-постоянными коэффициентами, причем при задании граничных условий используется специальное интегро-дифференциальное условие на искомую функцию. Физическая постановка задачи Физически прямую задачу ЭИТ можно описать следующим образом [4]. Предполагается, что исследуемый объект (рис. 1, область D) находится в воздухе, имеет достаточно гладкую границу Г и характеризуется кусочно-постоянными значениями коэффициента электропроводности с. Еі, ..., EL отмечают места крепления на поверхности объекта электродов, к которым подается электрический ток относительно низкой [1-4] частоты Һ, ..., IL. Из-за отсутствия источников тока внутри объекта сумма входящего и выходящего токов должна быть равна нулю по закону сохранения заряда. Электроды имеют одинаковые размеры, но могут обладать различными значениями электрического сопротивления z\\, ..., zL из-за различного качества крепления электродов к поверхности объекта. Требу-8 Афанасьева А.А., Старченко А.В. Численное решение прямой задачи электроимпедансной томографии ется найти распределение электрического потенциала внутри объекта и напряжение тока U\\, U2, ..., Ul на электродах. а b Рис. 1. Геометрия объекта с электродами (а); объект с указанием заданных и измеряемых величин (b) Fig. 1. Geometry of an object with electrodes (a); an object with indication of the set and measured values (b) Математическая постановка задачи Математическая постановка рассматриваемой задачи может быть получена с использованием уравнений Максвелла в проводниках [4, 8, 9]: 1 дИ 1 dE 4 п- (1) ---+ rotE = 0,--= rotH--1, с dt с dt c divE = p, divH = 0 и закона Ома для стационарных проводников [4]: j = aE. Здесь E и И - интенсивность электрического и магнитного поля соответственно, c - скорость света, t - время, j - плотность электрического тока, p - объемная плотность заряда электрического тока, с - коэффициент электрической проводимости. В соответствии с физической постановкой задачи можно считать [4, 8], что процесс стационарен I = 0 интенсивность магнитного поля невелика (л»о) и интенсивность источников заряженных частиц внутри ткани биологического объекта равна нулю (р = 0). Тогда имеем rot E = 0, или E = grad (и), где и - потенциал электрического поля. Из второго и третьего уравнений (1) можно записать уравнение, описывающее распределение электрического потенциала внутри исследуемой области: div (a grad и ) = 0. (2) Или с учетом допущения о двумерности изучаемого процесса (3) dX Н *, y )t И (a( *,y )=°' 9 Математика / Mathematics На границе области D , где нет электродов (см. рис. 1), плотность электрического тока j = ст grad(M) равняется 0 в силу отсутствия тока, здесь используются простые условия Неймана: ды дн = 0 (4) где n - внешняя нормаль к границе области D . На электродах рассматриваются граничные условия, соответствующие так называемой полной электродной модели [1-4, 8], в которой учитывается сопротивление электродов: (5) ды ы + Ziст- = Ul, l = 1,..., L . дн Здесь Ui - напряжение на l-м электроде, которое наряду с электрическим потенциалом u(x,y) является искомой величиной. Также на электроде известна сила подаваемого или отводимого тока Il, l = 1, ..., L, причем по закону сохранения L заряда ^ I = 0 (при принятом выше предположении об отсутствии внутренних i=і источников заряда внутри области). Величина силы тока l-м электроде определяется как интеграл от плотности электрического тока по контактной поверхности электрода Еі: I* I* ды Il =J inds =JCT -ds . (6) E Ei П Комбинируя граничные условия (5) и (6) и исключая пока неизвестные величины Ui, можно получить интегро-дифференциальные граничные условия для уравнения (3) на электродах - [Г uds ы + zlст- = - I I uds + ztIt I,l = 1,...,L. дн Ej Таким образом, с учетом принятых предположений и допущений математическая постановка прямой задачи ЭИТ в полной электродной модели имеет вид: д_ дх 0, (х, y) е D; ды ст- = дн L (7) /=1 стЕг=~^\\\\ ыds+ziIiI--, (ху)еEl, 1 = ^..^L; ^Il =0. дн ZlEl уіц J Zl j-f Как было доказано в работе [10], решение рассматриваемой математической по- L становки единственно, если ^ Ul = 0 . Также в [10] указывалось, что полная элекl=і тродная модель в (7) предсказывает экспериментальные данные с ошибкой менее 1%. 10 Афанасьева А.А., Старченко А.В. Численное решение прямой задачи электроимпедансной томографии Численный метод решения Для решения задачи (7) в силу сложности геометрической формы области исследуемого объекта будем использовать метод конечного объема [6] на неструктурированных сетках, в результате применения которого дифференциальная задача заменяется конечно-разностной (системой линейных уравнений). Решение конечно-разностной задачи даст приближенное решение дифференциальной задачи в узлах неструктурированной сетки. На первом этапе применения метода конечного объема область D = D (J Г триангулируется, т.е. покрывается конечным числом непересекающихся треугольных элементов, соседние из которых имеют общие вершины. Граница области исследования аппроксимируется ломаной, представляющей собой совокупность сторон треугольников, последовательно размещенных вдоль границы. Количество вершин треугольников обозначим за N, а количество треугольников -M. Внутри каждого треугольника о(х, у) имеет постоянное значение. В качестве внутреннего конечного объема, не имеющего общих границ с границей области, будем рассматривать барицентрическую ячейку (рис. 2), ограниченную ломаной линией, составленной из отрезков медиан треугольников с общей вершиной P0, которая является центром этой ячейки. Таким образом, множество всех вершин треугольников представляют узлы сетки ®h = jp JN, вокруг которых строятся барицентрические конечные объемы. Введем сеточную функцию vh, являющуюся приближенным решением дифференциальной задачи в узлах сетки, т.е. V и u(Pi) = u(xp,yp),i = 1,...N . Рис. 2. Барицентрическая ячейка внутри триангулированной области Fig. 2. Barycentric cell inside a triangulated area Для получения разностной схемы, связывающей искомые значения сеточной функции vh, предположим, что распределение приближенного решения внутри каждого m-го треугольника с вершинами Po, Pm, Pm+i (отсчет против часовой стрелки; см. рис. 2) можно представить в следующем виде: (8) ѵ^Чх, у) = Vp WP:\\x, y) + Vp VP:\\x, y) + v^, Wm (x, y) , 11 Математика / Mathematics линейные базисные функции, причем они где '¥(Pm\\x, y), x, y), (x, y) P0 P m P m+1 равны 1 только для своей вершины m-го треугольника, для остальных двух они равны 0. В связи с этим имеем i 2 S„ 1 x y , үО»)(x, y) = 1 Pm 2S„ 1 x P0 yP0 1 x Pm yP m 1 x y 1 x Pm+1 yP у Pm+1 1 x Pm+1 yP y Pm+1 т {pm )(x, y) = 1 x P0 yP0 , sm = -m 2 1 x P0 yP0 1 x Pm yPm 1 x Pm yPm 1 x y 1 x Pm+1 yPm+1 где Sm - площадь треугольника AP0PmPm+1. Заметим, что для одного треугольни тт (x, y) 1 m+1 ка базисные функции в сумме дают 1. Проинтегрируем уравнение (3) по одному внутреннему конечному объему V0 с границей L0 и, воспользовавшись формулой Грина [11], получим $ -А ди , ди -ст-dx + o - ду дх dy = 0. (9) При замене производных в (9) для каждого m-го треугольника будем испольdu dv(m) du dv(m) зовать - и-, - и- и с учетом (8) получим для производных постоdx dx dy dy янные значения, которые легко интегрируются [6] на участках границы L0. В итоге для внутренних узлов неструктурированной сетки разностная схема для нахождения сеточных значений будет выглядеть следующим образом: -^-*0 і CT m m=1 m L VP ((yPm - yPm+1 f +(xPm+1 - xPm ^ + VPm ((yPm+1 - yP0 ) (yPm - yPm+1 ) + (^0 - xPm+1 ) (xPm+1 ~ xPm )) + = 0, P0 Efflj . (10) +VPm+1 ((yP0 - yPm ) (yPm ~ yPm+1 ) + (Ч - ^0 ) (^ - xPm )) M0 - количество треугольников в барицентрической ячейке с общей вершиной Р0. Суммирование выполняется по всем треугольным элементам сетки с общей вершиной Р0, находящейся внутри области D, причем когда значение индекса m + 1 становится больше M0, нужно его взять равным 1. Получим разностную схему для значений сеточной функции в узлах, лежащих на границе области D. Пусть Р0 - такой узел, М1 и М2 - соседние узлы сетки, лежащие на границе слева и справа от узла Р0 соответственно (рис. 3). Рассмотрим конечный объем Vo, границу которого L0 составляют отрезки медиан треугольников APqM2М3 и APqM3М1: C2m2, m2C3, C3m1, m1C1, и участки границы области C1P0 и P0C2 (см. рис. 3). Интегрируя уравнение (3) по этому конечному объему и переходя к формуле Грина [11] с заменой производных, получим 12 Афанасьева А.А., Старченко А.В. Численное решение прямой задачи электроимпедансной томографии дп дп дп Рис. 3. Граничный конечный объем (M1P0 и P0M2 - часть границы области D) Fig. 3. Boundary finite volume (M1P0 and P0M2 are part of the boundary of area D) |%£*+ (%£*+ [>£*♦ =о.р0 «т, (id L. Первые два интеграла приближенно представляем, как и для внутренних конечных объемов, пользуясь линейным видом функции v(m)(x, у) (8). Для двух последних интегралов, которые рассматриваются вдоль границы, будем использовать принятые граничные условия (7). В случае, если отрезки MP и / или РМ2 являются границей контакта с воздухом, то соответствующие интегралы равны 0. Если же эти отрезки представляют границу контакта с l-м электродом, то fР ди 1 а ds = Ро Q| f fE uds + Zjlj j 1-1 uds, дп zA 1 zr C , Iе2 ди 1 а ds = |C2 Po| f f uds + ztIt ]-1 f uds. дп ZlEl ^ > zl fPo Здесь для вычисления интегралов воспользуемся заменой подынтегральной „( т-) функции сеточной ѵ(™2 ) (x, у) или v(m- ) (x, y), которые имеют линейную зависимость от x и у. Поэтому f uds и ГРо v(mi) II CO ^3 Co fPo L, +: C fc, Ci ^ 1 . 1 if '3 1 Аналогично uds и |PA| ( - Vp + - V fPo 44 Po 4 V - v. s - s L (s - SM, ) ds = |C1 P0 M v +3vP Для вычисления интеграла uds по контактной поверхности всего электроfE, да, которая образуется несколькими сторонами граничных треугольников, вос-13 Математика / Mathematics пользуемся формулой трапеций. В итоге получим разностную схему для граничных узлов сетки: Mo I + VPo (j.VPm - УҒт+1 ) +(XPm+1 ~ XPm ) +VPm ((yPm+! - УРО ) (yPm ~ ^+1 ) + (XPo ~ 4+1 ) (4+1 ~ XPm )) +VPm+1 ((yP, - yPm )(У Pm ~~ yPm+1 ) + (XPm - XPo )(4+1 - XPm )) \\( л K -1 Po m=1 m L U^I(% + % )WkNk+^+^-Vm^3p z I 2 E I Nt N“' +11 E 4 h |P0C2| 2 E, , , V 1 k=1 ( , K -1 (12) z K -1 У I (vNk + VNk+i ) INkNk+11 ZEl 3VPo + VM 2E^E Nk V 1 k=1 E = ° Po e Уй • Здесь ji = 1, если M1Po e Et, иначе ji = 0; j2 = 1, если PoM2 e E,, иначе j = 0; K - количество узлов сетки на электроде с номером l. Исследование аппроксимационных свойств полученной разностной схемы (11)-(12) с помощью формулы Тейлора затруднено из-за нерегулярной структуры расчетной сетки. Устойчивость (11)-(12) может быть исследована с помощью мажоранты Гершгорина [12]. Построение расчетной сетки В настоящее время существует ряд программ, разработанных для создания сеток. Они позволяют строить вычислительные сетки для произвольно сложных геометрических объектов. В качестве примера можно указать следующие генераторы сеток: Gambit, Tgrid, GMSH и т.д. В данной работе расчетные сетки построены с помощью программы Gambit. Идея работы в Gambit заключается в том, что моделирование объекта происходит постепенно: сначала определяются точки, затем строятся линии, после чего строится сетка и определяются граничные условия. А в дальнейшем получившуюся сетку можно будет экспортировать в текстовый файл и продолжить работу с данными. В качестве тестовой задачи ЭИТ рассмотрим смешанную граничную задачу для уравнения эллиптического типа с кусочно-постоянными коэффициентами (7) в круге единичного радиуса. Электроды полуширины w = 0.25 расположены при Ф = п / 2 и ф = 3 п / 2 . Построим для рассматриваемой области неструктурированную треугольную расчетную сетку со сгущением на электродах с помощью сеточного генератора Gambit (рис. 4). В статье [13] проведен анализ выбора расчетной сетки. В результате сравнения четырех сеток: без сгущения, со сгущением на электродах, со сгущением на границе с воздухом, со сгущением по всей границе, - был сделан вывод, что лучше использовать сетки со сгущением на электродах или, что еще лучше, со сгущением узлов ко всей границе. 14 Афанасьева А.А., Старченко А.В. Численное решение прямой задачи электроимпедансной томографии а b Рис. 4. Расчетная сетка: а - сгущение на электродах (1 256 ячейки); b - сгущение ко всей границе области (2 804 ячейки) Fig. 4. Calculation grid: (a) thickening on the electrodes (1256 cells); (b) thickening to the entire boundary of the region (2804 cells) Заметим, что граница круглой вставки радиуса р < 1 с другим значением электрической проводимости проходит через узлы сетки. Также в данной работе с целью оценки влияния на конечный результат были рассмотрены сетки, специально не подстраиваемые под границу вставки. Информация о влиянии такого способа задания электрической проводимости будет полезна при решении обратных задач ЭИТ. Решение задачи в постановке полной электродной модели Для решения тестовой задачи для двух электродов в постановке (7) воспользуемся разностной схемой (11)-(12). Сила тока на верхнем электроде (ф = л /2) равна 1 мА, на нижнем (ф = 3л / 2) -1 мА. Сопротивление на электродах имели одинаковые значения z1 = z2 = 1. Круглая вставка радиуса р = 0.5 имеет электрическую проводимость оі, а окружающее ее кольцо - 02. В расчетах рассматривались следующие значения о2/о1 = 0.01; 0.1; 1; 10; 100. Полуширина электродов во всех случаях была неизменной w = 0.25, и форма электрода совпадала с частью границы круга. Для данной тестовой задачи известно приближенное аналитическое решение в виде рядов Фурье [14], поэтому выполнялось сравнение численного решения с приближенным аналитическим, в котором учитывалось 700 парных членов ряда. Численное решение разностной схемы (11)-(12) выполнялось методом Гаусса с частичным выбором главного элемента. На рис. 5 представлены графики изменения электрического потенциала по границе круга, рассчитанные для различных отношений о2/о1. Из рисунка видно, что наибольшее по модулю значение потенциала имеет место на электродах, причем с ростом значения отношения о2/о1 абсолютная величина потенциала 15 Математика / Mathematics увеличивается, т.е. при уменьшении проводимости внутри области 01 для обеспечения одной и той же величины тока нужно увеличение значений напряжения на электродах. Визуально наибольшее абсолютное отклонение численного решения от аналитического наблюдается на электродах при с2/с1 > і. Рис. 5. Сравнение значений электрического потенциала по окружности для различных отношений 02/01. Значки - численный расчет в узлах сетки, сплошные линии -приближенное аналитическое решение Fig. 5. Comparison of values of the electric potential along the circle for different ratios 02/01. The icons are the numerical calculation in grid nodes; the solid lines show the approximate analytical solution Для получения объективных оценок воспользуемся величиной среднеквадратичного отклонения (RMSE), которая вычисляется по следующей формуле: N RMSE = І £(V - Uan (xi, Уi )) i=1 N где N1 - количество узлов сетки вдоль границы области. В таблице представлены рассчитанные значения RMSE для различных значений отношения электрических проводимостей о2/о1 и различных сеток: со сгущением узлов только на электродах, ко всей границе области, со специальным размещением узлов на границе круглой вставки и без (в этом случае принималось соглашение, что если центры треугольников находятся в круге радиусом р, то их проводимость равнялась 01). Кроме того, для анализа численной сходимости разностной схемы (11)-(12) использовались сетки, в которых каждый треугольник сеток рис. 4 разбивался по серединам сторон на четыре треугольника. 16 Афанасьева А.А., Старченко А.В. Численное решение прямой задачи электроимпедансной томографии Среднеквадратическое отклонение численного решения от аналитического для разных случаев 02/01 1 0.1 0.01 10 100 С выделением круглой вставки Сетка со сгущением на электродах 0.0093 0.0051 0.0048 0.0225 0.0273 Сетка со сгущением на электродах, с разбиением каждого треугольника на 4 0.0025 0.0022 0.0023 0.0078 0.01 Сетка со сгущением на всей границе 0.0065 0.0032 0.003 0.0173 0.0213 Сетка со сгущением на всей границе, с разбиением каждого треугольника на 4 0.0017 0.0017 0.0017 0.0063 0.0082 Без выделения круглой вставки Сетка со сгущением на электродах 0.0085 0.0151 0.0191 0.0136 0.0121 Сетка со сгущением на электродах, с разбиением каждого треугольника на 4 0.0025 0.0059 0.0087 0.0034 0.004 Сетка со сгущением на всей границе 0.0059 0.0088 0.0106 0.0167 0.0215 Сетка со сгущением на всей границе, с разбиением каждого треугольника на 4 0.0017 0.0039 0.0057 0.0057 0.0068 Анализируя представленные значения RMSE, можно заметить, что численные расчеты лучше согласуются с приближенным аналитическим решением при о2/о1 < 1, причем предпочтительнее использовать неструктурированную сетку со сгущением узлов ко всей границе, чем только к электродам. Кроме того, специальное построение сетки под имеющиеся включения с другими значениями электрической проводимости уменьшают ошибку численных расчетов. Заметим однако, что при решении обратной задачи ЭИТ границы неоднородностей по распределению o(x, y), как правило, неизвестны, и это, безусловно, будет вносить дополнительный шум в искомые значения электрической проводимости. Также значения RMSE из таблицы показывают, что при уменьшении размера сетки в два раза для всех рассматриваемых сеток и значений о2/о1 ошибки разностной схемы (11)-(12) уменьшаются, что свидетельствует о ее сходимости к точному решению. Быстрее всего RMSE уменьшается для однородного диска - почти в 4 раза. Это указывает, что порядок точности построенного численного метода близок ко второму. На рис. 6 в качестве иллюстрации согласования численного и приближенного аналитического решений [14] приведены графики изменения электрического потенциала по радиусу круга при ф = п/4 и ф = п/2. Рассматривалась сетка, представленная на рис. 4, b. Из рисунка видно, что величина электрической проводимости вставки оказывает существенное влияние на распределение потенциала. Если в круглой вставке оно практически линейно изменяется по радиусу, то в кольце зависимость потенциала от радиуса приобретает нелинейный характер, подчеркивая влияние проводимости вставки на эту величину. Также нужно заметить, что при низких значениях электрической проводимости включений потребуются специальные усилия по измельчению сетки для повышения качества численного решения задачи. 17 Математика / Mathematics Рис. 6. Сравнение значений электрического потенциала по радиусу для различных отношений 02/01. Значки - численный расчет в узлах сетки, сплошные линии -приближенное аналитическое решение Fig. 6. Comparison of the values of the electric potential along the radius for different ratios 02/01. The icons are the numerical calculation in grid nodes; the solid lines show the approximate analytical solution 18 Афанасьева А.А., Старченко А.В. Численное решение прямой задачи электроимпедансной томографии Заключение В данной работе представлена двумерная математическая постановка для прямой задачи ЭИТ в биологическом неоднородном токопроводящем объекте с учетом сопротивления приложенных к нему электродов, используемых для подачи и съема тока и измерения напряжения. Особенностью этой математической постановки для двумерного уравнения эллиптического типа с кусочно-постоянными коэффициентами является применение интегро-дифференциального граничного условия на контактной границе электродов. Для численного решения рассматриваемой задачи разработан эффективный численный метод, опирающийся на использование неструктурированных сеток, сгущающихся к границе области исследования, метод конечных объемов для барицентрических ячеек и метод исключения Гаусса с частичным выбором главного элемента для решения сеточных уравнений. Проведенные численные расчеты для тестовой задачи - неоднородного по электрической проводимости круглого диска с двумя электродами - сравнивались с приближенным аналитическим решением. Для различных отношений коэффициентов электрической проводимости (круглой вставки и кольца) получено хорошее согласование, причем при кратном по площади уменьшении размеров конечных объемов неструктурированной сетки расхождение между численным и приближенным аналитическим решениями уменьшалось, что свидетельствует о сходимости построенной конечно-разностной схемы.

Скачать электронную версию публикации