Метод ускоренного численного расчета переходных процессов в электрических цепях на основе аппроксимации решения алгебраическими полиномами Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК [621.3.011+621.3.013]::519.6

Тиховод С. М.1, Корнус Т. М.2, Паталах Д. Г.3

1Д-р техн. наук, доцент, Запорожский национальный технический университет, Украина 2Ст. преп., Запорожский национальный технический университет, Украина 3Магистр, Запорожский национальный технический университет, Украина

МЕТОД УСКОРЕННОГО ЧИСЛЕННОГО РАСЧЕТА ПЕРЕХОДНЫХ ПРОЦЕССОВ В ЭЛЕКТРИЧЕСКИХ ЦЕПЯХ НА ОСНОВЕ АППРОКСИМАЦИИ РЕШЕНИЯ АЛГЕБРАИЧЕСКИМИ ПОЛИНОМАМИ

Цель работы — разработка метода численного расчета переходных процессов в электрических цепях, приводящая к сокращению процессорного времени расчета, а также разработка схемной модели метода, создающая удобство при его практическом использовании. Научная новизна. Разработанный метод отличается тем, что в процессе расчета отсутствует операция вычисления производной, что повышает точность расчета. Разработанная схемная модель позволяет от электрической цепи, в которой переходные процессы описываются интегро-дифференциальными уравнениями, перейти к цепи с изображениями токов, для которых справедливы законы Кирхгофа особой цепи постоянного тока. Изображение тока, изменяющегося во времени, это вектор, содержащий коэффициенты полиномиальной аппроксимации. Практическая ценность. Разработанный метод открывает возможность использования всего многообразного аппарата теории цепей для работы с изображениями токов, на основании чего может быть разработан универсальный программный комплекс для расчета переходных процессов в сложных электрических цепях. Методы исследования. Использована полиномиальная аппроксимация решения и матричные методы.

Ключевые слова: переходный процесс, численные методы, схемная модель, полиномиальная аппроксимация.

ПОСТАНОВКА ЗАДАЧИ

Переходные электромагнитные процессы в электротехнических системах, приводят к броскам тока при ком -мутациях. Это представляет значительную опасность для оборудования, а также влияет на надежность релейной защиты, поэтому исследование этих процессов актуально. Для компьютерного моделирования переходных процессов в электрических цепях в настоящее время используется ряд готовых программных комплексов, таких как ЕМТР [1], Р8рюе [2], 81ши1шк [3] и др. Уравнения состояния, составляемые автоматически в этих программных комплексах по законам Кирхгофа для мгновенных значений, являются интегро-дифференциально-алгебраичес-кими уравнениями, которые во многих случаях бывают жесткими. Обычно системы таких уравнений преобразуют в системы дифференциальных уравнений первого порядка, но при этом порядок системы уравнений возрастает. Для решения таких систем в программных ком -плексах широко применяются многошаговые методы численного интегрирования дифференциальных уравнений с использованием полиномиальной аппроксимации решения. В результате получены различные разностные схемы [4], позволяющие вычислить значение искомой функции в одной временной точке по известным значениям функции (или ее производных) в нескольких предыдущих точках.

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

ические петли». Наличие таких петель отрицательно сказывается на выполнении расчетов с помощью указанных программных комплексов. Поэтому разработка усовершенствованного метода расчета переходных электромагнитных процессов более быстродействующего и устойчивого к «алгебраическим петлям» является актуальной задачей.

Инженерам-электрикам важен физический смысл математических действий. Поэтому, если какая-либо математическая операция заменяется схемной моделью, то это повышает наглядность этой операции. Схемная модель должна позволять переход от электрической цепи, в которой процессы описываются интегро-дифферен-циальными уравнениями, к цепи постоянного тока с изображениями токов. Это открывает возможность использования всего многообразного аппарата теории цепей для работы с изображениями токов. Поэтому, модификация численного метода, сопровождающаяся созданием адекватной схемной модели для инженеров-электриков, является ценной.

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

ИЗЛОЖЕНИЕ МАТЕРИАЛА

Рассмотрим одноконтурную цепь переменного тока, содержащую резистивный (К), индуктивный (Ь) и емкостный (С) элементы, включенные последовательно. Пусть до коммутации конденсатор был заряжен до напряжения и (0). При подключении при (=0 источника ЭДС е(() в цепи происходит пе-

© Тиховод С. М., Корнус Т. М., Паталах Д. Г., 2015

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

di 1 * L-t + Ri + — J i(t )dt + uC (0) = e(t)

(1)

Проинтегрируем выражение (2) от нуля до к-й точки при изменении номера к от 1 до N-1:

tk tk

i(tk) = Jp(t)dt + io = J (ao + ait + «2*2 +-----+ ■

0 0

Будем искать решение во временной области т, в которой выделим N узловых точек /0, /1, 12,... , следующих с шагом к.

Решение для производной тока, как функцию от времени, в интервале времени , %_1] аппроксимируем полиномом N-1 -ой степени: сИЦ)

dt

p(t) = ao + ait + «2t2 +— aNtN 1.

(2)

Для аппроксимирующего полинома (2) зададим условие, что в точках ?деления интервала изменения аргумента значения производной по времени тока совпадают со значениями аппроксимирующего полинома:

{(гк) = р{гк) для к = 0, 1, 2, ... N-1. (3)

Если условие (3) записать для каждой точки / то получим систему линейных алгебраических уравнений, если принять, что ?= 0:

a0 = i ' (t0) = i0

ao + aih + a2h 2 +-----+ aN-ihN-1 = i '(h)

a0 + ai (2h) + a2 (2h)2 +-----+ aN — (2h)N-1 :

: i(2h)

a0 + ai ((N - 1)h) + a2((N - 1)h)2 + - + an_i ((N - 1)h)N_ : = i ((N _ 1)h)

(4)

Вычтем из уравнений системы (4) первое уравнение и получим сокращенную систему, которая в матричной форме имеет вид:

+ • • ■ + aN _1tN 1)dt + i0 =

1 2 1 3 1 N

= a0tk + -a1tk + Ja2tk + - + NaN_1tk + l0 , (9)

где tk=kh.

Подстановка в интегралы выражения (9) значений k от 1 до N-1 дает:

J p(t)dt + i0 = i0 h + f h 2 + f h3 + -

+ aN_1 hN + . . , + —h + «0 = «1,

J p(t)dt + i0 = i0 (2h) + + i3L(2h)3 + -

0 2 3

aN _1 (2h)

N

■+¡0 = г2;

J p(t)dt + i0 = ¿0 ( N _ 1)h + OÏL (( N _ 1)h)2 +

0

3

+ a2((N _ 1)h)~ + + aN_1 (( N 1)h)N+1 + i i

+-3-+-■-+~n~ ((N _1)h) + i0 = ¿n _1 '

(10)

h 2h

(2h)

(N _ 1)h ((N _ 1)h)2 i'(h) _ i'0 i' (2h) _ i0

i ((N _ 1)h) _ i0

(2h)N

((N _ 1)h)N

a1 a2

aN _1

или

V A = I'- I'

(5)

(5')

где V - матрица Вандермонда без первой строки и первого столбца; А = [а1 а 2... а^]1 - вектор коэффициентов аппроксимирующего полинома;

I' = [1 (к) I (2к). I ((N-1 )к)]т - вектор значений производных тока в точках 1, 2,.. .N-1.

Будем считать, что номер к отрезка, на которые разделен интервал изменения аргумента, совпадает с номером точки деления ?, расположенной слева отрезка.

В результате получим следующую систему в матричной форме:

¿n _1

(2 h)2 (2h)3

N (2h) N

N

1 ? 1 3

-(( N _ 1)h)2 -((N _ 1)h)3 .2 3

"~((N _ 1)h)N

N

a1 " h '

a2 2h

+ i0

.aN _1, (N _ 1)h

+10

(11)

или

2

N _1

h

h

2

3

h

h

2

3

и

I = 8 А + I' ® Н + 10 , (12)

где I - вектор значений тока в опорных точках;

т

10 = ['0> '0, ¡0'■■■ '0] ;

Н -[к 2к ••• (Ы - 1)к]; (13)

символ ® означает поэлементное умножение вектора на каждый столбец матрицы.

Займемся теперь вычислением напряжения на конденсаторе (третье и четвертое слагаемые в выражении (1)).

UC = BWA + B

i0 ® н+2 i'0 ® н2

+ и

C 0, (16)

где обратная величина емкости конденсатора В - —;

W

h3 2 ■ 3 (2 h )3 2 ■ 3

h 4 3 ■ 4 (2h)4

3 ■ 4

— ((N - 1)h)3 ——((N - 1)h)4

2 ■ 3 3 ■ 4

C

h'

N ■ (N + 1)

(2h)

N+1

N ■ (N + 1)

1

N ■ (N + 1)

X ((N - 1)h)N

(17)

1 \

UC (t) =77Í i(t )dt + uc (0) =

C í

С 12 13

J (a0t + -+ -a2t + +

+ "' + "N aN-1tN + 'o)dt

+ uc (0).

(14)

Вычислим интеграл в выражении (14) и распишем это выражение для точек ( = (к = кк, к = 1, 2, „., N-1. Получим в матричной форме выражение:

UC1 UC 2

UCN-1

h3 h 4 hN+1

2 ■ 3 3 ■ 4 N ■ (N +1)

(2h)3 (2h)4 (2h)N+1

2 ■ 3 3 ■ 4 N ■ (N +1)

— ((N - 1)h)3 —((N - 1)h)4 2 ■ 3 3 ■ 4

aj a2

1 ((N - 1)h) N+1

N ■ (N +1)

aN-1

h 2h

(N - 1)h

+ i0

(2h)2

((N - 1)h)2

+ Uc

(15)

В компактном виде выражение (15) принимает вид:

H 2

(2h)2

((N - 1)h)2

(18)

иС0 = [иС0. иС0, ■ иС0]Т ; иС0 = ис(0).

В уравнение (1) для опорных точек подставим аппроксимацию всех слагаемых согласно выражений (5'), (12), (16) и получим матричное уравнение:

LVA + LI'0 + RSA + R(I0 +10 ® H) + B

WA +

+i0 ® н+-2 i'0 ® н2

+ UC0=e.

(19)

Преобразуем уравнение (19):

(LV + +RS + BW)A = e - LI'0 - R(I0 +1'0 ® H) -.

- B

i0 ® н+2 i'0 ® н2

- U

C 0,

(20)

где е - вектор значений ЭДС источника в точках 1, 2, ..., N-1 временного интервала;

Уравнение (20) можно интерпретировать следующим образом. Пусть в исходной ветви К-Ь-С протекает ток 1ф. Тогда согласно уравнению (20) исходной ветви соответствует ветвь замещения, в которой проходит операторный ток А, изображающий исходный ток При этом в ветви замещения резистивный элемент имеет операторное сопротивление Я8 и последовательно с ним навстречу току включается источник ЭДС

К(10 +10 ® Н) (см. рис. 1).

X

2

h

0

+

2

h

+

Рисунок 1 - Интерпретация уравнения (20): Ж, ¿V, BW - операторные сопротивления резистивного, индуктивного и

емкостного элементов;

Що +1 0 ®Н), ы'о, в

10 ® н+-21'0 ® н2

+ и

С о - дополнительные источники ЭДС; е -

вектор значений ЭДС в узловых точках.

Индуктивный элемент имеет операторное сопротивление ¿V и последовательно с ним навстречу току включается источник ЭДС ¿I 0, а емкостный элемент имеет операторное сопротивление - BW и последовательно с ним навстречу току включается источник ЭДС

в

10 ® н+-21'0 ® н2

+ и

С 0.

Докажем, что в узлах схемы замещения для изображений А соблюдается закон токов Кирхгофа. В любом узле электрической цепи для токов ветвей, принадлежащих узлу, в любой момент времени а, следовательно, токов в начале интервала ¡0 и для векторов токов I выполняется закон токов Кирхгофа:

I (I к - 10к) = 0, к=1

(21)

где Ь - количество ветвей, сходящихся к узлу, к - текущий номер ветви, сходящейся к узлу.

Продифференцируем уравнение (21) и получим:

I (I к - ^к) = 0. к=1

(22)

Согласно уравнению (5'): V А = I'- Г0 следовательно

I VAk=0.

к=1

(23)

Если уравнение (23) умножить на матрицу, обратную матрице V, то получим:

Ь

I Ак = 0.

к=1

(24)

Из изложенного материала сделаем выводы. Реальному току ¡(1) соответствует векторное изображение А в схеме замещения, показанной на рис. 1. Все изображения тока А в схеме замещения удовлетворяют законам Кирхгофа, если схема замещения составляется по правилам:

- источник ЭДС заменяется векторным источником е, содержащим значения ЭДС в N-1 опорной точке;

- резистивный элемент имеет операторное сопротивление ЯД и последовательно с ним навстречу току

включается источник ЭДС +1 0 ® Н);

- индуктивный элемент имеет операторное сопротивление ¿V и последовательно с ним навстречу току

включается источник ЭДС ¿I 0;

- емкостный элемент имеет операторное сопротивление BW и последовательно с ним навстречу току включается источник ЭДС

в

!0 ® Н + 2 !0 ® н2

+и

С 0.

Таким образом, в схеме замещения электрической цепи изображения А оригиналов токов ¡ (¿) удовлетво -ряют законам Кирхгофа. Следовательно, при известных значениях токов ветвей / и напряжений на конденсаторах пск в точке ¿0 начала интервала т система уравнений, составленная по законам Кирхгофа, для всех узлов без одного и для всех главных контуров имеет единственное решение. В результате решения системы линейных алгебраических уравнений получаем векторы полиномиальных коэффициентов А для всех ветвей. Зная для любой ветви коэффициенты полинома и значение /в начальной точке ¿0, мы можем получить значение тока и напряжения на конденсаторе во всех произвольных точках любого из N отрезков в интервале времени т .

Согласно [5] с ростом степени полинома погрешность интерполяции уменьшается, однако увеличивать степень полинома бесконечно нельзя. С ростом числа N матрицы становятся плохо обусловленными. Чтобы при малом шаге к матрицы V, W и 8 не стали плохо обусловленными необходимо шаг интегрирования умножить на нормирующий коэффициент такой, чтобы нормированный шаг был близок к единице. Тогда на нормирующий коэффициент нужно умножить значения всех индуктив-ностей и емкостей, входящих в схему, а значение частоты нужно разделить на нормирующий коэффициент. Это улучшит обусловленность матриц, но увеличение степени полинома больше десяти может в некоторых случаях все-таки привести к неадекватному решению. Погрешность интерполяции можно существенно снизить,

если выбирать положения опорных точек не равномерно, а по методу, разработанному Чебышёвым. В этом случае максимальная погрешность будет минимизирована. Если аргумент ( аппроксимируемой функции задан на отрезке [а, Ь], то опорные точки ( выбираются в нулях полиномов Чебышёва [6]:

tk =

где xk = - cos

2k +1

2 N

т Xk + a + b 2 :

, k=0,1,... N-1.

(25)

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

в следующем цикле. Если же

ные значения иС0 , —(0)

' Ж

достаточно значений тока только в узловых точках, то их можно вычислить не по общей формуле (9), а проще, согласно (12):

I = 8 А + 1' ® Н + 10 .

Значения напряжения на конденсаторе в узловых точках вычисляются согласно (16):

UC - BWA + B

i0 ® h+-2-i0 ® h2

+u,

C 0

Для испытания разработанной методики в системе МаНаЬ составлено несколько компьютерных программ для расчета переходных процессов в электрических цепях.

Приведем пример одного из расчетов. Исследуем переходный процесс в модельной одноконтурной цепи, содержащей резистивный (К), индуктивный (Ь) и емкостный (С) элементы, включенные последовательно. Рассчитаем процесс изменения тока в цепи после подключения источника постоянной ЭДС. Согласно изложенному правилу, составим схему замещения для изображения тока (рис. 1).

По программе, составленной согласно предложенной методике, выполнен расчет переходного процесса при следующих значениях исходных данных: е(()=Е= 100 В; К=0,7 Ом; Ь=0,005 Гн; С=500 мкФ; ¡(0)=0;

— Е

— (0) --; иС (0) - 0 .

— Ь

Алгоритм программы заключается в циклическом решении матричного уравнения:

А - (ЬУ + +К8 + BW)-1 • ] е - Ы0 - К(!0 +1'!) ® Н) -

- B

i0 ® h+-1'0 ® н2

- и

C0

(27)

и вычислению тока в опорных точках согласно (12) , напряжения на конденсаторе согласно (16) и производной тока согласно формуле: I' = У А + Г0 .

График тока ¡(() , полученный в результате расчета, представлен на рис. 2.

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

Полученное аналитическое выражение для тока 1(1) имеет вид:

i(t) - A1ePlt + A2eP2t,

(28)

где A1 = -15.81277i; A2 = 15,81277i; p1= -70 + 632,4002i; p2= -10 - 632,4002i.

Точки, соответствующие точному аналитическому выражению (28), на графике рис. 2 показаны звездочками. Шаг интегрирования при расчете выбран таким, чтобы отклонение значения тока от соответствующего точного значения в точках локальных максимумов не превышало 0,1% (при уменьшении шага интегрирования погрешность уменьшается). Выполнен также расчет модельной задачи при использовании метода Гира с максимальным постоянным шагом интегрирования таким, чтобы отклонение значения тока от соответствующего точного значения в точках локальных максимумов также не превышало 0,1%. С помощью операторов tic/ toc оценивалось время расчета. Сравнение процессорного времени расчета модельной задачи по предложенному методу и по методу Гира показало, что предложенный метод имеет быстродействие более чем в четыре раза лучшее, чем многошаговый метод Гира.

ВЫВОДЫ

1. Предложенная методика построения схем замещения электрической цепи позволяет непосредственно вычислять коэффициенты полиномиальной аппроксимации искомых токов в переходных режимах.

2. Векторы полиномиальных коэффициентов удов -летворяют законам токов и напряжений Кирхгофа для схемы замещения и являются решениями системы алгебраических уравнений.

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

4. Предложенная методика отличается тем, что в ней отсутствует операция вычисления производной решения, что позволяет значительно повысить точность расчетов.

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

71

0 10 20 30 40 50 60

тБ

Рисунок 2 - Временная зависимость тока ¡(1). Звездочками на графике показаны точные значения, вычисленные аналитически.

A

СПИСОК ЛИТЕРАТУРЫ

1. Electromagnetic Transient Program (EMTP) Application Guide // EPRI Report No: EL - 4650, Project 2149-1, Westinghouse Electric Corp., Pittsburgh, PA, 1986.

2. Кеоун Д. OrCAD Pspice. Анализ электрических цепей. / Дж. Кеоун. - сПб, : Питер, 2008. - 640 с.

3. Черных И. В. Simulink среда создания инженерных приложений. / И. В Черных. - М. : ДИАЛОГ-МИФИ, 2003. - 496 с.

4. Современные численне методы решения обыкновенных дифференциальных уравнений./ редакторы Дж. Холл, Дж. Уатт.-М.: Мир, 1979.- 312 с.

5. Shampine L. F. Fundamentals of numerical computing./ L. F. Shampine, R. C. Allen, S. Pruess.-New York Chichester Brisbane Toronto Singapore.: JOHN WILEY & SONS, INC., 1997. - 268 p.

6. Ильина В. А. Численные методы для физиков-теоретиков. 1. / В. А. Ильина, П. К. Силаев. - Москва-Ижевск: Институт компьютерных исследований. -2003. -132 с.

Статья поступила в редакцию 2.10.2015

Тиховод С. М.1, Корнус Т. М.2, Паталах Д. Г.3

'Д-р техн. наук, доцент, Запорiзький нацюнальний техшчний ушверситет, Украша

2Ст. викл., Запорiзький нацюнальний техшчний ушверситет, Украша

3Магистр, Запорiзький нацюнальний техшчний ушверситет, Украша

МЕТОД ПРИСКОРЕНОГО ЧИСЛОВОГО РОЗРАХУНКУ ПЕРЕХ1ДНИХ ПРОЦЕС1В В ЕЛЕКТРИЧ-НИХ КОЛАХ НА ОСНОВ1 АПРОКСИМАЦП РОЗВ'ЯЗКУ АЛГЕБРА1ЧНИМИ ПОЛ1НОМАМИ

Мета роботи — розробка методу числовогорозрахунку перехiдних проце^в в електричних колах, що призво-дить до скорочення процесорного часу розрахунку, а також розробка схемног моделi методу, що створюв зручтсть при його практичному використант. Наукова новизна. Розроблений метод вiдрiзняeться тим, що в процесрозрахунку вiдсутня операщя обчислення похiдно'i, що тдвищув точтсть розрахунку. Розроблена схемна модель дозволяв вiд електричного кола, вякому перехiднi процеси описуються ттегро-диференщальнимирiвнян-нями, перейти до кола 1з зображеннями струмiв, для яких справедливi закони Кирхгофа особливого кола пост-шного струму. Зображення струму, що змтюеться в чап, це вектор, що метить коефщвнти полiномiальноi апроксимаци.

Практична цттсть. Розроблений метод вiдкривавможлив^ть використання всьогорiзноманiття апара-ту теори кы для роботи з зображеннями струмiв, на пiдставi чого може бути розроблений утверсальний програмний комплекс для розрахунку перехiдних проце^в у складних електричних колах. Методи до^дження. Використана полiномiальна апроксимащя ршення та матричт методи.

Ключов1 слова: перехiдний процес, числовi методи, схемна модель, полiномiальна апроксимащя.

Tykhovod S.1, Kornus T.2, Patalakh D.3

'Doctor of science, Assoc. Prof., Zaporozhye national technical university, Ukraine 2Senior Lecturer, Zaporozhye national technical university, Ukraine 3Magistr, Zaporozhye national technical university, Ukraine

METHOD OF ACCELERATED NUMERICAL CALCULATION OF TRANSIENTS IN ELECTRICAL CIRCUITS BASED ON SOLUTION APPROXIMATION BY ALGEBRAIC POLYNOMIALS

The work purpose is modification of numerical calculation of transients in electric circuits that reduces the CPU time of calculation, and also development of circuit model method, that creates the convenience in its practical use. The circuit model of the modified numerical calculation of transients in electric circuits is developed. The circuit model created for the electric circuit, in which processes are described with the help of integro-differential equations, turn into the circuit with current images, for which Kirchhoff's rules, reducing to the algebraic equations, are valid. Time-varying current image is vector, containing the polynomial approximation coefficients. Practical importance. The developed method uncloses the possibility of use of the whole manifold of the circuit theory to work on current images. On the basis of that the universal software package for transient analysis in the complex electric circuits could be designed. Research methods. Polynomial approximation of solution and matrix calculations is used.

Key words: transient, numerical methods, circuit model, polynomial approximation

REFERENCES

1. Electromagnetic Transient Program (EMTP) Application Guide. EPRI Report No: EL. 4650, Project 2149-1, Westinghouse Electric Corp., Pittsburgh, PA, 1986.

2. Keown John. OrCAD Pspice and Circuit Analysis. John. Keown. sPb,: Piter, 2008, 640 c.

3. Chernyh I. V. Simulink environment of making of engineering applications. I. V. Chernyh, M, DIALOG-MIFI, 2003, 496 p.

4. Modern Numerical Methods for Ordinary Differential Equation, Edited by G. Hall and J. M. Watt, Clarendon Press Oxford 1976, 312 p.

5. Shampine L. F. Fundamentals of numerical computing, L. F. Shampine, R. C. Allen, S. Pruess, New York Chichester Brisbane Toronto Singapore, JOHN WILEY & SONS, INC., 1997, 268 p.

6. Il'ina V. A. Numerical methods for physicist-theorists. V. A. Il'ina, P. K. Silaev, Moskow-Izhevsk: Institut komp'juternyh issledovanij, 2003, 132 p.