авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 |

«ПЕРМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ УДК 517.97+517.688 На правах рукописи ...»

-- [ Страница 2 ] --

( ) l= Таким образом, оценка погрешности также, как и в случае = [a, b], свелась к простому (хотя и вычислительно дорогому) суммированию произведений.

Глава Программная реализация Описываются основные особенности программной реализации доказательного вычислительного эксперимента. Текст программы приведён в Приложении A.

4.1 Структура программы Исходный текст программы Fredholm состоит из трёх частей:

• вспомогательные процедуры (численное интегрирование, исследование системы линейных алгебраических уравнений, работа с многочленами и т.д.);

• постановка задачи;

• доказательный вычислительный эксперимент.

Бльшую часть программы составляют вспомогательные процедуры. Ис o пользование объектно–ориентированного программирования посзволяет суще ственно упростить собственно рабочую часть программы.

Для компьютерной реализации был выбран язык С++ (использовался на бор компиляторов GNU GCC 4.3.2 из Linux Fedora 10 x86_64).

4.2 Вспомогательные процедуры 4.2.1 Численное интегрирование Для оценки значения интеграла b1 bp... f (t1,..., tp ) dt1... dtp a1 ap применяется произведение одномерных формул Гаусса–Кронрода.

DOUBLE i n t e g r a l (FUNC f, const s t d : : v a l a r r a y DOUBLE& lb, const s t d : : v a l a r r a y DOUBLE& ub, DOUBLE& e r r o r, PARAMETERS p = NULL, s i z e _ t max_div = 1 0 0 0 ) Тип FUNC определён следующим образом:

typedef DOUBLE ( FUNC) ( const s t d : : v a l a r r a y DOUBLE& x, PARAMETERS p ) Здесь DOUBLE число с плавающей точкой (double или long double, в зависи мости от используемого компилятора (Microsoft Visual C++ или GNU GCC);

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

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

Определение числа процессоров (ядер) выполняется встраиваемой (inline) процедурой number_of_processors, определённой в файле datautil.hpp.

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

4.2.2 Многочлены и ортогональные функции Многочлены В шаблоне класса POLYNOM template c l a s s T c l a s s POLYNOM определяются параметризованные по типу многочлены с одной переменной вида c0 + c1 x + c2 x2 + · · · + cn xn.

Реализованы: вычисление значения многочлена (метод Горнера), арифмети ческие операции над многочленами, дифференцирование и интегрирование, построение последовательности Штурма и оценка числа корней на интервале (без учёта кратности).

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

Ортогональные многочлены Лежандра Класс LEGENDRE является потомком классов ORTHO и POLYNOMDOUBLE. В кон структоре первые шесть многочленов (со степенями от нуля до пяти) определя ются явно. Коэффициенты многочленов более высоких степеней вычисляются вычисляются по формуле (начиная с i = 5):

L0 (x) = 1, L1 (x) = x, 2i + 1 i Li+1 (x) = xLi (x) Li1 (x), i = 1, 2,...

i+1 i+ Сделаны отдельные реализации для многочленов Лежандра со значениями коэффициентов многочленов в виде рациональных чисел c l a s s LEGENDRE : public ORTHOmpq_class, public POLYNOMmpq_class и с коэффициентами в виде вещественных чисел двойной точности c l a s s LEGENDRE_DBL : public ORTHODOUBLE, public POLYNOMDOUBLE Это позволяет избежать ошибок округления при преобразованиях типов F Q.

4.2.3 Коэффициенты ряда Фурье Для оценки значения коэффициента fi1,...,ip,j1,...,jp ряда Фурье при разложении по заданной системе ортогональных многочленов mp p m f (x) ··· fi1,...,ip,j1,...,jp k ik (xk ), j1 =0 jp =0 k= ai xi bi.

применяется функция DOUBLE f o u r i e r _ c o e f f (DOUBLE ( f ) ( const s t d : : v a l a r r a y DOUBLE&, MATH: : PARAMETERS p ), const s t d : : v a l a r r a y DOUBLE& lb, const s t d : : v a l a r r a y DOUBLE& ub, & ortho, const s t d : : v e c t o r ORTHO DOUBLE e r r o r ) Здесь f раскладываемая в ряд функция, lb и lb нижняя и верхняя гра ницы области ((a1,..., ap ) и (b1,..., bp )), на которой выполняется разложение, вектор ортогональных функций.

ortho 4.2.4 Линейная алгебра Параметризованный класс MATRIX определён в файле linalg.hpp. В нём реали зованы арифметические операции с матрицами, процедуры вычисления следа и характеристического многочлена, а также решения системы линейных ал гебраических элементов и обращения матрицы методом Гаусса с частичным выбором главного элемента.

Для ускорения вычислений на многопроцессорных/многоядерных компью терах с общей памятью используется библиотека OpenMP. Распараллелены вы числения в процедурах перемножения матриц, а также в процедуре решения системы линейных алгебраических уравнений методом Гаусса–Жордана.

4.3 Постановка и решение задачи Общие для разных задач структуры данных определяются в заголовочном файле problem.hpp. Задание значения PROBLEM позволяет определить в файле problem.cpp пространство имён, соответствующее решаемой задаче.

В файле problem.hpp также определены структуры данных, которые ис пользуются в файле fredholm.cpp, где содержаться собственно процедуры ре шения интегрального уравнения Фредгольма второго рода.

Глава Примеры 5.1 Простейшая задача Лагранжа Для начала рассмотрим применение конструктивного подхода к решению про стейшей вариационной задачи Лагранжа:

b x(t)2 + p(t)x(t)2 + q(t)x(t) dt min, (5.1a) I(x) = a (5.1b) x(a) = 1, x(b) = Её выбор обусловлен тем, что она позволяет получить классическое реше ние, которое можно использовать для сравнения и верификации предлагаемого метода.

5.1.1 Классическое решение Необходимым условием существования решения задачи (5.1a)–(5.1b) является разрешимость уравнения Эйлера:

x (t) = p(t)x(t) + q(t), x(a) = 1, x(b) = Пусть a = 0, 1 = 1, b = 10, 2 = 1, и p(t) = sin(2t) exp(t/5), q(t) = 1 + cos(t).

Тогда d2 x(t) 1 + cos(t) = sin(2t) exp(t/5)x(t) +, dt x(0) = 1, x(10) = Для численного решения краевой задачи применим метод пристрелки, в котором задачу Коши будем решать методом Рунге–Кутта 4-го порядка. Для решения уравнения x 10 | x (a) = (определение значения производной в начале отрезка x (a)) будем использо вать метод деления отрезка пополам. Реализация для математического пакета GNU Octave приведена в Приложении B.

x (0) 2.6620.

- - - - - 0 2 4 6 8 5.1.2 Конструктивный подход Будем искать решение в классе абсолютно непрерывных функций (в простран стве X L2 [a, b] R). Следовательно, t (5.2) x(t) = 1 + z(s) ds.

a Подстановка (5.2) в (5.1a)–(5.1b) даёт b t I1 (z) = z(t) + p(t) 1 + z(s) ds + a a t b t + q(t) 1 + z(s) ds dt = z(t) + p(t) z(s) ds + a a a t z(s) ds + p(t)2 + q(t)1 dt.

+ 21 p(t) + q(t) a Переходя к сопряжённым интегральным операторам b t b b p(t) z(s) ds dt = z(t) p(s) ds dt, a a a t b t t p(t) z(s1 ) ds1 z(s2 ) ds2 dt = a a a b b b = z(t) p(s2 ) ds2 z(s1 ) ds1 dt, a a max(t,s) получаем:

b b z(t)2 + z(t) I1 (z) = 21 p(s) + q(s) ds + 1 p(t)1 + q(t) + a t b b p( ) d z(s) ds dt min, (5.3a) + z(t) a max(t,s) b (5.3b) 1 2 + z(t) dt = 0.

a Запишем функцию Лагранжа задачи (5.3a)–(5.3b):

b b z(t)2 + z(t) L(z, ) = 21 p(s) + q(s) ds + 1 p(t)1 + q(t) + a t b b b + z(t) p( ) d z(s) ds dt + 1 2 + z(t) dt.

a max(t,s) a Решение задачи z должно быть стационарной точкой функции Лагранжа, т.е.

b b b z Lz (, )[h] 2(t) + z p( ) d z (s) ds+ a a max(t,s) b + 21 p(s) + q(s) ds + h(t) dt = 0, t и, следовательно, корнем уравнения b z(t) K(t, s)z(s) ds = f (t), a где b b p( ) d и f (t) = K(t, s) = 1 p(s) + q(s) ds.

2 max(t,s) t Пусть, как и в предыдущем пункте, a = 0, 1 = 1, b = 10, 2 = 1, p(t) = sin(2t) exp(t/5), q(t) = 1 + cos(t).

def def Тогда, используя первообразные P (t) = p(t) dt и Q(t) = q(t) dt, функции = = ядра и правой части можно записать в виде (с учётом 1 = 1) K(t, s) = P max(t, s) P (b), Q(t) Q(b) f (t) = + 1 P (t) P (b) +, 2 где 5 t/ P (t) = e 10 cos(2t) + sin(2t), Q(t) = t + sin(t).

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

0. 0. 0. 0. -0. -0. -0. -0.4 0 2 10 0 2 4 6 8 Рис. 5.1: K(t, s) 0. 0. 0. -0. -0. 0 2 10 0 2 4 6 8 Рис. 5.2: K(t, s) K(t, s) Зададим на [a, b] сетку, узлы которой расположены с постоянным шагом, равным 2/3 (15 подынтервалов). В качестве базиса аппроксимации будем ис пользовать многочлены Лежандра до третьей степени включительно (аппрок симация f и K кубическими и бикубическими многочленами соответственно).

Аппроксимация ядра K(t, s) даёт кусочно–бикубический многочлен K(t, s), мало отличающийся от исходной функции. График поверхности функции по грешности K(t, s) = K(t, s)K(t, s) и проекция её линий уровня на плоскость (t, s) приведены на рис. 5.2. Видно, что погрешность аппроксимации больше там, где производная функции K(t, s) имеет разрыв, то есть в областях, лежа щих на диагонали t = s.

Графики функции правой части интегрального уравнения f (t) и погреш ности аппроксимации f (t) f (t) показаны рис. 5.3.

3 0. 2. 0. 0. 1. 0. -0. -0. -0. - -0. 0 2 4 6 8 10 0 2 4 6 8 Рис. 5.3: f (t) и f (t) f (t) После оценки коэффициентов Ki1 i2 j1 j2 и fij были вычислены оценки сверху для норм K K 0. L (L2 [a,b]) и f f 0.0012173432348182174419.

L2 [a,b] Так как отрезок [a, b] разбивается на 15 подотрезков, на каждом из кото рых применяется аппроксимация кубическими многочленами, то система ли нейных алгебраических уравнений (СЛАУ), к которой сводится интегральное уравнение с вырожденным ядром, будет иметь N = 60 уравнений и неизвест ных. Заметим, что если бы в качестве базиса аппроксимации использовались не ортогональные, а какие-либо другие функции (например, 1, t, t2, t3 ), то раз мерность СЛАУ увеличилась бы до 240 (15 · 4 · 4).

Оценка нормы резольвентного оператора R [8.1507713096520113538;

8.1507713407638782144].

L (L2 [a,b]) Так как оператор I K обратим и норма K K меньше нижней L (L2 [a,b]) границы интервальной оценки [0.10928040519877344888;

0.10928040557031780312] 1/(1 + R L (L2 [a,b]) ), то оператор I K также имеет обратный.

Оценка верхней границы спектра оператора K даёт max (K) 0.89020870876332758304, что меньше, чем 1 K K L (L2 [a,b]). Следовательно, экстремаль x доставляет функционалу минимум.

После решения СЛАУ с правой частью, зависящей от параметра, было восстановлено решение исходной задачи x(t). Значением параметра, приводя щим к выполнению второго краевого условия x(10) = 1, будет [5.9163981396704912186;

5.9163981378078460693].

На рис. 5.4 приведены графики функций z (t) = x (t) и x(t).

Норма производной решения x [7.3403467429287756474;

7.340346743943793939], L2 [a,b] и её погрешность x x 0.030268743295793779069.

L2 [a,b] - - - - - - - -4 - 0 2 4 6 8 10 0 2 4 6 8 Рис. 5.4: z (t) и x(t) Вывод программы:

Аппроксимация f... Ok Аппроксимация K... Ok Создание матрицы A... Ok Создание матрицы B... Ok Обращение матрицы A... Ok Решение СЛАУ... Ok Кусочно полиномиальное представление решения... Ok Вычисление множителей Лагранжа... Ok lambda = [-5.9163981396704912186;

-5.9163981378078460693] Собираем приближённое решение интегрального уравнения... Ok ||tilde z|| = [7.3403467429287756474;

7.340346743943793939] Собираем приближённое решение вариационной задачи... Ok ||tilde x|| = [20.677625794745821963;

20.677625797628390814] Оценка норм приближений... Ok ||tilde F|| = [5.2448921418123184424;

5.2448921425225893955] ||tilde K|| = [1.2210819144965860961;

1.2210819146872620156] Генерируем данные для графиков... Ok Оцениваем погрешность аппроксимации f...Ok max DeltaF = 0. ||DeltaF|| = 0. Оцениваем погрешность аппроксимации K...Ok max DeltaK = 0. ||DeltaK|| = 0. Оцениваем квадрат нормы резольвентного оператора... Ok ||R|| = [8.1507713096520113538;

8.1507713407638782144] 1/(1+||R||) = [0.10928040519877344888;

0.10928040557031780312] ||z - tilde z|| = 0. Вычисление характеристического многочлена для E-A... Ok Все корни внутри круга радиуса 1. Оцениваем верхнюю границу спектра... Ok max sigma = 0. Успешное завершение Затрачено 251.73сек.

5.2 Задача с сосредоточенным отклонением аргумента Усложним предыдущую задачу, добавив слагаемое с сосредоточенными откло нениями аргумента:

b x(t)2 + p(t)x(t)2 + q(t)x(t) + g(t)x h1 (t) x h2 (t) I(x) = dt min, a (5.4a) (5.4b) x() = 1 (), a, x() = 2 (), b, (5.4c) x(a) = 1, x(b) = 2.

Снова будем искать решение в классе абсолютно непрерывных функций (X L2 [a, b] R):

t (5.5) x(t) = 1 + z(s) ds.

a Введём обозначения def x h(t) = h (t) + Zh (t), (5.6) = где h (t) = 1 h(t) (,a) h(t) + 1 [a,b] h(t) + 2 h(t) (b,) h(t) и b Zh (t) = [a,b] h(t) [a,h(t)] (s)z(s) ds.

a Пусть b def x(t) = 0 (t) + Z0 (t), 0 (t) = 1, Z0 (t) = = [a,t] (s)z(s) ds.

a Подстановка (5.5) и (5.6) в (5.4a)–(5.4c) даёт b z(t)2 + p(t) 0 (t) + Z0 (t) + q(t) 0 (t) + Z0 (t) + I1 (z) = a + g(t) h1 (t) + Zh1 (t) h2 (t) + Zh2 (t) dt min, (5.7a) b (5.7b) 1 2 + z(t) dt = 0.

a Раскрыв в (5.7a) скобки и выполнив следующие преобразования b t b b 2p(t)1 + q(t) z(s) ds dt = z(s) 2p(t)1 + q(t) dt ds = a a a s b = 21 P (b) P (t) + Q(b) Q(t) z(t) dt, a f0 (t) b t t p(t) z(s1 ) ds1 z(s2 ) ds2 dt = a a a b b b = z(s1 ) p(t)[a,t] (s1 )[a,t] (s2 ) dt z(s2 ) ds1 ds2 = a a a b b = z(t) P (b) P max(t, s) z(s) ds dt, a a K0 (t,s) b h2 (t) g(t)h1 (t) [a,b] h2 (t) z(s) ds dt = a a b b g(s)h1 (s)[a,b] h2 (s) [a,h2 (s)] (t) ds z(t) dt, = a a f1 (t) b h1 (t) h g(t) (t) [a,b] h1 (t) z(s) ds dt = a a b b g(s)h2 (s)[a,b] h1 (s) [a,h1 (s)] (t) ds z(t) dt, = a a f2 (t) b h1 (t) h2 (t) g(t) [a,b] h1 (t) z(s1 ) ds1 [a,b] h2 (t) z(s2 ) ds2 dt = a a a b b = z(s1 ) K1 (s1, s2 )z(s2 ) ds2 ds1, a a b K1 (t, s) = g( )[a,b] h1 ( ) [a,b] h2 ( ) 2 a [a,h1 ( )] (t)[a,h2 ( )] (s) + [a,h2 ( )] (t)[a,h1 ( )] (s) d, получим b b I1 (z) = z(t) z(t) K(t, s)z(s) ds + f012 (t)z(t) + (t) dt, a a где K(t, s) = K0 (t, s) K1 (t, s), f012 (t) = f0 (t) + f1 (t) + f2 (t).

Функция Лагранжа задачи (5.7a)–(5.7b) имеет вид b L(z, ) = I1 (z) + 1 2 + z(t) dt.

a Условие стационарности функции Лагранжа по z в точке z ( решение за z дачи):

b b z Lz (, )[] 2(t) z K(t, s)(s) 2f (t) ds (t) dt = 0, z a a Таким образом, если z является экстремалью задачи (5.7a)–(5.7b), то оно долж но быть корнем уравнения b z (t) K(t, s)(s) = f (t).

z a Здесь f (t) = (f012 (t) + ).

5.2.1 Постоянное отклонение аргумента Пусть h1 (t) = t 1, h2 (t) = t + 2. Тогда b K1 (t, s) = g( )[a,b] ( 1 )[a,b] ( + 2 ) 2 a [a, 1 ] (t)[a, +2 ] (s) + [a, 1 ] (s)[a, +2 ] (t) d.

Заметим, что подынтегральная функция обращается в ноль вне диапазона [a + 1, b 2 ], поэтому b K1 (t, s) = g( ) [a, 1 ] (t)[a, +2 ] (s) + [a, 1 ] (s)[a, +2 ] (t) d.

2 a+ Первое слагаемое подынтегралной функции [a, 1 ] (t)[a, +2 ] (s) обращается в ноль при t + 1 или s 2. Следовательно, для t, s [a, b] получаем b2 b 1 K1 (t, s) = g( ) d + g( ) d 2 min(b2,max(t+1,s2 )) min(b2,max(s+1,t2 )) или, используя первообразные, K1 (t, s) = G(b 2 ) G(min(b 2, max(t + 1, s 2 )))+ + G(min(b 2, max(s + 1, t 2 ))).

Слагаемое f1 (t) в правой части имеет вид b g(s)h1 (s)[a,b] h2 (s) [a,h2 (s)] (t) ds = f1 (t) = a b = g(s)1 h1 (s) (,a) h1 (s) [a,b] h2 (s) [a,h2 (s)] (t) ds+ a b + g(s)1 [a,b] h1 (s) [a,b] h2 (s) [a,h2 (s)] (t) ds+ a b + g(s)2 h1 (s) (b,) h1 (s) [a,b] h2 (s) [a,h2 (s)] (t) ds = a b = g(s)1 (s 1 )(,a) (s 1 )[a,b] (s + 2 )[a,s+2 ] (t) ds+ a b + g(s)1 [a,b] (s 1 )[a,b] (s + 2 )[a,s+2 ] (t) ds+ a b + g(s)2 (s 1 )(b,) (s 1 )[a,b] (s + 2 )[a,s+2 ] (t) ds.

a Здесь b g(s)1 (s 1 )(,a) (s 1 )[a,b] (s + 2 )[a,s+2 ] (t) ds = a min(a+1,b2 ) = g(s)1 (s 1 ) ds, min(max(a,t2 ),a+1,b2 ) b b g(s)1 [a,b] (s 1 )[a,b] (s + 2 )[a,s+2 ] (t) ds = g(s)1 ds, a max(a+1,t2 ) b g(s)2 (s 1 )(b,) (s 1 )[a,b] (s + 2 )[a,s+2 ] (t) ds 0.

a def def Определив первообразные GF1 (t) = g(t)1 (t 1 ) dt и G(t) = = = g(t) dt, получим f1 (t) = GF1 min(a + 1, b 2 ) GF1 min(max(a, t 2 ), a + 1, b 2 ) + + 1 G b 2 1 G max(a + 1, t 2 ).

Аналогично, для f2 (t) имеем b g(s)h2 (s)[a,b] h1 (s) [a,h1 (s)] (t) ds = f2 (t) = a b = g(s)1 h2 (s) (,a) h2 (s) [a,b] h1 (s) [a,h1 (s)] (t) ds+ a b + g(s)1 [a,b] h2 (s) [a,b] h1 (s) [a,h1 (s)] (t) ds+ a b + g(s)2 h2 (s) (b,) h2 (s) [a,b] h1 (s) [a,h1 (s)] (t) ds = a b = g(s)1 (s + 2 )(,a) (s + 2 )[a,b] (s 1 )[a,s1 ] (t) ds+ a b + g(s)1 [a,b] (s + 2 )[a,b] (s 1 )[a,s1 ] (t) ds+ a b + g(s)2 (s + 2 )(b,) (s + 2 )[a,b] (s 1 )[a,s1 ] (t) ds.

a Где b g(s)1 (s + 2 )(,a) (s + 2 )[a,b] (s 1 )[a,s1 ] (t) ds 0, a b max(b2,t+1 ) g(s)1 [a,b] (s + 2 )[a,b] (s 1 )[a,s1 ] (t) ds = g(s)1 ds, a t+ b g(s)2 (s + 2 )(b,) (s + 2 )[a,b] (s 1 )[a,s1 ] (t) ds = a max(b,t+1 ) = g(s)2 (s + 2 ) ds.

max(b2,t+1 ) Обозначив через GF2 (t) первообразную функции g(t)2 (t + 2 ), получим f2 (t) = 1 G max(b 2, t + 1 ) 1 G(t + 1 )+ + GF2 max(b, t + 1 ) GF2 max(b 2, t + 1 ).

Таким образом значения функций K(t, s) и f012 (t) вычисляются без ис пользования численного интегрирования.

Снова возмём a = 0, b = 10, p(t) = sin(2t) exp(t/5), q(t) = 1 + cos(t), а также и 1 = 2, 2 = 3. 1 1, 2 1.

g(t) = t2 + Тогда K0 (t, s) = P (10) P (max(t, s)), K1 (t, s) = G(7) (G(a1 ) + G(a2 )), 5 t/ P (t) = e 10 cos(2t) + sin(2t), G(t) = 5 arctg(t), a1 = min(7, max(t + 2, s 3)), a2 = min(7, max(s + 2, t 3)) и f0 (t) = 2(P (10) P (t)) + Q(10) Q(t), f1 (t) = GF1 (2) GF1 (min(max(0, t 3), 2)) (G(7) + G(max(2, t 3))), f2 (t) = GF2 (max(10, t + 2)) GF2 (max(7, t + 2))+ + G(max(7, t + 2)) G(t + 2), Q(t) = t + sin(t), GF1 (t) = 5 arctg(t), GF2 (t) = 5 arctg(t).

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

Для аппроксимации введём на [a, b] сетку = {0, 1, 2,..., 10}, узлы ко 0. -0. -0. -0. -0. - -1. -1. -1. 0 2 10 0 2 4 6 8 Рис. 5.5: K(t, s) 0. 0.01 0. -0. -0. -0. 0 2 10 0 2 4 6 8 Рис. 5.6: K(t, s) K(t, s) торой расположены с постоянным шагом, равным единице. В качестве базиса аппроксимации будем использовать многочлены Лежандра до третьей степени включительно.

Аппроксимация ядра K(t, s) даёт кусочно–бикубический многочлен K(t, s), мало отличающийся от исходной функции. График поверхности функции по грешности K(t, s) = K(t, s)K(t, s) и проекция её линий уровня на плоскость (t, s) приведены на рис. 5.6.

Графики функции правой части интегрального уравнения f (t) и погреш ности аппроксимации f (t) f (t) показаны рис. 5.7.

0. 0. - 0. 0. - - -0. - -0. 0 2 4 6 8 10 0 2 4 6 8 Рис. 5.7: f (t) и f (t) f (t) После оценки коэффициентов Ki1 i2 j1 j2 и fij были вычислены оценки сверху для норм K K 0. L (L2 [a,b]) и f f 0.074413434498845386272.

L2 [a,b] Оценка нормы резольвентного оператора:

R [1.9771270453048617188;

1.9771270640409375208] L (L2 [a,b]) Так как оператор K обратим и норма K K меньше нижней L (L2 [a,b]) границы интервальной оценки [0.33589429624232166135;

0.33589429835621901951] 1/(1 + R L (L2 [a,b]) ), то оператор K также имеет обратный.

Оценка верхней границы спектра оператора K даёт max (K) 0.64087265934202963802, что меньше, чем 1 K K Следовательно, экстремаль x доставляет L (L2 [a,b]) функционалу минимум.

После решения СЛАУ с правой частью, зависящей от параметра, было восстановлено решение исходной задачи x(t). Значением параметра, приводя щим к выполнению второго краевого условия x(10) = 1, будет [2.2617755252867937088;

2.2617755243554711342].

На рис. 5.8 приведены графики функций z (t) = x (t) и x(t).

2 1.5 1 - 0. - - -0. - - - -1. - - 0 2 4 6 8 10 0 2 4 6 8 Рис. 5.8: z (t) и x(t) Норма производной приближённого решения x [4.6253698739152611097;

4.6253698747206657416], L2 [a,b] и погрешность производной x x 0.027093625371319605128.

L2 [a,b] Вывод программы:

Аппроксимация f... Ok Аппроксимация K... Ok Создание матрицы A... Ok Создание матрицы B... Ok Обращение матрицы A... Ok Решение СЛАУ... Ok Кусочно полиномиальное представление решения... Ok Вычисление множителей Лагранжа... Ok lambda = [-2.2617755252867937088;

-2.2617755243554711342] Собираем приближённое решение интегрального уравнения... Ok ||tilde z|| = [4.6253698739152611097;

4.6253698747206657416] Собираем приближённое решение вариационной задачи... Ok ||tilde x|| = [13.763986684880748612;

13.763986687045989044] Оценка норм приближений... Ok ||tilde F|| = [7.8466967066809925058;

7.8466967076305111917] ||tilde K|| = [3.0774280508222049413;

3.0774280514274656717] Генерируем данные для графиков... Ok Оцениваем погрешность аппроксимации f...Ok max DeltaF = 0. ||DeltaF|| = 0. Оцениваем погрешность аппроксимации K...Ok max DeltaK = 0. ||DeltaK|| = 0. Оцениваем квадрат нормы резольвентного оператора... Ok ||R|| = [1.9771270453048617188;

1.9771270640409375208] 1/(1+||R||) = [0.33589429624232166135;

0.33589429835621901951] ||z - tilde z|| = 0. Вычисление характеристического многочлена для E-A... Ok Все корни внутри круга радиуса 3. Оцениваем верхнюю границу спектра... Ok max sigma = 0. Успешное завершение Затрачено 47.07сек.

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

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

Предлагаемый метод не всегда позволяет получить решение задачи, даже если оно существует. Возможные причины:

• оператор Q в (2.4a) не может быть представлен в виде разности I K, где I тождественный оператор, K самосопряжённый вполне непре рывный оператор;

• оператор I K не является регулярным;

• спектр оператора I K не содержит нулевых собственных чисел, однако в ходе вычислительного эксперимента не удаётся построить такой при ближённый конечномерный оператор Kn, чтобы было выполнено усло вие (2.13);

• система линейных алгебраических уравнений (3.7) имеет настолько боль шую размерность, что с помощью имеющейся в наличии вычислитель ной техники не удаётся получить её решение, построить характеристи ческий многочлен матрицы коэффициентов, построить систему Штурма или оценить норму резольвентного оператора за приемлемое время.

Литература [1] Азбелев, Н. В. Функционально–дифференциальные уравнения и вариаци онные задачи / Н. В. Азбелев, С. Ю. Култышев, В. З. Цалюк. Москва– Ижевск: НИЦ Регулярная и хаотическая динамика, 2006. 122 с.

[2] Азбелев, Н. В. Элементы современной теории функционально– дифференциальных уравнений. Методы и приложения / Н. В. Азбелев, В. П. Максимов, Л. Ф. Рахматуллина. М.: Институт компьютерных исследований, 2002. 384 с.

[3] Алексеев, В. М. Оптимальное управление / В. М. Алексеев, В. М. Тихо миров, С. В. Фомин. М.: Наука. Гл. ред. физ.-мат. лит., 1979. 432 с.

[4] Алефельд, Г. Введение в интервальные вычисления / Г. Алефельд, Ю. Херцбергер. М.: Мир, 1987. 360 с.

[5] Алтунин, А. Е. Модели и алгоритмы принятия решений в нечетких усло виях / А. Е. Алтунин, М. В. Семухин. Тюмень: Изд-во Тюменского государственного ун-та, 2000. 352 с. www.plink.ru/tnm.

[6] Арнольд, В. И. Математические принципы классической механики / В. И. Арнольд. М.: Эдиториал УРСС, 2000. 408 с.

[7] Арнольд, В. И. Математические принципы классической и небесной ме ханики / В. И. Арнольд, В. В. Козлов, А. И. Нейштадт. 2-е, перераб. и доп. изд. М.: Эдиториал УРСС, 2002. 416 с.

[8] Афанасьев, В. Н. Математическая теория конструирования систем управ ления / В. Н. Афанасьев, В. Б. Колмановский, В. Р. Носов. 2-е, доп изд. М.: Высш. шк., 1998. 574 с.

[9] Бабенко, К. И. Основы численного анализа / К. И. Бабенко. Москва– Ижевск: НИЦ Регулярная и хаотическая динамика, 2002. 848 с.

[10] Гарантированная точность решения систем линейных уравнений в евкли довых пространствах / С. К. Годунов, А. Г. Антонов, О. П. Кирилюк, В. И. Костин. Новосибирск: Наука, Сиб. отделение, 1988. 465 с.

[11] Голуб, Д. Х. Матричные вычисления / Д. Х. Голуб, Ч. Ф. Ван Лоун.

М.: Мир, 1999. 548 с.

[12] Груздев, А. А. О редукции экстремальных задач к линейным уравнениям в гильбертовом пространстве / А. А. Груздев // Изв. ВУЗов. Матема тика. 1993. № 5 (372). С. 36–42.

[13] Гусаренко, С. А. Оптимальное управление: экстремальные и вариацион ные задачи / С. А. Гусаренко. Пермь: Перм. ун-т, Перм. техн. ун-т, 2001. 86 с.

[14] Данфорд, Н. Линейные операторы. Общая теория / Н. Данфорд, Д. Шварц. М.: Издательство иностранной литературы, 1962. 896 с.

[15] Каменский, Г. А. Экстремумы функционалов с отклоняющимися аргу ментами / Г. А. Каменский, А. Л. Скубачевский. М.: МАИ, 1979.

54 с.

[16] Канторович, Л. В. Функциональный анализ / Л. В. Канторович, Г. П. Акилов. М.: Наука. Гл. ред. физ.-мат. лит., 1977. 744 с.

[17] Кнут, Д. Э. Искусство программирования / Д. Э. Кнут. 3-е изд. М.:

Издательский дом Вильямс, 2000. Т. 1. Основные алгоритмы. 720 с.

[18] Кнут, Д. Э. Искусство программирования / Д. Э. Кнут. 3-е изд. М.:

Издательский дом Вильямс, 2000. Т. 2. Получисленные алгоритмы.

832 с.

[19] Колмогоров, А. Н. Элементы теории функций и функционального анали за / А. Н. Колмогоров, С. В. Фомин. 4-е, перераб изд. М.: Наука. Гл.

ред. физ.-мат. лит., 1976. 545 с.

[20] Ландау, Л. Теоретическая физика: учебное пособие в 10 т. / Л. Ландау, Е. М. Лифшиц. 4-е, испр. изд. М.: Наука. Гл. ред. физ.-мат. лит., 1988. Т. I. Механика. 216 с.

[21] Максимов, В. П. Вычислительный эксперимент при оптимизации процес сов с сегрегацией / В. П. Максимов, А. Н. Румянцев, В. А. Шишкин // Журнал физической химии. 1997. Т. 71, № 10. С. 1913–1916.

[22] Прасолов, В. В. Многочлены / В. В. Прасолов. 2-е, стереотипное изд.

М.: МЦНМО, 2001. 336 с.

[23] Пу, Т. Нелинейная экономическая динамика / Т. Пу. Ижевск: Изда тельский дом Удмуртский университет, 2000. 200 с.

[24] Рисс, Ф. Лекции по функциональному анализу / Ф. Рисс, Б. Сёкефальви Надь. М.: Мир, 1979. 569 с.

[25] Треногин, В. А. Функциональный анализ / В. А. Треногин. М.: ФИЗ МАТЛИТ, 2002. 488 с.

[26] Фаддеев, Д. Вычислительные методы линейной алгебры / Д. Фаддеев, В. Фаддеева. СПб.: Лань, 2002. 736 с.

[27] Шарый, С. П. Интервальные алгебраические задачи и их численное реше ние: Дисс... докт. физ.–мат. наук: 01.01.07. Новосибирск, 2000. 327 с.

[28] Шишкин, В. А. Конструктивный подход к исследованию вариационных задач для квадратичных функционалов и его компьютерная реализация / В. А. Шишкин // Материалы докладов. Международная конференция “Информационные технологии в инновационных проектах”. Ижевск:

1999. С. 155–157.

[29] Шишкин, В. А. Конструктивный подход к исследованию вариационных задач для квадратичных функционалов / В. А. Шишкин // Эконо мическая кибернетика: методы и средства эффективного управления.

Пермь, 2000. С. 90–94.

[30] Шишкин, В. А. Конструктивный подход к исследованию вариационных задач для квадратичных функционалов / В. А. Шишкин // Материалы Всероссийского семинара “Теория сеточных методов для нелинейных кра евых задач”. Казань: 2000. С. 135–137.

[31] Шишкин, В. А. Использование принципа расширения при решении задач оптимизации в условиях неопределённости / В. А. Шишкин // Экономи ческая кибернетика: математические и инструментальные методы анали за, прогнозирования и управления. Пермь, 2004. С. 153–156.

[32] Шишкин, В. А. Вариационные задачи для квадратичных функционалов.

доказательный вычислительный эксперимент / В. А. Шишкин // Совре менные проблемы прикладной математики и математического моделиро вания: Материалы конференции. Воронеж: Воронежская государствен ная академия, 2005. С. 246.

[33] Шишкин, В. А. Вариационные задачи для квадратичных функционалов.

доказательный вычислительный эксперимент / В. А. Шишкин // Вест ник Тамбовского университета. Серия Естественные и технические на уки. 2006. Т. 11, № 3. С. 268–269.

[34] Chukwu, E. N. Stability and time–optimal cointrol of hereditary systems / E. N. Chukwu. USA: Academic Press, Inc, 1992. 509 pp.

[35] Higham, N. J. Accuracy and stability of numerical algorithms / N. J. High am. USA, Philadelphia: SIAM, 1996. 688 pp.

[36] Kaucher, E. Computer arithmetic, scientic computation and mathematical modelling / E. Kaucher, S. M. Markov, G. Mayer // IMACS Annals on Com puting and Appl. Math. 1992. no. 12.

[37] Kearfott, R. B. INTLIB: A Portable FORTRAN 77 Interval Standard Function Library / R. B. Kearfott. USA, University of Southwestern Louisiana, 1995.

[38] Kearfott, R. B. INTERVAL ARITHMETIC: A Fortran 90 Module for an Inter val Data Type / R. B. Kearfott. USA, University of Southwestern Louisiana, 1996.

[39] Kulisch, U. Computer Arithmetic in Theory and Practice / U. Kulisch, W. L. Miranker. New York: Academic Press, 1981.

[40] Kushner, B. A. The constructive mathematics of A. A. Markov / B. A. Kushn er // The Mathematical Association of America. Monthly. Vol. 113, no. 6.

Pp. 559–566.

[41] Maksimov, V. P. On constructing solutions of functional dierential sys tems with a guaranteed precision / V. P. Maksimov, A. N. Rumyantsev, V. A. Shishkin // Functional Dierential Equations. 1995. Vol. 3, no.

1–2. Pp. 135–144.

[42] Moore, R. E. Interval arithmetic and automatic analysis in digital computing / R. E. Moore / Stanford Univercity. Stanford: 1962. 134 pp.

[43] Moore, R. E. The automatic analysis and control of error in digital computa tion based on the use of interval numbers / R. E. Moore // Error in Digital Computation. Vol. 1. New York: John Wiley & Sons, Inc., 1965. Pp. 61– 130.

[44] Moore, R. E. Interval Analysis / R. E. Moore, C. T. Yang. Sunnyvale, Cal ifornia: Lockheed Aircraft Corporation, Missiles and Space Division, 1959.

Vol. 1. 46 pp.

[45] Neumaier, A. Interval Methods for Systems of Equations / A. Neumaier.

Cambridge: Cambridge Univercity Press, 1990.

[46] Schulman, L. S. Some dierential–dierence equations containing both ad vance and retardation / L. S. Schulman // J. Math. Phys. 1974. Vol. 15, no. 3. Pp. 295–298.

[47] Shishkin, V. A. Computer-assisted study of variational problems with quadrat ic functionals / V. A. Shishkin // Book of Abstracts. VI International Congress of Mathematical Modeling. University of Nizhny Novgorod, 2004. P. 123.

[48] Wheeler, J. A. Classical electrodinamics in terms of direct interparticle ac tion / J. A. Wheeler, R. P. Feynman // Reviews of Modern Physics. 1945.

Vol. 21, no. 3. Pp. 425–433.

[49] Wheeler, J. A. Interaction with the absorber as the mechanism of radiation / J. A. Wheeler, R. P. Feynman // Reviews of Modern Physics. 1945.

Vol. 17, no. 2 and 3. Pp. 157–179.

Приложение A Программная реализация Программная реализация расчётной части предлагаемого конструктивного под хода выполнена на языке C++. Использовался набор компиляторов GNU (GNU GCC 4.3.2) операционной системы GNU/Linux Fedora 101.

При программировании использовались процедуры и функции из следую щих внешних библиотек программ:

1. GNU MP (арифметика произвольной точности) и MPFR (вещественные числа большой точности);

2. PTHREAD (потоки POSIX);

3. OpenMP (параллельные вычисления).

A.1 Makele GOAL=freddy OBJS=main.o integral.o legendre.o fourier.o interval.o problem.o fredholm.o LIBS=-lpthread -lgmpxx -lgmp -lgomp -lmpfr CPPFLAGS=-O2 -Wall -march=core2 -msse4.1 -fopenmp #CPPFLAGS=-ggdb %.o:%.cpp g++ -c $(CPPFLAGS) -o $@ $ $(GOAL): $(OBJS) g++ -o $(GOAL) $(OBJS) $(LIBS) clean:

rm -f $(GOAL) $(OBJS) *~ http://fedoraproject.org/ A.2 Базовые данные datautil.hpp #i f n d e f DATAUTIL #define DATAUTIL #include v a l a r r a y #i f d e f _WIN #include windows. h #e l s e #include p t h r e a d. h #include s y s / s y s i n f o. h #e ndif namespace MATH { #i f d e f _WIN typedef double DOUBLE;

i n l i n e s i z e _ t number_ of_ processors ( ) { SYSTEM_INFO s i ;

GetSy stemInfo (& s i ) ;

return s i. dwNumberOfProcessors ;

} #e l s e typedef long double DOUBLE;

i n l i n e s i z e _ t number_ of_ processors ( ) { return get_ nprocs ( ) ;

} #e ndif st ruct PARAMETERS { } ;

typedef DOUBLE ( FUNC) ( const s t d : : v a l a r r a y DOUBLE& x, PARAMETERS p ) ;

template c l a s s T T ABS ( const T& x ) { i f ( x 0) return x ;

else return x ;

} template c l a s s T T MIN ( const T& x, const T& y ) { if (x y) return x ;

else return y ;

} template c l a s s T T MAX ( const T& x, const T& y ) { if (x y) return x ;

else return y ;

} template c l a s s T void SWAP (T& x, T& y ) { T z = x;

x = y;

y = z;

} // линейноеотображение $ [ a, b ][\ a l p h a, \ b e t a ] $ template c l a s s T T p h i ( const T& t, const T& a, const T& b, const T& alpha, const T& b e t a ) { return a l p h a + ( ta ) ( betaa l p h a ) / ( ba ) ;

} };

#e ndif A.3 Математические процедуры A.3.1 Оценка значения интеграла integral.hpp #i f n d e f INTEGRAL #define INTEGRAL #include " d a t a u t i l. hpp" namespace MATH { // Численное интегрирование функции одной переменной // f подынтегральная функция // lb, ub границы области интегрирования // error оценка погрешности DOUBLE i n t e g r a l (DOUBLE ( f ) (DOUBLE), DOUBLE lb, DOUBLE ub, DOUBLE& e r r o r ) ;

// Численное интегрирование функции нескольких переменных // f подынтегральная функция // lb, ub границы области интегрирования // error оценка погрешности // p параметр, передаваемый в подынтегральную функцию // max_div максимально допустимое число итераций DOUBLE i n t e g r a l (FUNC f, const s t d : : v a l a r r a y DOUBLE& lb, const s t d : : v a l a r r a y DOUBLE& ub, DOUBLE& e r r o r, PARAMETERS p = NULL, s i z e _ t max_div = 1 0 0 0 ) ;

// правило с каким числом точек испльзовать #define USE_GK61 #define USE_GK21 };

#e ndif integral.cpp cmath #include v e c t o r #include l i s t #include a l g o r i t h m #include #include " i n t e g r a l. hpp" namespace MATH { namespace INTEGRAL { #i f USE_GK // Узлы и веса (30-61)-точечного метода Гаусса–Кронрода // gkf61.f (quadpack, www.netlib.org) const DOUBLE wg [ ] = { 0.007968192496166605615465883474674,0.018466468311090959142302131912047, 0.028784707883323369349719179611292,0.038799192569627049596801936446348, 0.048402672830594052902938140422808,0.057493156217619066481721689402056, 0.065974229882180495128128515115962,0.073755974737705206268243850022191, 0.080755895229420215354694938460530,0.086899787201082979802387530715126, 0.092122522237786128717632707087619,0.096368737174644259639468626351810, 0.099593420586795267062780282103569,0.101762389748405504596428952168554, 0. };

const DOUBLE xgk [ ] = { 0.999484410050490637571325895705811,0.996893484074649540271630050918695, 0.991630996870404594858628366109486,0.983668123279747209970032581605663, 0.973116322501126268374693868423707,0.960021864968307512216871025581798, 0.944374444748559979415831324037439,0.926200047429274325879324277080474, 0.905573307699907798546522558925958,0.882560535792052681543116462530226, 0.857205233546061098958658510658944,0.829565762382768397442898119732502, 0.799727835821839083013668942322683,0.767777432104826194917977340974503, 0.733790062453226804726171131369528,0.697850494793315796932292388026640, 0.660061064126626961370053668149271,0.620526182989242861140477556431189, 0.579345235826361691756024932172540,0.536624148142019899264169793311073, 0.492480467861778574993693061207709,0.447033769538089176780609900322854, 0.400401254830394392535476211542661,0.352704725530878113471037207089374, 0.304073202273625077372677107199257,0.254636926167889846439805129817805, 0.204525116682309891438957671002025,0.153869913608583546963794672743256, 0.102806937966737030147096751318001,0.051471842555317695833025213166723, 0. };

const DOUBLE wgk [ ] = { 0.001389013698677007624551591226760,0.003890461127099884051267201844516, 0.006630703915931292173319826369750,0.009273279659517763428441146892024, 0.011823015253496341742232898853251,0.014369729507045804812451432443580, 0.016920889189053272627572289420322,0.019414141193942381173408951050128, 0.021828035821609192297167485738339,0.024191162078080601365686370725232, 0.026509954882333101610601709335075,0.028754048765041292843978785354334, 0.030907257562387762472884252943092,0.032981447057483726031814191016854, 0.034979338028060024137499670731468,0.036882364651821229223911065617136, 0.038678945624727592950348651532281,0.040374538951535959111995279752468, 0.041969810215164246147147541285970,0.043452539701356069316831728117073, 0.044814800133162663192355551616723,0.046059238271006988116271735559374, 0.047185546569299153945261478181099,0.048185861757087129140779492298305, 0.049055434555029778887528165367238,0.049795683427074206357811569379942, 0.050405921402782346840893085653585,0.050881795898749606492297473049805, 0.051221547849258772170656282604944,0.051426128537459025933862879215781, 0. };

#e l s e // Узлы и веса (10-21)-точечного метода Гаусса–Кронрода // gkf21.f (quadpack, www.netlib.org) const DOUBLE wg [ ] = { 0.066671344308688137593568809893332,0.149451349150580593145776339657697, 0.219086362515982043995534934228163,0.269266719309996355091226921569469, 0. };

const DOUBLE xgk [ ] = { 0.995657163025808080735527280689003,0.973906528517171720077964012084452, 0.930157491355708226001207180059508,0.865063366688984510732096688423493, 0.780817726586416897063717578345042,0.679409568299024406234327365114874, 0.562757134668604683339000099272694,0.433395394129247190799265943165784, 0.294392862701460198131126603103866,0.148874338981631210884826001129720, 0. };

const DOUBLE wgk [ ] = { 0.011694638867371874278064396062192,0.032558162307964727478818972459390, 0.054755896574351996031381300244580,0.075039674810919952767043140916190, 0.093125454583697605535065465083366,0.109387158802297641899210590325805, 0.123491976262065851077958109831074,0.134709217311473325928054001771707, 0.142775938577060080797094273138717,0.147739104901338491374841515972068, 0. };

#e ndif const s i z e _ t knum = s i z e o f ( wgk ) / s i z e o f ( wgk [ 0 ] ) ;

// Версия метода Гаусса–Кронрода для функции одной переменной // (не адаптивная) DOUBLE gk (DOUBLE ( f ) (DOUBLE), DOUBLE lb, DOUBLE ub, DOUBLE& e r r o r ) { DOUBLE c e n t e r = ( ub+l b ) / DOUBLE( 2 ) ;

DOUBLE h l e n = ( ubl b ) / DOUBLE( 2 ) ;

DOUBLE r e s g = 0 ;

DOUBLE r e s k = wgk [ knum1] f ( c e n t e r ) ;

f o r ( s i z e _ t i = 0 ;

i ( knum 1);

i ++) { DOUBLE s h i f t = h l e n xgk [ i ] ;

DOUBLE fsum = ( f ( c e n t e r s h i f t )+ f ( c e n t e r+ s h i f t ) ) ;

r e s k += fsum wgk [ i ] ;

// если i нечётно, то этот узел также // применяется в формуле Гаусса i f ( ( i &1) != 0 ) r e s g += fsum wg [ i 1];

} DOUBLE r e s u l t = r e s k h l e n ;

// очень грубая оценка погрешности e r r o r = f a b s ( ( r e s k r e s g ) h l e n ) ;

return r e s u l t ;

} // Версия метода Гаусса–Кронрода для функции нескольких переменных // Глобально адаптивная. Рекурсивная.

DOUBLE gk (FUNC f, s t d : : v a l a r r a y DOUBLE & x, const s t d : : v a l a r r a y DOUBLE& c e n t e r, const s t d : : v a l a r r a y DOUBLE& hlen, DOUBLE& e r r o r, PARAMETERS p = NULL, size_t l e v e l = 0) { bool use_ func = ( l e v e l +1) == x. s i z e ( ) ;

DOUBLE r e s g = 0 ;

// если level + 1 равно длине x, то // следует вычислять значение подынтегральной функции // (достигнута максимальная глубина рекурсии), в противном случае // рекурсивно вызывается сама процедура интегрирования gk x [ le ve l ] = center [ l ev el ] ;

DOUBLE r e s k = wgk [ knum1] ( use_ func ? f ( x, p ) : gk ( f, x, c e n t e r, hlen, e r r o r, p, l e v e l + 1 ) ) ;

f o r ( s i z e _ t i = 0 ;

i ( knum 1);

i ++) { DOUBLE s h i f t = h l e n [ l e v e l ] xgk [ i ] ;

x [ l e v e l ] = c e n t e r [ l e v e l ] s h i f t ;

DOUBLE f v a l 1 = use_ func ? f ( x, p ) : gk ( f, x, c e n t e r, hlen, e r r o r, p, l e v e l + 1);

x [ l e v e l ] = c e n t e r [ l e v e l ]+ s h i f t ;

DOUBLE f v a l 2 = use_ func ? f ( x, p ) : gk ( f, x, c e n t e r, hlen, e r r o r, p, l e v e l + 1);

DOUBLE fsum = f v a l 1 + f v a l 2 ;

r e s k += fsum wgk [ i ] ;

// если i нечётно, то этот узел используется // и в формуле Гаусса i f ( ( i &1) != 0 ) r e s g += fsum wg [ i 1];

} DOUBLE r e s u l t = r e s k h l e n [ l e v e l ] ;

// очень грубая оценка погрешности // (только на наименьшем уровне вложенности) i f ( l e v e l == 0 ) { e r r o r = f a b s ( ( r e s k r e s g ) h l e n [ l e v e l ] ) ;

} return r e s u l t ;

} // Область, по которой вычисляется оценка интеграла st ruct REGION { s t d : : v a l a r r a y DOUBLE lb, ub ;

DOUBLE v a l u e, e r r o r ;

REGION ( ) : v a l u e ( 0 ), e r r o r ( 0 ) {} REGION ( s t d : : v a l a r r a y DOUBLE _lb, s t d : : v a l a r r a y DOUBLE _ub) : l b ( _lb ), ub (_ub ), v a l u e ( 0 ), e r r o r ( 0 ) {} REGION ( const REGION& d ) : l b ( d. l b ), ub ( d. ub ), v a l u e ( d. v a l u e ), e r r o r ( d. e r r o r ) {} REGION& operator= ( const REGION& d ) ;

// номер измерения (ребра области) с максимальной длиной s i z e _ t get_max_dim (DOUBLE& l e n ) const ;

};

REGION& REGION : : operator= ( const REGION& d ) { i f ( t h i s != &d ) { lb. r e s i z e (d. lb. s i z e ( ) ) ;

lb = d. lb ;

ub. r e s i z e ( d. ub. s i z e ( ) ) ;

ub = d. ub ;

} return t h i s ;

} s i z e _ t REGION : : get_max_dim (DOUBLE& l e n ) const { s t d : : v a l a r r a y DOUBLE L = ubl b ;

s i z e _ t dim = 0 ;

for ( size_t i = 1 ;

i L. s i z e ( ) ;

i ++) { i f (L [ i ] L [ dim ] ) dim = i ;

} l e n = L [ dim ] ;

return dim ;

} bool compare ( const REGION d1, const REGION d2 ) { return ( d1e r r o r ) ( d2e r r o r ) ;

} void d e l (REGION d ) { de le t e d ;

} // структура DATA нужна потому, что при создании потока // в функцию можно передавать указатель только на один // параметр st ruct DATA { FUNC f ;

PARAMETERS p ;

REGION r e g i o n ;

// список всех подобласте,й отсортированный // в порядке убывания погрешностей s t d : : l i s t REGION l s t ;

};

#i f d e f _WIN DWORD WINAPI e s t i m a t e (LPVOID p ) #e l s e void e s t i m a t e ( void p ) #e ndif { DATA data = static_cast & DATA(p ) ;

// определяем центр и полудлины рёбер области интегрирования s t d : : v a l a r r a y DOUBLE c e n t e r = ( data. r e g i o n ub+data. r e g i o n l b ) / DOUBLE( 2 ) ;

s t d : : v a l a r r a y DOUBLE h l e n = ( data. r e g i o n ubdata. r e g i o n l b ) / DOUBLE( 2 ) ;

// оцениваем значение интеграла // (по умолчанию значение уровня 0) s t d : : v a l a r r a y DOUBLE x ( c e n t e r. s i z e ( ) ) ;

data. r e g i o n v a l u e = gk ( data. f, x, c e n t e r, hlen, data. r e g i o n e r r o r, data. p ) ;

// вставляем данные в список // (для каждого потока список отдельный, поэтому // не требуется никаких блокировок) s t d : : l i s t REGION :: i t e r a t o r i t = data. l s t. b e g i n ( ) ;

f o r ( ;

i t != data. l s t. end ( ) ;

i t ++) { i f ( data. r e g i o n e r r o r ( i t ) e r r o r ) break ;

} data. l s t. i n s e r t ( i t, data. r e g i o n ) ;

#i f d e f _WIN return NULL;

#e l s e return 0 ;

#e ndif } };

// оценка значения интеграла для функции одной переменной DOUBLE i n t e g r a l (DOUBLE ( f ) (DOUBLE), DOUBLE lb, DOUBLE ub, DOUBLE& e r r o r ) { return INTEGRAL : : gk ( f, lb, ub, e r r o r ) ;

} // оценка значения интеграла для функции многих переменных DOUBLE i n t e g r a l (FUNC f, const s t d : : v a l a r r a y DOUBLE& lb, const s t d : : v a l a r r a y DOUBLE& ub, DOUBLE& e r r o r, PARAMETERS p, s i z e _ t max_div ) { using namespace INTEGRAL;

// количество процессов должно быть не меньше двух, // но лучше, если оно будет равно числу процессоров/ядер s i z e _ t n = number_ of_ processors ( ) ;

i f ( n 2) n = 2 ;

// для каждого потока создаётся свой набор данных s t d : : v e c t o r DATA i d ( n ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { i d [ i ]. r e g i o n = new INTEGRAL : : REGION( lb, ub ) ;

id [ i ]. l s t. c le a r ( ) ;

id [ i ]. f = f ;

id [ i ]. p = p ;

} DOUBLE v a l = 0, e r r = 0 ;

// пока не достигли максимально допустимого числа итераций // или требуемой точности, разбиваем область и // оцениваем значения интегралов на подобластях f o r ( s i z e _ t c o u n t e r = 0 ;

c o u n t e r max_div ;

c o u n t e r ++) { // определяем, по какому измерению будет // область будет разбиваться на подобласти...

DOUBLE l e n, s t e p ;

s i z e _ t max_dim = i d [ 0 ]. r e g i o n get_max_dim ( l e n ) ;

s t e p = l e n /n ;

//... и производим разбиение i d [ 0 ]. r e g i o n ub [ max_dim ] = i d [ 0 ]. r e g i o n l b [ max_dim]+ s t e p ;

f o r ( s i z e _ t i = 1 ;

i n ;

i ++) { i d [ i ]. r e g i o n l b [ max_dim ] = i d [ i 1 ]. r e g i o n ub [ max_dim ] ;

i d [ i ]. r e g i o n ub [ max_dim ] = i d [ i ]. r e g i o n l b [ max_dim]+ s t e p ;

} // так как область будет обрабатываться, уменьшаем // val и err на соответствующие значения v a l = i d [ 0 ]. r e g i o n v a l u e ;

e r r = i d [ 0 ]. r e g i o n e r r o r ;

// запускаем процедуры интегрирования в несколько потоков #i f d e f _WIN HANDLE t h r e a d = new HANDLE[ n ] ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { t h r e a d [ i ] = C reateT hread (NULL, 0, e s t i m a t e,& ( i d [ i ] ), 0, NULL ) ;

} W a i t F o r M u l t i p l e O b j e c t s ( n, thread,TRUE, INFINITE ) ;

de le t e [ ] t h r e a d ;

#e l s e pthread_t t h r e a d = new pthread_t [ n ] ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { p t h r e a d _ c r e a t e(& t h r e a d [ i ], 0, e s t i m a t e,& ( i d [ i ] ) ) ;

} f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { p t h r e a d _ j o i n ( t h r e a d [ i ],NULL ) ;

} de le t e [ ] t h r e a d ;

#e ndif // обновляем значения интеграла и погрешности f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { v a l += i d [ i ]. r e g i o n v a l u e ;

e r r += i d [ i ]. r e g i o n e r r o r ;

} // определяем, в какой подобласти наибольшая погрешность...

s i z e _ t max_idx = 0 ;

f o r ( s i z e _ t i = 1 ;

i n ;

i ++) { i f ( ( i d [ i ]. l s t. b e g i n ()) e r r o r ( i d [ max_idx ]. l s t. b e g i n ()) e r r o r ) max_idx = i ;

} //... перемещаем её в id[i].region...

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { i d [ i ]. r e g i o n = new REGION( i d [ max_idx ]. l s t. b e g i n ( ) ) ;

} //... и удаляем из списка подобластей de le t e i d [ max_idx ]. l s t. b e g i n ( ) ;

i d [ max_idx ]. l s t. pop_ front ( ) ;

// если достигли требуемой точности, заканчиваем i f ( e r r = e r r o r ) break ;

} // очищаем списки f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { de le t e i d [ i ]. r e g i o n ;

f o r ( s t d : : l i s t REGION :: i t e r a t o r i t = i d [ i ]. l s t. b e g i n ( ) ;

i t != i d [ i ]. l s t. end ( ) ;

i t ++) { de le t e i t ;

} } error = err ;

return v a l ;

} };

A.3.2 Многочлены polynom.hpp #i f n d e f POLYNOM #define POLYNOM v e c t o r #include v a l a r r a y #include i o s t r e a m #include omp. h #include #include " d a t a u t i l. hpp" namespace MATH { // Простейшее представление многочленов // c0 + c1 x + c2 x2 + · · · + cn xn template c l a s s T c l a s s POLYNOM { // коэффициенты s t d : : v e c t o r T c ;


// символ переменной, который используется // при выводе на экран char varname ;

public :

POLYNOM ( ) : c ( 1 ), varname ( ’ x ’ ) { c [ 0 ] = T ( 0 ) ;

} POLYNOM ( s i z e _ t deg, char vn = ’ x ’ ) : c ( deg +1), varname ( vn ) {} POLYNOM ( const POLYNOM & p ) : c ( p. c ), varname ( p. varname ) {} T void r e s i z e ( s i z e _ t deg ) { c. r e s i z e ( deg + 1);

} T operator [ ] ( s i z e _ t i d x ) const { i f ( i d x = c. s i z e ( ) ) return T ( 0 ) ;

return c [ i d x ] ;

} T& operator [ ] ( s i z e _ t i d x ) { i f ( i d x = c. s i z e ( ) ) throw "POLYNOM: выходзаграницымассива" ;

return c [ i d x ] ;

} // Вычисление значения многочлена. Метод Горнера T operator ( ) (T x ) const ;

// Степень многочлена s i z e _ t d e g r e e ( ) const { return c. s i z e () 1;

} // Унарный минус // c0 c1 x c2 x2 · · · cn xn POLYNOM T operator ( ) const ;

// Нормализация многочлена // (если степень многочлена больше 0, то коэффициент в мономе // максимальной степени должен быть ненулевой) void n o r m a l i z a t i o n ( ) ;

// Арифметические операции POLYNOM & operator += ( const POLYNOM & p ) ;

T T POLYNOM & operator = ( const POLYNOM & p ) ;

T T };

template c l a s s T T POLYNOM T : : operator ( ) (T x ) const { T s = c [ c. s i z e () 1];

f o r ( s i z e _ t i = d e g r e e ( ) ;

i 0 ;

i ) { s = x ;

s += c [ i 1 ] ;

} return s ;

} template c l a s s T POLYNOM T POLYNOM T : : operator ( ) const { POLYNOM T r e s ( d e g r e e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i = d e g r e e ( ) ;

i ++) r e s [ i ] = c [ i ] ;

return r e s ;

} template c l a s s T void POLYNOM T : : n o r m a l i z a t i o n ( ) { size_t i = degree ( ) ;

f o r ( ;

i 0 ;

i ) { i f ( c [ i ] ! =T ( 0 ) ) break ;

} i f ( i != d e g r e e ( ) ) c. r e s i z e ( i + 1);

} template c l a s s T POLYNOM & POLYNOM T T : : operator += ( const POLYNOM & p ) T { s i z e _ t s z = MAX( c. s i z e ( ), p. c. s i z e ( ) ) ;

size_t prev_size = c. s i z e ( ) ;

c. r e s i z e ( sz ) ;

f o r ( s i z e _ t i = p r e v _ s i z e ;

i p. c. s i z e ( ) ;

i ++) c [ i ] = T ( 0 ) ;

f o r ( s i z e _ t i = 0 ;

i p. c. s i z e ( ) ;

i ++) c [ i ] += p. c [ i ] ;

normalization ( ) ;

return t h i s ;

} template c l a s s T T : : operator = ( const POLYNOM & p ) POLYNOM & POLYNOM T T { s i z e _ t s z = MAX( c. s i z e ( ), p. c. s i z e ( ) ) ;

size_t prev_size = c. s i z e ( ) ;

c. r e s i z e ( sz ) ;

f o r ( s i z e _ t i = p r e v _ s i z e ;

i p. c. s i z e ( ) ;

i ++) c [ i ] = T ( 0 ) ;

f o r ( s i z e _ t i = 0 ;

i p. c. s i z e ( ) ;

i ++) c [ i ] = p. c [ i ] ;

normalization ( ) ;

return t h i s ;

} // n-я производная многочлена template c l a s s T POLYNOM T d i f f ( const POLYNOM & p, s i z e _ t n = 1 ) T { i f ( n==0) return p ;

s i z e _ t deg = p. d e g r e e ( ) ;

i f ( degn ) return POLYNOM T ( ) ;

deg = n ;

POLYNOM T dp ( deg ) ;

f o r ( s i z e _ t i = 0 ;

i = deg ;

i ++) { dp [ i ] = p [ i+n ] ;

f o r ( s i z e _ t j = i +1;

j = i+n ;

j ++) dp [ i ] = T( j ) ;

} return dp ;

} // Неопределённый интеграл от многочлена template c l a s s T POLYNOM T i n t e g r a l ( const POLYNOM & p )T { POLYNOM T i p ( p. d e g r e e ( ) + 1 ) ;

ip [ 0 ] = 0;

f o r ( s i z e _ t i = 1 ;

i = i p. d e g r e e ( ) ;

i ++) i p [ i ] = p [ i 1]/T( i ) ;

return i p ;

} };

// Сложение многочленов template c l a s s T MATH: : POLYNOMT operator+ ( const MATH: :POLYNOM & x, T const MATH: :POLYNOM & y ) T { using namespace MATH;

s i z e _ t deg = MAX( x. d e g r e e ( ), y. d e g r e e ( ) ) ;

POLYNOM T z ( deg ) ;

f o r ( s i z e _ t i = 0 ;

i = deg ;

i ++) { z [ i ] = 0;

i f ( i =x. d e g r e e ( ) ) z [ i ] += x [ i ] ;

i f ( i =y. d e g r e e ( ) ) z [ i ] += y [ i ] ;

} z. normalization ( ) ;

return z ;

} template c l a s s T MATH: : POLYNOM T operator+ ( const MATH: :POLYNOM & p, const T& x ) T { MATH: : POLYNOM T r e s = p ;

r e s [ 0 ] += x ;

return r e s ;

} template c l a s s T i n l i n e MATH: : POLYNOM T operator+ ( const T& x, const MATH: :POLYNOM & p ) T { return p+x ;

} // Вычитание многочленов template c l a s s T MATH: : POLYNOMT operator ( const MATH: :POLYNOM & x, T const MATH: :POLYNOM & y ) T { using namespace MATH;

s i z e _ t deg = MAX( x. d e g r e e ( ), y. d e g r e e ( ) ) ;

POLYNOM T z ( deg ) ;

f o r ( s i z e _ t i = 0 ;

i = deg ;

i ++) { z [ i ] = 0;

i f ( i =x. d e g r e e ( ) ) z [ i ] += x [ i ] ;

i f ( i =y. d e g r e e ( ) ) z [ i ] = y [ i ] ;

} z. normalization ( ) ;

return z ;

} template c l a s s T MATH: : POLYNOM T operator ( const MATH: :POLYNOM & p, const T& x ) T { MATH: : POLYNOM T r e s = p ;

r e s [ 0 ] = x ;

return r e s ;

} template c l a s s T MATH: : POLYNOM T operator ( const T& x, const MATH: :POLYNOM & p ) T { T r e s = p ;

MATH: : POLYNOM r e s [ 0 ] += x ;

return r e s ;

} // Перемножение многочленов template c l a s s T T operator ( const MATH: :POLYNOM & x, MATH: : POLYNOM T const MATH: :POLYNOM & y ) T { MATH: : POLYNOM T z ( x. d e g r e e ()+ y. d e g r e e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i = z. d e g r e e ( ) ;

i ++) z [ i ] = 0 ;

f o r ( s i z e _ t i = 0 ;

i = x. d e g r e e ( ) ;

i ++) f o r ( s i z e _ t j = 0 ;

j = y. d e g r e e ( ) ;

j ++) { z [ i+j ] += x [ i ] y [ j ] ;

} z. normalization ( ) ;

return z ;

} template c l a s s T T operator ( const MATH: :POLYNOM & p, const T& x ) MATH: : POLYNOM T { MATH: : POLYNOM T r e s = p ;

f o r ( s i z e _ t i = 0 ;

i = r e s. d e g r e e ( ) ;

i ++) r e s [ i ] = x ;

return r e s ;

} template c l a s s T i n l i n e T operator ( const T& x, const MATH: :POLYNOM & p ) MATH: : POLYNOM T { return px ;

} // Деление многочленов // x, y делимое и делитель // quo частное // rem остаток template c l a s s T void d i v ( const MATH: : POLYNOM & x, const MATH: : POLYNOM & y, T T MATH: : POLYNOM & quo, MATH: : POLYNOM & rem ) T T { i f ( y. d e g r e e ()==0 && y [0]==T( 0 ) ) throw " d i v : делениенаноль" ;

// если делимое меньше делителя, то частное - ноль // и остаток равен делимому i f ( x. d e g r e e () y. d e g r e e ( ) ) { quo. r e s i z e ( 0 ) ;

quo [ 0 ] = 0 ;

rem = x ;

return ;

} size_t n = x. degree ( ) ;

// временное хранилище (делимое, частное и остаток) s t d : : v a l a r r a y T q ( n + 1);

f o r ( s i z e _ t i = 0 ;

i = n ;

i ++) q [ i ] = x [ i ] ;

// деление столбиком f o r ( s i z e _ t i = 0 ;

i = x. d e g r e e () y. d e g r e e ( ) ;

i ++) { q [ n ] = q [ n ]/ y [ y. degree ( ) ] ;

pragma omp p a r a l l e l f o r # f o r ( s i z e _ t j = 1 ;

j = y. d e g r e e ( ) ;

j ++) { q [ nj ] = q [ nj ] y [ y. d e g r e e () j ] q [ n ] ;

} n;

} n = x. d e g r e e () y. d e g r e e ( ) ;

quo. r e s i z e ( n ) ;

f o r ( s i z e _ t i = 0 ;

i = n ;

i ++) quo [ i ] = q [ i+y. d e g r e e ( ) ] ;

quo. n o r m a l i z a t i o n ( ) ;

rem. r e s i z e ( y. d e g r e e ( ) 1 ) ;

for ( size_t i = 0 ;

i y. degree ( ) ;

i ++) rem [ i ] = q [ i ] ;

rem. n o r m a l i z a t i o n ( ) ;

} template c l a s s T s t d : : ostream& operator ( s t d : : ostream& os, const MATH: : POLYNOM & p ) T { bool f i r s t = true ;

i f ( p [ 0 ] ! = T( 0 ) ) { o s p [ 0 ] ;

f i r s t = false ;

} f o r ( s i z e _ t i = 1 ;

i = p. d e g r e e ( ) ;

i ++) { i f ( p [ i ]==T( 0 ) ) continue ;

if ( first ) { i f ( p [ i ] T( 0 ) ) o s "" ;

f i r s t = false ;

} else { i f ( p [ i ] T( 0 ) ) o s " " ;

else o s " + " ;

} i f (MATH: : ABS( p [ i ] ) ! = T( 1 ) ) o s MATH: : ABS( p [ i ] ) " x" ;

else o s "x" ;

i f ( i 1) o s "^" i ;

} return o s ;

} namespace MATH { // Вычисление последовательности Штурма s для многочлена p template c l a s s T void sturm ( const MATH: : POLYNOM & p, T std : : vectorMATH: : POLYNOMT & s ) { s. r e s i z e (p. degree ()+1);

s [0] = p;

s [ 1 ] = d i f f (p ) ;

POLYNOM T q, r ;

for ( size_t i = 2 ;

i s. s i z e ( ) ;

i ++) { d i v ( s [ i 2], s [ i 1], q, r ) ;

s [ i ] = r ;

i f ( r. d e g r e e ()==0 && r [0]==T( 0 ) ) break ;

} } // Число корней многочлена на отрезке [a, b] (без учёта кратности) // lval и rval значения элементов последовательности Штурма // на границах отрезка template c l a s s T s i z e _ t root_num ( const s t d : : v e c t o r POLYNOM T & s, T a, T b, s t d : : v a l a r r a y T l v a l = NULL, s t d : : v a l a r r a y T r v a l = NULL) { s t d : : v a l a r r a y int sgn ( s. s i z e ( ) ) ;

i f ( l v a l !=NULL) l v a l r e s i z e ( s. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i s. s i z e ( ) ;

i ++) { T val = s [ i ] ( a ) ;

sgn [ i ] = SIGN ( v a l ) ;

i f ( l v a l !=NULL) ( l v a l ) [ i ] = v a l ;

} s i z e _ t sum1 = 0 ;

f o r ( s i z e _ t i = 1 ;


i s. s i z e ( ) ;

i ++) { i f ( sgn [ i ]==0) sgn [ i ] = sgn [ i 1 ] ;

else sum1 += ( sgn [ i ]==sgn [ i 1]) ? 0 : 1 ;

} i f ( r v a l !=NULL) r v a l r e s i z e ( s. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i s. s i z e ( ) ;

i ++) { T val = s [ i ] ( b ) ;

sgn [ i ] = SIGN ( v a l ) ;

i f ( r v a l !=NULL) ( r v a l ) [ i ] = v a l ;

} s i z e _ t sum2 = 0 ;

f o r ( s i z e _ t i = 1 ;

i s. s i z e ( ) ;

i ++) { i f ( sgn [ i ]==0) sgn [ i ] = sgn [ i 1 ] ;

else sum2 += ( sgn [ i ]==sgn [ i 1]) ? 0 : 1 ;

} return sum1sum2 ;

} // Число корней на отрезке [a, b] // Если абсолютное значение многочлена не больше eps, // то оно заменяется нулём template c l a s s T s i z e _ t root_num ( const s t d : : v e c t o r POLYNOMT & s, T a, T b, T eps = 0) { s t d : : v a l a r r a y int sgn ( s. s i z e ( ) ) ;

for ( size_t i = 0 ;

i s. s i z e ( ) ;

i ++) { T val = s [ i ] ( a ) ;

i f (ABS( v a l)=e p s ) v a l = 0 ;

sgn [ i ] = SIGN ( v a l ) ;

} s i z e _ t sum1 = 0 ;

f o r ( s i z e _ t i = 1 ;

i s. s i z e ( ) ;

i ++) { i f ( sgn [ i ]==0) sgn [ i ] = sgn [ i 1 ] ;

else sum1 += ( sgn [ i ]==sgn [ i 1]) ? 0 : 1 ;

} for ( size_t i = 0 ;

i s. s i z e ( ) ;

i ++) { T val = s [ i ] ( b ) ;

i f (ABS( v a l)=e p s ) v a l = 0 ;

sgn [ i ] = SIGN ( v a l ) ;

} s i z e _ t sum2 = 0 ;

f o r ( s i z e _ t i = 1 ;

i s. s i z e ( ) ;

i ++) { i f ( sgn [ i ]==0) sgn [ i ] = sgn [ i 1 ] ;

else sum2 += ( sgn [ i ]==sgn [ i 1]) ? 0 : 1 ;

} return sum1sum2 ;

} };

#e ndif A.3.3 Ортогональные функции ortho.hpp /Ортогональныефункции / #i f n d e f ORTHO #define ORTHO #include gmpxx. h #include " d a t a u t i l. hpp" #include " polynom. hpp" namespace MATH { // абстрактный базовый класс для функций, // ортогональных на отрезке [l, u] template c l a s s T c l a s s ORTHO { protected :

T l,u;

public :

ORTHO (T _l, T _u) : l ( _l ), u (_u) {} v i r t u a l ~ORTHO ( ) {} // вес v i r t u a l T w e i g h t (T x ) = 0 ;

// квадрат нормы v i r t u a l T norm2 ( ) = 0 ;

// нижняя и верхняя границы области, на которой // функции являются ортогональными T l o ( ) const { return l ;

} T h i ( ) const { return u ;

} // значение функции в точке x v i r t u a l T operator ( ) (T x ) const = 0 ;

};

// ортогональные многочлены Лежандра (точное представление) // ортогональны на отрезке [1, 1] с единичным весом c l a s s LEGENDRE : public ORTHO mpq_class, public POLYNOMmpq_class { public :

LEGENDRE ( s i z e _ t deg ) ;

v i r t u a l mpq_class norm2 ( ) ;

v i r t u a l mpq_class w e i g h t ( mpq_class x ) { return mpq_class ( 1 ) ;

} v i r t u a l mpq_class operator ( ) ( mpq_class x ) const { return ( dynamic_castconst POLYNOMmpq_class( t h i s ) ) ( x ) ;

} };

// ортогональные многочлены Лежандра (DOUBLE–представление) c l a s s LEGENDRE_DBL : public ORTHO DOUBLE, public POLYNOMDOUBLE { public :

LEGENDRE_DBL ( s i z e _ t deg ) ;

v i r t u a l DOUBLE norm2 ( ) ;

v i r t u a l DOUBLE w e i g h t (DOUBLE x ) { return 1 ;

} v i r t u a l DOUBLE operator ( ) (DOUBLE x ) const { return ( dynamic_castconst POLYNOM DOUBLE(t h i s ) ) ( x ) ;

} };

};

#e ndif legendre.cpp #include cmath #include " o r t h o. hpp" namespace MATH { // pn 2 = 2/(2n + 1) mpq_class LEGENDRE : : norm2 ( ) { return mpq_class ( 2, 2 d e g r e e ( ) + 1 ) ;

} mpq_class (mpq_class ( 1), mpq_class ( 1 ) ) LEGENDRE : : LEGENDRE ( s i z e _ t deg ) : ORTHO { r e s i z e ( deg ) ;

const mpq_class ZERO = mpq_class ( 0 ) ;

switch ( deg ) { case 0 :

( this ) [ 0 ] = mpq_class ( 1 ) ;

break ;

case 1 :

( this ) [ 0 ] = ZERO;

( this ) [ 1 ] = mpq_class ( 1 ) ;

break ;

case 2 :

( this ) [ 1 ] = ZERO;

( t h i s ) [ 0 ] = mpq_class ( 1, 2 ) ;

( t h i s ) [ 2 ] = mpq_class ( 3, 2 ) ;

break ;

case 3 :

( t h i s ) [ 0 ] = ( t h i s ) [ 2 ] = ZERO;

( t h i s ) [ 1 ] = mpq_class ( 3, 2 ) ;

( t h i s ) [ 3 ] = mpq_class ( 5, 2 ) ;

break ;

case 4 :

( t h i s ) [ 1 ] = ( t h i s ) [ 3 ] = ZERO;

( this ) [ 0 ] = mpq_class ( 3, 8 ) ;

( this ) [ 2 ] = mpq_class ( 1 5, 4 ) ;

( this ) [ 4 ] = mpq_class ( 3 5, 8 ) ;

break ;

case 5 :

( this ) [ 0 ] = ( t h i s ) [ 2 ] = ( t h i s ) [ 4 ] = ZERO;

( t h i s ) [ 1 ] = mpq_class ( 1 5, 8 ) ;

( t h i s ) [ 3 ] = mpq_class ( 3 5, 4 ) ;

( t h i s ) [ 5 ] = mpq_class ( 6 3, 8 ) ;

break ;

de fault :

{ s t d : : v a l a r r a y mpq_class p1 ( deg +1), p2 ( deg +1), p ( deg + 1);

// P p1 [ 0 ] = p1 [ 2 ] = p1 [ 4 ] = ZERO;

p1 [ 1 ] = mpq_class ( 1 5, 8 ) ;

mpq_class ( 3 5, 4 ) ;

p1 [ 3 ] = p1 [ 5 ] = mpq_class ( 6 3, 8 ) ;

// P p2 [ 1 ] = p2 [ 3 ] = ZERO;

p2 [ 0 ] = mpq_class ( 3, 8 ) ;

mpq_class ( 1 5, 4 ) ;

p2 [ 2 ] = p2 [ 4 ] = mpq_class ( 3 5, 8 ) ;

// Pn+1 = 2n+1 xPn n+1 Pn n n+ f o r ( s i z e _ t n = 5 ;

n deg ;

n++) { mpq_class c1 = mpq_class ( 2 n+1,n + 1);

mpq_class c2 = mpq_class ( n, n + 1);

p [ 0 ] = c2 p2 [ 0 ] ;

f o r ( s i z e _ t j = 1 ;

j = n ;

j ++) { p [ j ] = c1 p1 [ j 1] c2 p2 [ j ] ;

} p [ n+1] = c1 p1 [ n ] ;

p2 = p1 ;

p1 = p ;

} f o r ( s i z e _ t i = 0 ;

i = d e g r e e ( ) ;

i ++) { ( this ) [ i ] = p [ i ] ;

} } break ;

} } //============================================= DOUBLE LEGENDRE_DBL: : norm2 ( ) { return 2. 0 / ( 2. 0 d e g r e e ( ) + 1. 0 ) ;

} DOUBLE ( 1. 0, 1. 0 ) LEGENDRE_DBL: : LEGENDRE_DBL ( s i z e _ t deg ) : ORTHO { r e s i z e ( deg ) ;

switch ( deg ) { case 0 :

( this ) [ 0 ] = 1.0;

break ;

case 1 :

( this ) [ 0 ] = 0.0;

( this ) [ 1 ] = 1.0;

break ;

case 2 :

( this ) [ 1 ] = 0.0;

( t h i s ) [ 0 ] = 0.5;

( this ) [ 2 ] = 1. 5 ;

break ;

case 3 :

( this ) [ 0 ] = ( this ) [ 2 ] = 0. 0 ;

( t h i s ) [ 1 ] = 1.5;

( this ) [ 3 ] = 2. 5 ;

break ;

case 4 :

( this ) [ 1 ] = ( this ) [ 3 ] = 0. 0 ;

( this ) [ 0 ] = 0.375;

( this ) [ 2 ] = 3.75;

( this ) [ 4 ] = 4.375;

break ;

case 5 :

( this ) [ 0 ] = ( this ) [ 2 ] = ( this ) [ 4 ] = 0. 0 ;

( this ) [ 1 ] = 1. 8 7 5 ;

( t h i s ) [ 3 ] = 8.75;

( this ) [ 5 ] = 7. 8 7 5 ;

break ;

de fault :

{ s t d : : v a l a r r a y DOUBLE p1 ( deg +1), p2 ( deg +1), p ( deg + 1);

// P p1 [ 0 ] = p1 [ 2 ] = p1 [ 4 ] = 0. 0 ;

p1 [ 1 ] = 1.875;

8.75;

p1 [ 3 ] = p1 [ 5 ] = 7.875;

// P p2 [ 1 ] = p2 [ 3 ] = 0. 0 ;

p2 [ 0 ] = 0.375;

3.75;

p2 [ 2 ] = p2 [ 4 ] = 4.375;

// Pn+1 = 2n+1 xPn n n+1 Pn n+ for ( size_t n = 5 ;

n deg ;

n++) { (2.0 n+1.0)/(n +1.0);

DOUBLE c1 = DOUBLE c2 = n/( n +1.0);

p [ 0 ] = c2 p2 [ 0 ] ;

f o r ( s i z e _ t j = 1 ;

j = n ;

j ++) { p [ j ] = c1 p1 [ j 1] c2 p2 [ j ] ;

} p [ n+1] = c1 p1 [ n ] ;

p2 = p1 ;

p1 = p ;

} f o r ( s i z e _ t i = 0 ;

i = d e g r e e ( ) ;

i ++) { ( this ) [ i ] = p [ i ] ;

} } break ;

} } };

A.3.4 Операции линейной алгебры linalg.hpp #i f n d e f LINALG #define LINALG v e c t o r #include v a l a r r a y #include i o s t r e a m #include omp. h #include #include " d a t a u t i l. hpp" #include " polynom. hpp" namespace MATH { // Простейшее представление матриц template c l a s s T c l a s s MATRIX { s t d : : v e c t o r s t d : : v a l a r r a y T m;

public :

MATRIX ( ) {} MATRIX ( s i z e _ t r, s i z e _ t c ) { r e s i z e ( r, c ) ;

} MATRIX ( const MATRIX& x ) { t h i s = x ;

} void r e s i z e ( s i z e _ t r, s i z e _ t c ) ;

MATRIX & operator= ( const MATRIX & x ) ;

T T s i z e _ t rows ( ) const { return m. s i z e ( ) ;

} s i z e _ t c o l s ( ) const { return m [ 0 ]. s i z e ( ) ;

} T operator ( ) ( s i z e _ t r, s i z e _ t c ) const { return m[ r ] [ c ] ;

} T& operator ( ) ( s i z e _ t r, s i z e _ t c ) { return m[ r ] [ c ] ;

} // матрица преобразуется в диагональную // с элементами d на главной диагонали void d i a g ( const T& d = T ( 1 ) ) ;

// след матрицы T t r a c e ( ) const ;

// характеристический многочлен матрицы POLYNOM T c h a r p o l y ( ) const ;

MATRIX & operator+= ( const MATRIX & x ) ;

T T MATRIX & operator= ( const MATRIX & x ) ;

T T MATRIX & operator= ( const MATRIX & x ) ;

T T };

template c l a s s T void MATRIX T : : r e s i z e ( s i z e _ t r, s i z e _ t c ) { m. r e s i z e ( r ) ;

f o r ( s i z e _ t i = 0 ;

i r ;

i ++) m[ i ]. r e s i z e ( c ) ;

} template c l a s s T MATRIX & MATRIX T T : : operator= ( const MATRIX & x ) T { i f ( t h i s != &x ) { r e s i z e ( x. rows ( ), x. c o l s ( ) ) ;

m = x.m;

} return t h i s ;

} template c l a s s T void MATRIX T : : d i a g ( const T& d ) { f o r ( s i z e _ t i = 0 ;

i rows ( ) ;

i ++) { m[ i ] = T ( 0 ) ;

i f ( i c o l s ( ) ) m[ i ] [ i ] = d ;

} } template c l a s s T T MATRIX T : : t r a c e ( ) const { T s = T( 0 ) ;

f o r ( s i z e _ t i = 0 ;

i MIN( rows ( ), c o l s ( ) ) ;

i ++) { s += m[ i ] [ i ] ;

} return s ;

} // Вычисление характеристического многочлена // методом Леверье–Фаддеева template c l a s s T POLYNOM T MATRIX T : : c h a r p o l y ( ) const { i f ( rows ( ) ! = c o l s ( ) ) throw " c h a r p o l y : матрицанеквадратная" ;

// степень характеристического многочлена s i z e _ t n = rows ( ) ;

POLYNOM T p ( n ) ;

p [ n ] = T( 1 ) ;

MATRIXT A=t his, B ;

f o r ( s i z e _ t i = 1 ;

i = n ;

i ++) { p [ ni ] = A. t r a c e ( ) /T( i ) ;

// std::cout p[n-i] std::endl;

i f ( i == n ) break ;

B = A;

f o r ( s i z e _ t j = 0 ;

j n ;

j ++) { B( j, j ) += p [ ni ] ;

} A = ( t h i s ) B ;

} return p ;

} template c l a s s T MATRIX & MATRIX T T : : operator+= ( const MATRIX & x ) T { i f ( rows ( ) ! = x. rows ( ) != c o l s ( ) ! = x. c o l s ( ) ) throw "MATRIX T+= : несовпадающиеразмерности" ;

f o r ( s i z e _ t i = 0 ;

i rows ( ) ;

i ++) { i f ( t h i s != &x ) m[ i ] += x.m[ i ] ;

else m[ i ] = T ( 2 ) ;

} return t h i s ;

} template c l a s s T MATRIX & MATRIX T T : : operator= ( const MATRIX & x ) T { i f ( ( rows ( ) ! = x. rows ( ) ) | | ( c o l s ( ) ! = x. c o l s ( ) ) ) throw "MATRIX T= : несовпадающиеразмерности" ;

f o r ( s i z e _ t i = 0 ;

i rows ( ) ;

i ++) { i f ( t h i s != &x ) m[ i ] = x.m[ i ] ;

else m[ i ] = T ( 0 ) ;

} return t h i s ;

} template c l a s s T MATRIX & MATRIX T T : : operator= ( const MATRIX & x ) T { i f ( c o l s ( ) != x. rows ( ) ) throw "MATRIX T= : несовпадающиеразмерности" ;

MATRIXT r e s ( rows ( ), x. c o l s ( ) ) ;

size_t j, k ;

pragma omp p a r a l l e l f o r private ( j, k ) # f o r ( s i z e _ t i = 0 ;

i rows ( ) ;

i ++) f o r ( j = 0 ;

j x. c o l s ( ) ;

j ++) { r e s.m[ i ] [ j ] = T ( 0 ) ;

f o r ( k = 0 ;

k c o l s ( ) ;

k++) { r e s.m[ i ] [ j ] += m[ i ] [ k ] x.m[ k ] [ j ] ;

} } f o r ( s i z e _ t i = 0 ;

i rows ( ) ;

i ++) { m[ i ]. r e s i z e ( x. c o l s ( ) ) ;

m[ i ] = r e s.m[ i ] ;

} return t h i s ;

} // Решение системы линейных алгебраических уравнений.

// Метод Гаусса–Жордана с частичным выбором главного элемента template c l a s s T void SLAE_gauss ( const MATH: : MATRIX & A, T const MATH: : MATRIX & B, T MATH: : MATRIX & x ) T { i f (A. rows ( ) ! =A. c o l s ( ) ) throw "SLAE_gauss : неквадратнаяматрица" ;

i f (A. rows ( ) ! =B. rows ( ) ) throw "SLAE_gauss : несовпадающеечислостроквAиB" ;

// временная матрица, чтобы не портить A MATRIX T AA(A ) ;

x = B;

size_t i, j, k ;

f o r ( i = 0 ;

i AA. rows ( ) ;

i ++) { // выбор главного элемента size_t idx = i ;

f o r ( j = i ;

j AA. rows ( ) ;

j ++) { i f (AA( idx, i )AA( j, i ) ) i d x = j ;

} // обмениваем строки, чтобы потом // не пришлось возиться с матрицей перестановок i f ( i d x != i ) { f o r ( j = i ;

j AA. c o l s ( ) ;

j ++) { SWAP(AA( i, j ),AA( idx, j ) ) ;

} f o r ( j = 0 ;

j x. c o l s ( ) ;

j ++) { SWAP( x ( i, j ), x ( idx, j ) ) ;

} } // нормализация строки f o r ( j = i +1;

j AA. c o l s ( ) ;

j ++) { AA( i, j ) /= AA( i, i ) ;

} f o r ( j = 0 ;

j x. c o l s ( ) ;

j ++) { x ( i, j ) /= AA( i, i ) ;

} AA( i, i ) = T ( 1 ) ;

// обнуление столбца pragma omp p a r a l l e l f o r private ( k ) # f o r ( j = 0 ;

j AA. rows ( ) ;

j ++) { i f ( j==i ) continue ;

f o r ( k = i +1;

k AA. c o l s ( ) ;

k++) { AA( j, k ) = AA( i, k ) AA( j, i ) ;

} f o r ( k = 0 ;

k x. c o l s ( ) ;

k++) { x ( j, k ) = x ( i, k ) AA( j, i ) ;

} AA( j, i ) = T ( 0 ) ;

} } } // Обращение матрицы методом Гаусса template c l a s s T MATRIX T i n v e r s e ( const MATRIX & A) T { i f (A. rows ( ) ! =A. c o l s ( ) ) throw " i n v e r s e : неквадратнаяматрица" ;

MATRIX T B(A. rows ( ),A. c o l s ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i A. rows ( ) ;

i ++) f o r ( s i z e _ t j = 0 ;

j A. c o l s ( ) ;

j ++) { B( i, j ) = ( i==j ) ? T( 1 ) : T ( 0 ) ;

} MATRIXT x ;

SLAE_gauss (A, B, x ) ;

return x ;

} };

template c l a s s T MATH: : MATRIX T operator+ ( const MATH: : MATRIX & x, T const MATH: : MATRIX & y ) T { i f ( x. rows ( ) ! = y. rows ( ) | | x. c o l s ( ) ! = y. c o l s ( ) ) throw "MATRIX( + ) : несовпадающиеразмеры" ;

MATH: : MATRIX T res (x );

for ( size_t i = 0 ;

i y. rows ( ) ;

i ++) for ( size_t j = 0 ;

j y. c o l s ( ) ;

j ++) r e s ( i, j ) += y( i, j );

return r e s ;

} template c l a s s T MATH: : MATRIX T operator ( const MATH: : MATRIX & x, T const MATH: : MATRIX & y ) T { i f ( x. rows ( ) ! = y. rows ( ) | | x. c o l s ( ) ! = y. c o l s ( ) ) throw "MATRIX( ): несовпадающиеразмеры" ;

MATH: : MATRIX T res (x );

for ( size_t i = 0 ;

i y. rows ( ) ;

i ++) for ( size_t j = 0 ;

j y. c o l s ( ) ;

j ++) r e s ( i, j ) = y( i, j );

return r e s ;

} template c l a s s T T operator ( const MATH: : MATRIX & x, MATH: : MATRIX T const MATH: : MATRIX & y ) T { i f ( x. c o l s ( ) ! = y. rows ( ) ) throw "MATRIX( ) : числостолбцоввx неравночислустрокв y" ;

MATH: : MATRIXT r e s ( x. rows ( ), y. c o l s ( ) ) ;

size_t i, j, k ;

pragma omp p a r a l l e l f o r private ( j, k ) # f o r ( i = 0 ;

i x. rows ( ) ;

i ++) f o r ( j = 0 ;

j y. c o l s ( ) ;

j ++) { res ( i, j ) = T( 0 ) ;

f o r ( k = 0 ;

k x. c o l s ( ) ;

k++) r e s ( i, j ) += x ( i, k ) y ( k, j ) ;

} return r e s ;

} #e ndif A.3.5 Аппроксимация approx.hpp #i f n d e f APPROXIMATION #define APPROXIMATION #include " d a t a u t i l. hpp" #include " o r t h o. hpp" namespace MATH { // Вычисление коэффициента ряда Фурье при разложении функции f // в области d [lbi, ubi ] по системе ортогональных 10 i= // функций ortho.

// Значение интеграла оценивается с точностью error DOUBLE f o u r i e r _ c o e f f (DOUBLE ( f ) ( const s t d : : v a l a r r a y DOUBLE&, PARAMETERS p ), const s t d : : v a l a r r a y DOUBLE& lb, const s t d : : v a l a r r a y DOUBLE& ub, const s t d : : v e c t o r ORTHO DOUBLE & ortho, DOUBLE e r r o r = DOUBLE( 1 e 5), MATH: : PARAMETERS param = NULL ) ;

};

#e ndif fourier.cpp #include " approx. hpp" #include " i n t e g r a l. hpp" namespace MATH { namespace APPROXIMATION { // f аппроксимируемая функция // c1, c2 коэффициенты линейных функций // для преобразования [, ] [a, b] // ortho - базисные функции st ruct APPROX_PARAMETERS : public PARAMETERS { DOUBLE ( f ) ( const s t d : : v a l a r r a y DOUBLE&, PARAMETERS ) ;

s t d : : v a l a r r a y DOUBLE c1, c2 ;

const s t d : : v e c t o r ORTHODOUBLE o r t h o ;

PARAMETERS p ;

};

// подынтегральная функция в числителе DOUBLE app_func ( const s t d : : v a l a r r a y DOUBLE& x, PARAMETERS p ) { APPROX_PARAMETERS& param = static_cast APPROX_PARAMETERS(p ) ;

// отображаем [, ] на [ai, bi ] s t d : : v a l a r r a y DOUBLE f x ( x. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i x. s i z e ( ) ;

i ++) { f x [ i ] = param. c1 [ i ] + param. c2 [ i ] x [ i ] ;

} // f (x) n i (j (xk )) k= DOUBLE v a l = param. f ( fx, param. p ) ;

f o r ( s i z e _ t i = 0 ;

i x. s i z e ( ) ;

i ++) { v a l = ( ( param. o r t h o ) [ i ] ) ( x [ i ] ) ( param. o r t h o ) [ i ] w e i g h t ( x [ i ] ) ;

} return v a l ;

} } DOUBLE f o u r i e r _ c o e f f (DOUBLE ( f ) ( const s t d : : v a l a r r a y DOUBLE&, PARAMETERS ), const s t d : : v a l a r r a y DOUBLE& lb, const s t d : : v a l a r r a y DOUBLE& ub, const s t d : : v e c t o r ORTHO DOUBLE & ortho, DOUBLE e r r o r, MATH: : PARAMETERS param ) { APPROXIMATION : : APPROX_PARAMETERS p ;

p. f = f ;

p. o r t h o = &o r t h o ;

p. p = param ;

DOUBLE a l p h a = o r t h o [0] l o ( ) ;

DOUBLE b e t a = o r t h o [0] h i ( ) ;

// коэффицинты c1 и c2 определяются из условий // c1 + c2 = a и c1 + c2 = b size_t n = lb. s i z e ( ) ;

p. c1. r e s i z e ( n ) ;

p. c2. r e s i z e ( n ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { p. c1 [ i ] = ( a l p h a ub [ i ] b e t a l b [ i ] ) / ( alphab e t a ) ;

p. c2 [ i ] = ( l b [ i ]ub [ i ] ) / ( alphab e t a ) ;

} // интегрирование ведётся по области [, ]n s t d : : v a l a r r a y DOUBLE a ( alpha, n ), b ( beta, n ) ;

DOUBLE v a l u e = i n t e g r a l (APPROXIMATION : : app_func, a, b, e r r o r,&p ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { v a l u e /= o r t h o [ i ]norm2 ( ) ;

} return v a l u e ;

} };

A.3.6 Интервальная арифметика interval.hpp #i f n d e f INTERVAL #define INTERVAL #include i o s t r e a m #include gmpxx. h #include " d a t a u t i l. hpp" namespace MATH { // Простейшая реализация интервальных чисел.

// При выполнении арифметических действий используется аппаратная // поддержка вычислений с направленным округлением (Intel, AMD) c l a s s INTERVAL { DOUBLE l, u ;

public :

INTERVAL ( ) : l ( 0 ), u ( 0 ) {} INTERVAL (DOUBLE x ) : l ( x ), u ( x ) {} INTERVAL ( long numer, long denom ) ;

INTERVAL ( mpq_class& x ) ;

INTERVAL& operator= ( const INTERVAL& x ) { i f ( t h i s != &x ) { l = x. l ;

u = x.u;

} return t h i s ;

} DOUBLE l o ( ) const { return l ;

} DOUBLE& l o ( ) { return l ;

} DOUBLE h i ( ) const { return u ;



Pages:     | 1 || 3 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.