Исследование схемы с симметризованными производными с двумерной переменной по пространству для численного решения уравнения переноса с запаздыванием Сергей Свиридов sergey.sviridov@urfu.ru УрФУ (Екатеринбург) Аннотация Рассматривается уравнение переноса с запаздыванием с простран- ственной переменной второго порядка. Для такого уравнения, с по- зиций принципа разделения конечномерной и бесконечномерной со- ставляющих состояния, строится сеточный метод симметризован- ных производных. Для учета эффекта наследственности приме- няются одномерная и двойная кусочно-линейная интерполяции и экстраполяция продолжением. Доказывается, что рассмотренный метод имеет порядок локальной погрешности O(hx 2 + hy 2 + ∆2 ), где hx , hy – шаги дискретизации по пространственным перемен- ным, ∆ – шаг дискретизации по временной переменной. Иссле- дуются свойства двойной кусочно-линейной интерполяции. Ис- пользуя результаты общей теории разностных схем, установле- ны условия устойчивости предложенного метода. С помощью вло- жения в общую схему численных методов для функционально- дифференциальных уравнений получена теорема о порядке сходи- мости сконструированного алгоритма. Приведен тестовый пример. 1 Введение Уравнение переноса – это уравнение в частных производных первого порядка. Оно описывает переход сохраняющейся скалярной величины в пространстве, обусловленный диффузией. В моделях физики такое уравнение часто называют уравнением конвекции, в моделях биологии – уравнением адвекции. Этот эф- фект может осложняться различными видами запаздывания [1]. В силу сложности объектов, на первый план выходит конструирование численных алгоритмов их решения. Численные методы решения уравнения переноса без запаздывания достаточно хорошо изучены, напри- мер, в работах [3, 2, 7]. Общие вопросы для функционально-дифференциальных уравнений первого порядка рассматривались в ряде работ, к примеру, в [8]. В данной статье численный метод конструируется с позиции принципа разделения состояния системы на конечномерную и бесконечномерную составляющие. По конечномерной составляющей конструируются Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org 326 аналоги методов, известных для объектов без запаздывания. Для предыстории (бесконечномерная со- ставляющая) используется интерполяция с заданными свойствами. Этот подход позволяет эффективно конструировать новые методы и создавать соответствующее программное обеспечение. Отметим, что ранее этот подход применялся к простейшим сеточным алгоритмам, анонсированным в работе [9], и реализован в [10, 11]. Метод симметризованных производных был изучен для одномер- ной пространственной переменной в [6], он показал лучшие результаты среди всех рассмотренных ранее методов. Данная работа продолжает эти исследования, в ней предложена общая методика и приводят- ся доказательства некоторых теорем о порядках сходимости метода симметризованных производных для двухмерной переменной по пространству. Рассмотрим уравнение переноса с запаздыванием по времени:   ∂u ∂u ∂u +a + = f (x, y, t, u(x, y, t), ut (x, y, ·)), (1) ∂t ∂x ∂y 0 6 t 6 T, 0 6 x 6 X, 06y6Y с граничными u(0, y, t) = g1 (y, t), u(x, 0, t) = g2 (x, t), 0 6 t 6 T, 0 6 x 6 X, 06y6Y (2) и начальным условиями u(x, y, s) = ϕ(x, y, s), 0 6 x 6 X, 0 6 y 6 Y, −τ 6 s 6 0. (3) Здесь x, y, t – независимые переменные, u(x, y, t) – искомая функция, ut (x, y, ·) = {u(x, y, t + s), −τ 6 s < 0} – предыстория искомой функции по времени t, τ > 0 – величина запаздывания. Положим, что функционал f , функции g1 (y, t), g2 (x, t), ϕ(x, y, s) и константа a > 0 таковы, что задача (1)–(3) имеет единственное решение. 2 Дискретизация задачи Пусть hx = X/Nx , hy = Y /Ny , введем точки xi = ihx , yk = khy , i = 0, . . . , Nx , k = 0, . . . , Ny . Пусть ∆ = T /M, введем следующие обозначения: tj = j∆, j = −K, . . . , M. Величина τ /∆ = K – положительное число. Обозначим за uji,k приближенные значения функций u(xi , yk , tj ) в узлах получившейся сетки. Для каждого фиксированного индекса i = 0, . . . , Nx , k = 0, . . . , Ny введем дискретную предысторию по времени tj , j = 0, . . . , M : {uli,k }j = {uli,k , j − m 6 l 6 j}. Оператор интерполяции-экстраполяции – это оператор, введенный на наборе всех допустимых предыс- j торий и действующий по правилу I : {uli,k }j → vi,k (·) ∈ Q[tj − τ, tj + ∆]. Здесь Q[α, β] – это набор функций u(s), кусочно-непрерывных на [α, β] с конечным числом точек разрыва первого рода. В точках разрыва будем считать фукнцию непрерывной справа. Определим норму Q = Q[α, β] как ku(·)kQ = max |u(s)|. s∈[α,β] Будем полагать, что, во-первых, оператор интерполяции-экстраполяции липшицев, то есть существует такая константа LI , что для всех предысторий дискретной модели {uli,k }j и {yi,k l }j выполняется j j max |vi,k (t) − wi,k (t)| 6 LI max |uli,k − yi,k l |, tj −τ 6t6tj +∆ j−m6l6j j j где vi,k (·) = I({uli,k }j ), wi,k l (·) = I({yi,k }j ). Во-вторых, положим, что оператор интерполяции-экстраполяции согласован, то есть j vi,k (tl ) = uli,k , l = j − m, . . . , j. Будем говорить, что оператор интерполяции-экстраполяции имеет порядок p, если существуют констан- j ты C1 и C2 такие, что |vi,k (t) − u(xi , yk , t)| 6 C1 maxj−m6l6j |uli,k − u(xi , yk , tl )| + C2 ∆p для всех i, k, j и t ∈ [tj − τ, tj+1 ]. Простейший способ интерполяции – кусочно-линейная функция. Простейший способ экстраполяции – экстраполяция продолжением (см. [12]). Будем использовать кусочно-линейную функцию с экстраполяцией продолжением, которая является липшицевым оператором (LI = 2), согласованна и имеет порядок p = 2 (см. [12]). 327 3 Метод симметризованных производных Рассмотрим метод: 1 j (u − uji,k−1 + uji−1,k − uji−1,k−1 + uj−1 j−1 j−1 j−1 i,k − ui,k−1 + ui−1,k − ui−1,k−1 )+ 4∆ i,k a + (uj − uji−1,k + uj−1 j−1 j j j−1 j−1 i,k − ui−1,k + ui,k−1 − ui−1,k−1 + ui,k−1 − ui−1,k−1 )+ (4) 4hx i,k a j− 12 + (uji,k − uj−1 j j−1 j j−1 j j−1 i,k + ui−1,k − ui−1,k + ui,k−1 − ui,k−1 + ui−1,k−1 − ui−1,k−1 ) = fi− 12 ,k− 21 . 4hy Каждая производная в (4) представлена как 41 суммы приближенных значений функции в узлах сетки, показанных на рис. 1. Величина uji,k может быть выражена из ранее подсчитанных значений, значит, метод симметризованных производных (4) является явным для уравнения (1)–(3). Рис. 1: Узлы для подсчета очередного значения искомой функции j− 1 hx hy ∆ hx hy ∆ hx hy fi− 12,k− 1 = f (xi − , yk − tj − , υ(xi − , yk − , tj − ), υtj − ∆ (xi − , yk − , ·)). 2 2 2 2 2 2 2 2 2 2 2 Здесь используется расширенный оператор интерполяции-экстраполяции дискретной предыстории Iˆ : ({uji,k , uji−1,k , uji,k−1 , uji−1,k−1 }) → υi− j 1 ,k− 1 (·) ∈ Q[tj − τ, tj + ∆]. 2 2 4 Сходимость метода Будем говорить, что метод сходится, если εji,k → 0 при hx → 0, hy → 0 и ∆ → 0 для всех i = 0, . . . , Nx , p k = 0, . . . , Ny и j = 0, . . . , M . Будем говорить, что метод сходится с порядком hpxx +hyy +∆q , если существует p константа C такая, что |εji,k | 6 C(hpxx + hyy + ∆q ) для всех i = 0, . . . , Nx , k = 0, . . . , Ny и j = 0, . . . , M. Порядок сходимости метода зависит, в первую очередь, от порядка локальной погрешности или невязки. Невязка без интерполяции метода (4) – это, по определению, сеточная функция j 1 ψi,k = (u(xi , yk , tj ) − u(xi , yk−1 , tj ) + u(xi−1 , yk , tj ) − u(xi−1 , yk−1 , tj ) + u(xi , yk , tj−1 )− 4∆ −u(xi , yk−1 , tj−1 ) + u(xi−1 , yk , tj − 1) − u(xi−1 , yk−1 , tj−1 ))+ a + (u(xi , yk , tj ) − u(xi−1 , yk , tj ) + u(xi , yk , tj − 1) − u(xi−1 , yk , tj−1 ) + u(xi , yk−1 , tj )− 4hx −u(xi−1 , yk−1 , tj ) + u(xi , yk−1 , tj−1 ) − u(xi−1 , yk−1 , tj−1 ))+ (5) a + (u(xi , yk , tj ) − u(xi , yk , tj−1 ) + u(xi−1 , yk , tj ) − u(xi−1 , yk , tj−1 ) + u(xi , yk−1 , tj )− 4hy −u(xi , yk−1 , tj−1 ) + u(xi−1 , yk−1 , tj ) − u(xi−1 , yk−1 , tj−1 ))− −f (xi− 12 , yk− 21 , tj− 12 , u(xi− 21 , yk− 21 , tj− 12 ), utj− 1 (xi− 21 , yk− 21 , ·)). 2 p Будем говорить, что невязка имеет порядок hpxx + hyy + ∆p∆ , если существует константа C такая, что j px py p∆ |ψi,k | 6 C(hx + hy + ∆ ) для всех i = 1, . . . , Nx , k = 1, . . . , Ny , j = 0, . . . , M − 1. 328 Лемма 1. Если все частные производные для точного решения задачи (1) – (3) существуют и непре- рывны вплоть до 3 порядка включительно, то метод (4) имеет порядок невязки без интерполяции O(h2x + h2y + ∆2 ). Лемма 1 может быть доказана при помощи формулы Тейлора для каждого узла uji,k на сетке. Сведем схему к общей схеме, рассмотренной в [5]. Без ограничения общности, будем рассматривать нулевые граничные условия. Введем следующие обозначения: yj = (uj1,1 , uj2,1 , . . . , ujNx ,1 , uj1,2 , uj2,2 , . . . , ujNx ,2 , . . . , ujNx ,Ny ) ∈ Y = RNx Ny , σx = a∆ a∆ hx , σy = hy , σx > 0, σy > 0, A : Ayj = (uj1,1 − uj0,1 , . . . , uji,k − uji−1,k , . . . , ujNx ,Ny − ujNx −1,Ny ), B : Byj = (uj1,1 − uj1,0 , . . . , uji,k − uji,k−1 , . . . , ujNx ,Ny − ujNx ,Ny −1 ), C : Cyj = (uj1,1 − uj0,0 , . . . , uji,k − uji−1,k−1 , . . . , ujNx ,Ny − ujNx −1,Ny −1 ), Тогда наш метод (4) в векторной форме может быть записан как 1  2yj − Ayj − (2yj−1 − Ayj−1 + yj − Byj − yj−1 + Byj−1 + yj − Cyj − yj−1 + Cyj−1 ) + 4∆ a a + (Ayj + Ayj−1 − Byj + Cyj − Byj−1 + Cyj−1 ) + (Byj − Ayj + Cyj + Byj−1 − Ayj − 1 + Cyj−1 ) = 4hx 4hy = Fj (υ(·)), N ,N N ,N где Fj (υ(·)) = (Fj1,1 (υj1,1 (·)), Fj2,1 (υj2,1 (·)), . . . , Fj x y (υj x y (·))), υ(·) = I({yi,k }j ) ∈ QNx Ny [−τ, ∆]. Здесь V = QNx Ny [−τ, ∆] – интерполяционное пространство, пространство Nx Ny -мерных вектор-функций, где каждый компонент принадлежит пространству Q[−τ, ∆]. Таким образом, после упрощений наша схема принимает следующий вид: (4E − A − B − C + σx (A − B + C) + σy (B − A + C))yj = (6) = −(−4E + A + B + C + σx (A − B + C) + σy (B − A + C))yj−1 + 4∆Fj (υ(·)). Изучим устойчивость метода. Матрица A может быть записана как Nx − 1 Nx Nx + 1 Nx + 2 1 0 0 ... 0 0 0 0 ... 0 0 −1 1 0 ... 0 0 0 0 ... 0 0 0  −1 1 ... 0 0 0 0 ... 0 0   . .. .. .. .. .. .. .. .. .. ..   .  . . . . . . . . . . .  0 0 0 ... −1 1 0 0 ... 0 0 Nx   0 0 0 ... 0 0 1 0 ... 0 0 Nx + 1   0 0 0 ... 0 0 −1 1 ... 0 0 Nx + 2    . .. .. .. .. .. .. .. .. .. ..   . . . . . . . . . . . .  0 0 0 ... 0 0 0 0 ... −1 1 Матрица B: Nx Nx Ny − Nx 1 0 ... 0 0 0 ... 0 ... 0 0 1 ... 0 0 0 ... 0 ... 0  . .. .. .. .. .. .. .. .. ..   . .  . . . . . . . . .  0 0 ... 1 0 0 ... 0 ... 0 Nx   −1 0 ... 0 1 0 ... 0 ... 0 Nx + 1   0 −1 ... 0 0 1 ... 0 ... 0 Nx + 2    . .. .. .. .. .. .. .. ..   . .. . . . . . . . . . .  0 0 ... 0 0 0 ... −1 ... 1 329 Матрица C: Nx + 1 Nx Ny − Nx − 1 1 0 ... 0 0 0 ... 0 ... 0 0 1 ... 0 0 0 ... 0 ... 0  . .. .. .. .. .. .. .. .. ..   . .  . . . . . . . . .  0 0 ... 1 0 0 ... 0 ... 0 Nx + 1   −1 0 ... 0 1 0 ... 0 ... 0 Nx + 2   0 −1 ... 0 0 1 ... 0 ... 0 Nx + 3    . .. .. .. .. .. .. .. ..   . .. . . . . . . . . . .  0 0 ... 0 0 0 ... −1 ... 1 В левой части (6) видим нижнетреугольную матрицу (как линейную комбинацию нижнетреугольных матриц):   1 + σx + σy 0 0 ... 0 0   · 1 + σx + σy 0 ... 0 0     · · 1 + σ x + σ y . . . 0 0   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 + σx + σy 0  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 + σx + σy Обратная матрица для нее может быть записана в виде: 1 0 0 ... 0 0   1+σx +σy 1   · 1+σx +σy 0 ... 0 0   1   · · 1+σx +σy ... 0 0   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1+σx1+σy 0    . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1+σx1+σy Таким образом, схема (6) может быть представлена в виде: yj+1 = Syj + ∆Φ(υ(·)), (7) где матрица S:  1−σx −σy  1+σx +σy 0 0 ... 0 0 0 1−σx −σy · 0 ... 0 0 0    1+σx +σy   1−σx −σy    · · 1+σx +σy . . . 0 0 0   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   1−σx −σy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1+σx +σy 0    1−σx −σy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1+σx +σy 1−σ −σ Собственные числа матрицы S равны λ(S) = 1+σxx +σyy , |λ(S)| < 1. Таким образом, любая степень S n ограничена в любой подчиненной норме. Это означает, что доказана следующая лемма. Лемма 2. Метод (4) устойчив. Определим функцию точных значений для схемы: zj = (u(x1 , y1 , tj ), u(x2 , y1 , tj ), . . . , u(xNx −1 , yNy , tj ), u(xNx , yNy , tj )). Начальные значения модели равны функции точных значений: yj = zj = (ϕ(x1 , y1 , tj ), ϕ(x2 , y1 , tj ), . . . , ϕ(xNx , yNy , tj )), j = −m, . . . , 0. 330 j p Пусть невязка в смысле (5) имеет порядок ∆q + hpxx + hyy , функции Fi,k липшицевы, оператор интерполяции-экстраполяции I имеет порядок погрешности p0 на точном решении. Тогда невязка с ин- p терполяцией в смысле общей схемы [4] имеет порядок погрешности ∆min{p0 ,q} + hpxx + hyy . Определим невязку с интерполяцией как сеточную функцию dn = (zn+1 − Szn )/∆ − Φ(tn , I(zin ), ∆), n = 0, . . . , M − 1. Лемма 3. Пусть выполняются условия Леммы 1. Если используется двойная кусочно-линейная ин- терполяция с экстраполяцией продолжением, то невязка с интерполяцией в общей схеме имеет порядок O(∆2 + h2x + h2y ) Метод сходится с порядком ∆p1 + hpx2 + hpy3 , если существует константа C: kzj − yj kY 6 C(∆p1 + hpx2 + hpy3 ) для всех j = −m, . . . , M . Будем использовать следующий результат. Теорема о сходимости в общей схеме [4]. Пусть метод устойчив, функция Φ липшицева по вто- рому аргументу, оператор интерполяции I удовлетворяет условию липшицевости, погрешность апрок- симации с интерполяцией имеет порядок ∆p1 +hpx2 +hpy3 , где p1 > 0, p2 > 0, p3 > 0, тогда метод сходится с порядком, не меньшим ∆p1 + hpx2 + hpy3 . Таким образом, доказана следующая Теорема. Метод (4) сходится с порядком O(∆2 + h2x + h2y ). 5 Пример Проверим устойчивость и сходимость на следующем примере.   ∂u ∂u ∂u +a + = (1 + τ )(sin πx + cos πy) + aπt(cos πx − sin πy) + u(x, y, t − τ ) − u(x, y, t), ∂t ∂x ∂y где x ∈ [0, 2], y ∈ [0, 2], t ∈ [0, 3], a = 1, τ = 14 с начальным и граничными условиями u(x, y, t) = t(sin πx + cos πy), x ∈ [0, 2], y ∈ [0, 2], t ∈ [−τ, 0], u(0, y, t) = t cos πy, t ∈ [0, 3], y ∈ [0, 2]. u(x, 0, t) = t sin πx, t ∈ [0, 3], x ∈ [0, 2]. Это уравнение имеет точное решение: u(x, y, t) = t(sin πx + cos πy). Максимум абсолютной погрешности подсчитан для различного количества шагов по x, y и t: T Nx = 5, 10, 15, ..., 100, Ny = Nx , K = ∆ = 5, 10, 15, ..., 100. Результаты экспериментов можно видеть на рис. 2-4. Список литературы [1] J. Wu. Theory and Application of Partial Functional Differential Equations. New York: Springer-Verlag, 1996. [2] N. N. Kalitkin. Numerical methods. Second edition. Sankt-Peterburg: BHV-Pererburg, 2011 (in Russian). = Н. Н. Калиткин. Численные методы. 2-е издание. Санкт-Петербург: БХВ-Петербург, 2011. [3] A. A. Samarsky. Theory of differencial schemes. Third edition. M.: Nauka, 1989 (in Russian). = А. А. Са- марский. Теория разностных схем. 3-е издание. М.: Наука, 1989. [4] V. G. Pimenov, A. B. Lozhnikov. Differencial schemes for numerical solution of a heat equation with delay. Trudy IMM, 17(1):178–189, 2011 (in Russian). = В. Г. Пименов, А. Б. Ложников. Разностные схемы численного решения уравнений теплопроводности с последействием. Труды ИММ, 17(1):178–189, 2011. [5] V. G. Pimenov. Differencial methods for solution of functional equations with delay. Ekaterinburg: Izdatelstvo uralskogo universiteta. 2014 (in Russian). = В. Г. Пименов. Разностные методы решения уравнений в частных производных с наследственностью. Екатеринбург: Издательство уральского университета. 2014. 331 Рис. 2: Максимум абсолютной погрешности для рассмотренного примера для различных значений шагов (Nx = 5, ..., 100, Ny = Nx , K = 5, ..., 100) [6] S. V. Sviridov. Investigation of the scheme with symmetrized derivatives for numerical solution of an advection equation with delay. Vestnik Tambovskogo universiteta. Seriya: Estestvennie i tehnicheskie nauki, 20(5):1420—1421, 2015 (in Russian). = С. В. Свиридов. Исследование схемы с симметризованными производными для численного решения уравнения переноса с запаздыванием. Вестник Тамбовского университета. Серия: Естественные и технические науки, 20(5):1420—1421, 2015. [7] I. B. Petrov, A. I. Lobanov. Lectures on computational mathematics. M.: Binom, 2006 (in Russian). = И. Б. Петров, А. И. Лобанов. Лекции по вычислительной математике. М.: Бином, 2006. [8] Z. Kamont, V. Chernous. Implicit differencial methods for functional differencial equations of Gamilton– Jacobi. Sib. zhurn. vychisl. matematiki, 12(1):57–70, 2009 (in Russian). = З. Камонт, В. Черноус. Неявные разностные методы для функциональных дифференциальных уравнений Гамильтона–Якоби. Сиб. журн. вычисл. математики, 12(1):57–70, 2009. [9] L. S. Volkanin. Numerical solution of an advection equation with delay. Teoriya upravlenija i matematicheskoe modelirovanie. Izhevsk, 12–13, 2012 (in Russian). = Л. С. Волканин. Численное ре- шение уравнения переноса с эффектом наследственности. Теория управления и математическое моделирование. Ижевск, 12–13, 2012. [10] S. I. Solodushkin. The differencial scheme for numerical solution of an advection equation with delay. Izvestiya vyshikh uchebnykh zavedenij. Matematika, 10:77–82, 2013 (in Russian). = С. И. Солодушкин. Разностная схема для численного решения уравнения переноса с последействием. Известия высших учебных заведений. Математика, 10:77–82, 2013. [11] V. G. Pimenov, S. V. Sviridov. Grid methods for solving of an advection equation with delay. Vestnik Udmurtskogo Universiteta. Matematika, Mekhanika. Komputernye nauki, 3:59–74, 2014 (in Russian). = В. Г. Пименов, С. В. Свиридов. Сеточные методы решения уравнения переноса с запаздыванием. Вестник Удмуртского университета. Математика. Механика. Компьютерные науки, 3:59–74, 2014. [12] A. V. Kim, V. G. Pimenov. i-smooth analysis and a numerical methods for solving of a functional-differential equations). M.-Izhevsk: RHD, 2004 (in Russian). = А. В. Ким, В. Г. Пименов. i-гладкий анализ и численные методы решения функционально-дифференциальных уравнений. М.-Ижевск: РХД, 2004. 332 Рис. 3: Максимум абсолютной погрешности для рассмотренного примера для различных значений шагов (Nx = 15, ..., 100, Ny = Nx , K = 15, ..., 100) Рис. 4: Максимум абсолютной погрешности для рассмотренного примера для различных значений шагов (Nx = 55, ..., 100, Ny = Nx , K = 55, ..., 100) 333 Investigation of the scheme of symmetrized derivatives for numerical solution of an advection equation with delay and 2-dimensional spatial variable Sergey V. Sviridov Ural Federal University (Yekaterinburg, Russia) Keywords: advection equation, delay, grid schemes, interpolation, extrapolation, stability, convergence order. We consider an equation in partial derivatives of the second order by space variable with delay. For such equation, based on the principle of separation of finite-dimensional and infinite-dimensional parts of state, we construct a grid method with symmetrized derivative. The piecewise linear two-dimensional and double interpolation and extrapolation by continuation are applied to the accounting of effect of heredity. It is proved that the considered method has order of a local error O(hx 2 + hy 2 + ∆2 ), where hx , hy are steps of discretization of the space variables, ∆ is a step of discretization of the time variable. Properties of double piecewise linear interpolation, conditions of stability and convergence are investigated by using the general theory of differential schemes. A test example is given. 334