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

№ 4 (32), 2014

Физико-математические науки. Математика

УДК 519.622

E. A. Новиков

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

Аннотация.

Актуальность и цели. При схемотехническом проектировании радиоэлектронных схем, моделировании электроэнергетических систем и других важных приложениях возникает необходимость решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений, не разрешенных относительно производной. Традиционные методы в основном ориентированы на решение разрешенных задач. Даже в простейшем случае приведение неявной системы к разрешенному виду связано с решением на каждом шаге линейной системы алгебраических уравнений. Матрица при производных, как правило, плохо обусловлена и часто является вырожденной, а разрешенная задача жесткая. Для ее решения требуется применение Z-устойчивых методов, в которых тоже требуется обращение матрицы. Эффективность расчетов можно повысить за счет одновременного разрешения системы и обеспечения Z-устойчивости численной схемы с применением одной матрицы.

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

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

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

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

E. A. Novikov

NUMERICAL MODELLING OF THE RING MODULATOR

BY THE METHOD FOR IMPLICIT SYSTEMS SOLUTION

Abstract.

Background. At schematic designing of radioelectronic circuits and other important applications there occurs a necessity to solve the Cauchy problem for stiff systems of ordinary differential equations, unsolved for derivatives. The known methods are mainly aimed at solving explicit problems. Even in a basic case, reduction of the implicit system to an explicit form is associated with solution of a linear system of algebraic equations at each step. The matrix at derivatives usually is poor conditioned and often degenerate, and the problem in an explicit form is stiff. For its solution one needs applying the L-stability methods, which also require decomposi-

1 Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 14-01-00047).

Physical and mathematical sciences. Mathematics

17

Известия высших учебных заведений. Поволжский регион

tion of the matrix. Efficiency of calculations can be increased by contemporary solving the system and meeting requirements of L-stability for a numerical scheme applying the same matrix.

Materials and methods. The decision function and its derivative were calculated approximately. For calculation accuracy control inequalities were applied. The first inequality provided accuracy of calculations and the second one was used for accuracy control of solution derivative calculations.

Results. The author developed an algorithm based on the L-stability method for solution of implicit problems. The method differs from the classic Rosenbrock schemes by approximate calculation of a solution derivative. The author constructed inequalities for stability control and gave the results of modelling the Ring Modulator.

Conclusion. The algorithm is designed for solving the Cauchy problem for implicit stiff systems of ordinary differential equations. The calculation results confirm the effectiveness of the constructed algorithm.

Key words: implicit system, Rosenbrock methods, accuracy control, ring modulator.

Введение

При моделировании динамических процессов в электрических цепях, электронных схемах и механических системах, при решении некоторых уравнений химической кинетики и других важных приложениях возникает проблема решения задачи Коши для жестких систем дифференциальных уравнений, не разрешенных относительно производной вида [1-6]

F(x,x,t) = 0, x(t0) = x0, to < t < tk , (1)

где x и F - вещественные At-мерные вектор-функции; t - независимая переменная. Современные численные методы обычно предполагают задание явной зависимости производной [1, 7-8], т.е.

x= f(t, x), x(to) = xo, to < t < tk . (2)

Приведение задачи (1) к виду (2) требует дополнительных затрат на шаг интегрирования, которые связаны с обращением матрицы Fy =dF (x, y, t) / dy , причем Fy часто вырожденная матрица. Разрешенная задача (2), как правило, жесткая. Это приводит к необходимости применения L -устойчивых методов, которые, в свою очередь, нуждаются в обращении матрицы Якоби df / dx . Возникает естественный вопрос об одновременном разрешении задачи (1) и обеспечении L -устойчивости численной схемы с целью повышения эффективности расчетов. Наиболее известные алгоритмы решения задачи (1) основаны на многошаговых численных формулах [1, 4]. Однако многие практические задачи описываются жесткими непрерывнодискретными системами, которые в современной литературе называют гибридными задачами [8]. Для них характерным является непрерывное описание режимов, смена которых происходит дискретно в некоторых точках. Такие точки задаются событийными функциями g(x), которые зависят от неизвестного решения x(t). В такой ситуации многошаговые методы могут

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

18

University proceedings. Volga region

№ 4 (32), 2014

Физико-математические науки. Математика

При решении жестких задач широкое распространение получили методы типа Розенброка [9] благодаря простоте реализации и достаточно хорошим свойствам точности и устойчивости. Здесь построен двухстадийный L-устойчивый метод второго порядка точности решения неявных задач. Метод отличается от классических схем типа Розенброка приближенным нахождением производной решения. Построены неравенства для контроля точности вычисления решения и его производной. Приведены результаты численного моделирования кольцевого модулятора.

1. Методы типа Розенброка для разрешенных задач

Для решения задачи (2) методы типа Розенброка записываются в виде

i-1

xn+1 = xn + Z Piki , Dnki = hf fn + cA xn + Z hkj

(3)

i=1

j=1

где ki, 1 < i < N, - стадии метода; h - шаг интегрирования; Dn = E - ah/П, E - единичная матрица; f„=df (tn, xn) / Эх - матрица Якоби системы (2), a , ci, Pi, Pj , 1 < i < N , 1 < j < i -1, - числовые коэффициенты, определяющие

свойства точности и устойчивости (3). В настоящее время методы типа Розенброка трактуются более широко - всякий метод, в котором матрица Якоби включена в численную схему, относится к схемам типа Розенброка или безы-терационным методам [1]. Для автономной задачи двухстадийный метод (3) имеет вид

xn+1 = xn + P1k1 + P2 k2 ,

Dnk1 = hf txn ) , Dnk2 = hf (n +P21k1) . (4)

При значениях коэффициентов

P1 = P21 =a = 1 -0,W2;, P2 = 1 - a = 0,s/2

численная схема (4) имеет второй порядок точности и является L -устойчивой. Для контроля точности численной формулы (4) можно применять неравенство [10-11]

||k2 - k4 <£,

где £ - требуемая точность интегрирования; Ц-Ц - некоторая норма в RN .

В случае большой размерности задачи (2) основные вычислительные затраты связаны с обращением матрицы Dn . При вычислении стадий ki , 1 < i < 2, обычно вместо обращения матрицы Dn дважды решается линейная система алгебраических уравнений с применением LU -разложения матрицы Dn с выбором главного элемента по строке или столбцу, а иногда и по всей матрице. Это приводит к тому, что при большой размерности задачи (2) общие вычислительные затраты фактически полностью определяются временем декомпозиции матрицы Dn .

Physical and mathematical sciences. Mathematics

19

Известия высших учебных заведений. Поволжский регион

2. Методы типа Розенброка для неявных систем

Используя обозначение x = у, задачу (1) можно переписать в виде задачи Коши для системы дифференциально-алгебраических уравнений вида

х = у, F(х, у, t) = 0, t0 < t < tk , (5)

с начальными условиями х (to ) = хо и у (to ) = уо . Дополнительное условие У (to ) = уо можно вычислить, например, решая задачу F (о, у, to ) = 0 методом установления [11]. Ниже будем предполагать существование и единственность решения задачи (1) или (5) [3]. Будем также предполагать, что функция F нужное число раз дифференцируема на шаге, а матрица

Dn = Fny + ahFnx

не вырожденная, где a - числовой коэффициент; h - шаг интегрирования. При записи Dn использованы обозначения

F = dF ixn,уп, tn ) f = dF (n,уп, tn )

дх ду

Тогда m -стадийный метод типа Розенброка для численного решения задачи (1) или (5) имеет следующий вид [5]:

хп+1 = хп + Z Piki , уп+1 = уп + Z Piki ,

i=1

i=1

Dnki = №пу

i—1

уп+Z

j=1

— ah Fnt —

—hF

i—1

i—1

i—1

+Z j, уп+Z j, tn +h • Z P/j

j=1

j=1

j=1

(6)

kf = —

ah

kf — h •

i—1

уп+Z htf

j=1

где a , Pi, $ij , 1 < i < N, 1 < j < i — 1, - числовые коэффициенты, Fnt =dF(хn,уп,tn)/ dt. Очевидно, что при применении численных формул

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

выполненным соотношение |^о||< Chp , где Fo = F(хо,уо,^); |(Ц - некоторая

nN /л

норма в R ; C - константа, не зависящая от величины шага интегрирования; p - порядок точности метода. Ниже будем предполагать невырожденность матрицы Fпу на каждом промежутке [tn, tn+1 ].

20

University proceedings. Volga region

№ 4 (32), 2014

Физико-математические науки. Математика

3. Двухстадийный метод второго порядка

Для решения (5) рассмотрим двухстадийную формулу вида

Xn+1 = xn + P\k\ + P2k2 , Уп+\ = Уп + Plk\ + Р2k2 ,

Dnk1 — h Fny ' yn - ahFnt - Fn (n, yn, tn ,

k{ — a-1 (x -hyn)/h, (7)

Dnk2 = hFny ■ (Уп + p21kiy ) - ahFnt -

-hFn (xn + P21k1X , yn + P21kiy , tn + P21h) ,

kiy = a 1

/ h .

Разложим стадии k1X, k1y, k2 и ty в ряды Тейлора по степеням h и

подставим в первую формулу (7). Сравнивая полученное представление приближенного решения с рядом Тейлора для точного решения, получим условия второго порядка точности схемы (7), т.е.

Р + Р2 = U Р21Р2 = 2 - a . (8)

Разлагая Fn+1 = F (xn+1, yn+1, tn+1) в ряд Тейлора и учитывая условия порядка (8), получим

Fn+1 -■

■2a + 0,5 F + 2а + (2a -1) P21 h 2 a2 n 4a

h 29(tn) + O (h 3)

где функция 9(t) не зависит от размера шага интегрирования. Требование Fn+1 - O(h2) приводит к уравнению

■2a + — - 0 . 2

(9)

Данное соотношение обеспечивает L -устойчивость схемы (7) при решении скалярного тестового уравнения х' — Хх, где X есть произвольное комплексное число, Re(Х)< 0 . Уравнение (9) имеет два вещественных корня

a1 — 1 - 0,W2 и a2 — 1 + 0,W2 . При решении разрешенных задач обычно выбирают корень a — a 1, потому что в этом случае меньше коэффициент в главном члене локальной ошибки. Как показывают многочисленные расчеты, при решении неявных систем выбор коэффициента a — a 1 тоже предпочтительнее.

Physical and mathematical sciences. Mathematics

21

Известия высших учебных заведений. Поволжский регион

В работе [12] показано, что требование L -устойчивости внутренних или промежуточных численных схем приводит к повышению эффективности расчетов и надежности вычислений. Требование L -устойчивости схемы

Хи+в = xn + Р2А

приводит к соотношению Р21 = а . Отсюда и условий порядка (8) получим коэффициенты L -устойчивого метода (7) второго порядка точности, т.е.

р1 =р21 =а = 1 -0,W2;, p2 = 1 - а = 0,W2. (10)

4. Контроль точности вычислений

Контроль точности вычисления решения схемы (7) можно построить по аналогии [5, 11], т.е. на каждом шаге следует проверять неравенство

jjkf-kf||<£, (11)

где k* и k2 заданы при описании численной схемы (7); £ - требуемая точность расчетов; Ц-Ц - некоторая норма RN .

В расчетах использовалась норма ||^n || вида

Ik

\tn

max

1<i<N I | хП | +r

Если

х

< r, то по i -й компоненте решения контролируется абсолютная

ошибка r - £; в противном случае контролируется относительная ошибка £ .

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

D~lF

1 П

<£ .

(12)

В расчетах использовалась норма ||^Щ вида

|| t,n ||= max | & | .

1<i< N

5. Результаты расчетов

В качестве примера выбрана модель кольцевого модулятора [13], схема которого приведена на рис. 1. Получая на входе низкочастотный сигнал Uin1 и высокочастотный сигнал Uin2, кольцевой модулятор генерирует на выходе смешанный сигнал U2 .

Каждый конденсатор, входящий в радиоэлектронную схему, приводит к появлению дифференциально уравнения CU' = I. Применение первого

22

University proceedings. Volga region

№ 4 (32), 2014

Физико-математические науки. Математика

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

CU[ = I1 - 0,513 + 0,5I4 +I7 - R-1U1, CU'2= I2 - 0,5I5 + 0,5I6 +I8 - R-1U2 , CSU'3= I3 -q(Udi) + q(UD4), CSU'4=-I4 + q(Um)-q(Um),

CsU5 = I5 + q (Udi )-q (Ud 3), CsU'6=-I6-q (Ud2 ) + q Ud4 ),

CpU7 =-R-U7 + q (Udi ) + q (Ud 2 )-q (Ud 3 )-q (Ud 4), где Udi , Ud2, Ud3 и Ud4 определены следующими соотношениями:

Udi = U3 -U5 -U7 -Uin2 , Ud2 = -U4 +U6 -U7 -Uin2 ,

Ud3 = U4 + U5 + U7 + Uin2, Ud4 =~U3-U6 +U7 +Un2.

Функция, определяющая поведение диода, задается в виде q (U ) = Y [ exp (SU)-1] , где у и 8 - константы.

Каждый индуктор также приводит к дифференциальному уравнению вида LI’ = U . Применение второго правила Кирхгофа к замкнутой цепи, содержащей индуктор, приводит к восьми дифференциальным уравнениям:

LhI{ = -Ui, LhI2 =-U2,

Ls2 I3 = 0,5Ui - U3 - Rg213, Ls3 Г4=-0,5Щ +U4 -Rg 314,

Ls 215 = 0,5U2 - U5 - Rg 215, Ls 316 =-0,5U2 +U6-Rg 3 Ib,

Ls117 = _U1 + Uin1 - (Ri + Rg1)17 ,

LslIs=~Ui-(Rc +Rgl) 18.

Physical and mathematical sciences. Mathematics

23

Известия высших учебных заведений. Поволжский регион

В начальный момент времени напряжения и силы токов отсутствуют, т.е. начальные условия по всем переменным равны нулю. Отождествляя напряжения с у-, 1 < i < 7, и силы токов с у-, 8 < i < 15, т.е. полагая

T

у = {Ui,...,U7,Ii,...,I8) , получим систему пятнадцати обыкновенных дифференциальных уравнений:

у1= C 1 ^ у8 - 2 у10 + 2 у11 + y14 - R 1 у1

у2 = C 1 [ у9 - 2 У12 + 2 У13 + У15 - R 1У2

у3 = Csl [ у10 - q (UD1 И q (UD4 , у4 = Cs 1 [-у11 + q (UD 2 )-q (UD3 ,

у5 = Csl [ у12 + q (UD1 )- q (UD3 , у6 = Csl [-у13 - q (UD2 ) + q (UD4 )\ ,

/ п-1

у7 = Ср

Rp у7 + q(UD1) + q(UD2 )- q(UD3 )- q(UD4 ) ,

y8 = ~Lhy1 , y9 = ~Lhy2 , y10 = Ls2 ^2У1 - y3 - Rg2y10 j , y11 = Ls3 ^ - 2 y1 + y4 - Rg3y11 j , y12 = Ls2 ^2 У2 р y5 - Rg2 y12 j ,

y13 = Ls3 ( -2y2 + y6 -Rg3y13 j, y14 = Ls1 -y1 + Uin1 (t)-{Ri + Rg1 )y14

У15 = L-1 [-У2-( + Rg1 )У15 ], У e R15, У(0) = 0, 0 < t < 10-3. (13)

Вспомогательные функции Ud1 , Ud2 , Ud3 и Ud4, q , Uin1 и Uin2 за-

даются формулами

UD1 = y3 - y5 - y7 - Uin2 (t) , UD2 =-y4 + y6 -y7 - Uin2 (t) ,

UD3 = y4 + y5 + y7 + Uin2 (t) , UD4 = -y3 -y6 + y7 +Uin2 (t),

q (U ) = y(e5U -1), Uin1 (t )= -2-sin (2000nt), Uin2 (t ) = 2sin (20000nt).

При вычислениях применялись параметры:

C = 1,6 • 10-8, Cp = 10“8, Cs = 10-12, Lh = 4,45, Ls1 = 0,002,

Ls2 = Ls3 = 5 -10-4, y = 40,67286402 U0-9, R = 25000,

Rp = R. = 50, Rg1 = 36,3, Rg 2 = Rg 3 = 17.3, Rc = 600, 5 = 17,7493332.

Расчеты проводились на Intel Core 2 Quad CPU с численной матрицей

-2 -2

Якоби с требуемой точностью е = 10 . Фактическая точность 10 достига-

24

University proceedings. Volga region

№ 4 (32), 2014

Физико-математические науки. Математика

—3

ется при задаваемой точности 10 . Время счета методами (4) и (7) практиче-

ски совпадает. На рис. 2 слева приведен выходной сигнал U2, а справа поведение решения для переменной U4 .

Рис. 2. Зависимость от времени U2 (а) и U4 (б)

Заключение

Максимальная эффективность построенного алгоритм достигается при небольшой точности расчетов - порядка 1 % и ниже, что связано с низким порядком точности применяемых численных формул. Задача (13) является достаточно серьезным испытанием для метода интегрирования вследствие большого разброса собственных чисел матрицы Якоби - коэффициент жест-12

кости порядка 10 . Кроме того, по ряду компонент наблюдаются сложные

колебания, которые систематически становятся высокочастотными (см. пове-

Physical and mathematical sciences. Mathematics

25

Известия высших учебных заведений. Поволжский регион

дение функции U4 на рис. 2,б). Это накладывает дополнительные требования

как на численную схему, так и на алгоритм управления величиной шага интегрирования.

Список литературы

1. Hairer, E. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems / E. Hairer, G. Wanner. - Berlin : Springer-Verlag, 1996. -614 p.

2. Brayton, R. K. A new efficient algorithm for solving differential-algebraic systems using implicit backward differentiation formulae / R. K. Brayton, F. G. Gustavson, G. D. Hachtel // Proceedings of IEEE, 1972. - P. 98-108.

3. Бояринцев, Ю. Е. Численные методы решения сингулярных систем / Ю. А. Бояринцев, В. А. Данилов, А. А. Логинов, В. Ф. Чистяков. - Новосибирск : Наука, 1989. - 223 с.

4. Gear, C. W. Differential-algebraic equations index transformations / C.W. Gear // SIAM J. Sci. Stat. Comput. - 1988. - Vol. 9, № 1. - P. 39-47.

5. Novikov, E. A. Some methods of solving ordinary differential equations unsolved for the derivative / E. A. Novikov, L. A. Yumatova // Soviet Math. Dokl. - 1988. -Vol. 36, № 1. - P. 138-141.

6. Левыкин, А. И. Класс (даД)-методов решения неявных систем / А. И. Левы-кин, Е. А. Новиков // Доклады РАН. - 1996. - Т. 348, № 4. - С. 442-445.

7. Dekker, K. Stability of Runge-Kutta methods for stiff nonlinear differential equations / K. Dekker, J. G. Verwer. - North-Holland, 1984. - 332 p.

8. Новиков, Е. А. Компьютерное моделирование жестких гибридных систем / Е. А. Новиков, Ю. В. Шорников. - Новосибирск : Наука, 2013. - 451 с.

9. Rosenbrock, H. H. Some general implicit processes for the numerical solution of differential equations / H. H. Rosenbrock // Computer. - 1963. - № 5. - P. 329-330.

10. Артемьев, С. С. Минимизация овражных функций численным методом для решения жестких систем уравнений [Препринт № 74] / С. С. Артемьев, Г. В. Демидов, Е. А. Новиков. - Новосибирск : ВЦ СО АН СССР, 1980. - 13 c.

11. Новиков, Е. А. Явные методы для жестких систем / Е. А. Новиков. - Новосибирск : Наука, 1997. - 194 с.

12. Новиков, Е. А. Замораживание матрицы Якоби в (3,2)-методе решения жестких систем / Е. А. Новиков, А. Л. Двинский // Вычислительные технологии. -2003. - Т. 8, № 3. - С. 272-278.

13. Mazzia, F. Test Set for Initial Value Problem Solvers / F. Mazzia, F. Iavernaro. -University of Bari : Department of Mathematics, 2003. - 295 p.

References

1. Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problem. Berlin: Springer-Verlag, 1996, 614 p.

2. Brayton R. K., Gustavson F. G., Hachtel G. D. Proceedings of IEEE, 1972, pp. 98-108.

3. Boyarintsev Yu. E., Danilov V. A., Loginov A. A., Chistyakov V. F. Chislennye meto-dy resheniya singulyarnykh system [Numerical methods for solving singular systems]. Novosibirsk: Nauka, 1989, 223 p.

4. Gear C. W. SIAM J. Sci. Stat. Comput. 1988, vol. 9, no. 1, pp. 39-47.

5. Novikov E. A., Yumatova L. A. Soviet Math. Dokl. 1988, vol. 36, no. 1, pp. 138-141.

6. Levykin A. I., Novikov E. A. Doklady RAN [Reports of the Russian Academy of Sciences]. 1996, vol. 348, no. 4, pp. 442-445.

7. Dekker K., Verwer J. G. Stability of Runge-Kutta methods for stiff nonlinear differential equations. North-Holland, 1984, 332 p.

26

University proceedings. Volga region

№ 4 (32), 2014

Физико-математические науки. Математика

8. Novikov E. A., Shomikov Yu. V. Komp’yuternoe modelirovanie zhestkikh gibridnykh sistem [Computer modeling of stiff hybrid systems]. Novosibirsk: Nauka, 2013, 451 p.

9. Rosenbrock H. H. Computer. 1963, no. 5, pp. 329-330.

10. Artem'ev S. S., Demidov G. V., Novikov E. A. Minimizatsiya ovrazhnykh funktsiy chislennym metodom dlya resheniya zhestkikh sistem uravneniy [Minimization of ravine functions by the numerical method for solving stiff systems of equations]. Preprint № 74. Novosibirsk: VTs SO AN SSSR, 1980, 13 p.

11. Novikov E. A. Yavnye metody dlya zhestkikh sistem [Explicit methods for stiff systems]. Novosibirsk: Nauka, 1997, 194 p.

12. Novikov E. A., Dvinskiy A. L. Vychislitel’nye tekhnologii [Computing technologies]. 2003, vol. 8, no. 3, pp. 272-278.

13. Mazzia F., Iavernaro F. Test Set for Initial Value Problem Solver. University of Bari: Department of Mathematics, 2003, 295 p.

Новиков Евгений Александрович

доктор физико-математических наук, главный научный сотрудник, Институт вычислительного моделирования Сибирского отделения Российской академии наук (Россия, Красноярск, Академгородок, 50, стр. 44)

E-mail: novikov@icm.krasn.ru

Novikov Evgeniy Aleksandrovich Doctor of physical and mathematical sciences, chief researcher, Institute of computational modeling of the Siberian Branch of the Russian Academy of Sciences (building 44, 50 Academgorodok campus, Krasnoyarsk, Russia)

УДК 519.622 Новиков, Е. А.

Численное моделирование кольцевого модулятора методом решения неявных систем / E. A. Новиков // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2014. - № 4 (32). -С.17-27.

Physical and mathematical sciences. Mathematics

27