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

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

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


Pages:   || 2 | 3 |
-- [ Страница 1 ] --

ПЕРМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

УДК 517.97+517.688

На правах рукописи

Шишкин Владимир Андреевич

ДОКАЗАТЕЛЬНЫЙ ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ

В ИССЛЕДОВАНИИ ВАРИАЦИОННЫХ ЗАДАЧ

ДЛЯ КВАДРАТИЧНЫХ ФУНКЦИОНАЛОВ

05.13.18 Математическое моделирование, численные методы

и комплексы программ

Диссертация на соискание учёной степени кандидата

физико–математических наук

Научный Доктор физико–математических наук, руководитель профессор Максимов В. П.

Пермь 2009 Оглавление Используемые обозначения 4 Введение 7 1 Предварительные сведения 15 1.1 Примеры вариационных задач для квадратичных функци оналов............................... 1.2 Минимизация квадратичного функционала......... 1.3 Доказательные вычисления................... 2 Абстрактная вариационная задача для квадратичного функционала 2.1 Постановка задачи........................ 2.2 Условия существования решения............... 2.3 Проекционный метод исследования.............. 3 Доказательный вычислительный эксперимент 3.1 Алгоритм доказательного вычислительного эксперимента. 3.2 Приближённое интегральное уравнение............ 3.3 Обобщение на случай функций многих переменных..... 4 Программная реализация 4.1 Структура программы..................... 4.2 Вспомогательные процедуры.................. 4.3 Постановка и решение задачи................. 5 Примеры 5.1 Простейшая задача Лагранжа................. 5.2 Задача с сосредоточенным отклонением аргумента..... Заключение A Программная реализация A.1 Makele.............................. A.2 Базовые данные......................... A.3 Математические процедуры.................. A.4 Постановка задачи........................ A.5 Доказательный вычислительный эксперимент........ B GNU Octave. Вспомогательные процедуры B.1 Решение скалярного уравнения................ B.2 Решение задачи Коши...................... B.3 Решение краевой двухточечной задачи............ B.4 Графики............................. Используемые обозначения множество натуральных чисел.

N множество целых чисел.

Z множество рациональных чисел.

Q евклидово пространство n-мерных вещественных векторов.

Rn множество интервалов, границы которых принадлежат множеству IX множество интервалов, границы которых могут быть точ X. IR но (без погрешности) представлены используемой вычислительной техникой.

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

f 1 (t) функция, обратная к функции f (t). Для однозначных функций f 1 f (t) = t.

· X. X пространство, сопря банахово пространство с нормой X значение функционала f X в жённое пространству X. f, x X точке x X.

абстрактное гильбертово пространство со скалярным произведени H ем ·, · H.

вещественное гильбертово пространство функций, суммируе L2 () мых на с квадратом. При [a, b] будем писать L2 [a, b].

x, y = x()y() d.

L2 () множество линейных ограниченных операторов, действую L (X, Y) щих из банахова пространства X в банахово пространство Y. При X = Y вместо L (X, X) будем писать L (X).

T L (Y, X ) оператор, сопряжённый оператору T L (X, Y).

ядро оператора T : X Y:

ker T ker T = {x X | T x = 0}.

производная Фреше отображения F L (X, Y) [3, Fx () L (X, Y) x стр. 138]:

F ( + h) F () = Fx ()[h] + o( h ).

x x x Если x скаляр, то может использоваться обозначение F ().

x вторая производная Фреше [3, стр. 154].

Fxx () L (X, L (X, Y)) x спектр оператора T.

(T ) интервальные числа (элементы множества IX).

x, y, z нижняя и верхняя границы интервала x:

x, x x = [x, x].

середина интервала x:

mid x mid x = (x + x)/2.

суперпозиция функций f : Y Z и g : X Y:

f g f g : X Z.

мера множества.

mes характеристическая функция множества X:

X 1, x X, X (x) = 0, x X.

n след n n-матрицы M : Sp M = Mii.

Sp M i= Введение Актуальность исследования. При изучении динамических систем в естественных, технических и экономических науках часто будущее состо яние исследуемой системы зависит не только от её настоящего, но также и от предыстории развития [34]. При построении достаточно точных и адекватных моделей сложных систем в большинстве случаев приходит ся учитывать запаздывание, возникающее вследствие конечных скоро стей при распространении информации и перемещении материальных объектов. В некоторых задачах требуется учитывать также и будущее состояние системы [46, 49]. Всё это приводит к тому, что при исследова нии и оптимизации деятельности таких систем приходится иметь дело с функционалами, содержащими не только локальные, но и с нелокальные операторы [15].

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

Классический подход к управлению линейными системами с квадра тичным функционалом описан, например, в [8]. В [1, 2] рассматривает ся решение задачи минимизации квадратичного функционала на осно ве идей и результатов теории абстрактного функционально–дифферен циального уравнения, развитых в работах Пермского семинара. Редук ция исходной задачи к задаче минимизации в гильбертовом пространстве описана в [12, 13].

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

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

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

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

Задачи:

1. разработать метод конструктивного исследования вариационных задач для квадратичных функционалов;

2. теоретически доказать применимость предлагаемого метода к ре шению вариационных задач с уравнением Эйлера в виде интеграль ного уравнения Фредгольма второго рода;

3. теоретически обосновать и разработать схему доказательного вы числительного эксперимента, ориентированного на получение ре зультатов с гарантированной оценкой точности решения;

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

5. оценить практическую применимость разработанного метода к ис следованию вариационных задач для квадратичных функционалов на основе его применения к тестовым модельным примерам;

под хода;

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

Объект исследования. В диссертационной работе исследуется зада ча минимизации вида N (0.1a) Ix = T1i x, T2i x + F0, x min, H X i= (0.1b) p(x) = 0, q(x) 0.

Решение ищется в банаховом пространстве X, изоморфном прямому про изведению вещественного сепарабельного гильбертова пространства H и пространства m-мерных вещественнх векторов Rm (X H Rm ). Пред полагается, что T1i, T2i, i = 1,..., N, линейные ограниченные опера торы, действующие из X в H;

F0 X ;

p и q аффинные вектор– функционалы.

В квадратичном функционале (0.1a) могут находиться, кроме самой функции x и её производных любого порядка, интегральные операторы, а также операторы с сосредоточенным или/и распределённым запазды ванием (интегральные операторы Стилтьеса).

Исследуемая задача (0.1a)–(0.1b) сводится к задаче минимизации в подходящем гильбертовом пространстве (0.2a) I1 z = Qz, z H f0, z H + 0 min, (0.2b) g(z) = 0, h(z) 0, где Q ограниченный самосопряжённый оператор. В работе рассмат риваются только те задачи, в которых оператор Q можно записать в виде разности I K, где I тождественный оператор, K вполне непрерывный и самосопряжённый оператор. В этом случае исследование необходимых и достаточных условий существования решения сводится к исследованию обратимости и положительной определённости оператора I K интегрального уравнения Фредгольма второго рода (0.3) z Kz = f.

Методологическая и теоретическая основа исследования. В хо де проверки необходимых и достаточных существования решения задачи (0.2a)–(0.2b), исследуемое уравнение Фредгольма z Kz = f заменяет ся близким уравнением z Kz = f с конечномерным оператором K.

Для этого строятся проекции K и f в конечномерные подпространства исходных пространств L (H) и H соответственно.

Если оператор I K обратим и положительно определён, то для про верки обратимости и положительной определённости исходного опера тора I K используются теорема о взаимной обратимости близких по норме операторов и теорема о спектре вполне непрерывного самосопря жённого оператора с вполне непрерывным самосопряжённым возмуще нием.

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

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

Структура диссертации.

Работа состоит из пяти глав.

В первой главе приводятся примеры вариационных задач для квад ратичных функционалов. Кратко описывается разработанный Пермским семинаром по функционально–дифференциальным уравнениям подход к исследованию задачи минимизации квадратичного функционала. Даётся описание методов конструктивной математики.

Вторая глава посвящена общему описанию предлагаемого метода исследования применительно к абстрактной вариационной задаче для квадратичного функционала. На основе изоморфизма X H Rm ис ходная задача редуцируется к задаче минимизации в гильбертовом про странстве.

Для проверки необходимых и достаточных условий существования решения предлагается заменить исходный объект (0.3) более простым уравнением Фредгольма второго рода с конечномерным ядром z Kn z = yn. Обосновывается выбор метода аппроксимации K и f на основе проекции в подпространства, натянутые на подмножества {1,..., n } заданной полной системы функций {i }, ортогональных с единичным весом.

Если уравнение z Kn z = yn имеет решение, то проверка разреши мости исходного уравнения осуществляется на основе теоремы об обрат ном операторе с помощью сравнения значений K Kn L (Hn ) и (I Kn )1 1(Hn ). Для проверки достаточного условия существования реше L ния исходной вариационной задачи оценивается нижняя граница спектра оператора I K, для чего используется теорема о спектре возмущённого ограниченного самосопряжённого компактного оператора.

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

В третьей главе описывается конструктивная реализация предла гаемого подхода для пространства H = L2 [a, b].

Даётся обоснование выбора в качестве базиса аппроксимации ортого нальных многочленов Лежандра или системы ортогональных функций Радемахера–Уолша.

Описаны два подхода к решению системы линейных алгебраических уравнений с учётом погрешностей аппроксимации.

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

В конце главы рассматривается случай H = L2 [a1, b1 ] · · · [ad, bd ].

Четвёртая глава посвящена описанию особенностей программной реализации.

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

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

Вторая задача также является задачей Лагранжа, но в ней уже ис пользуются элементы с отклоняющимся аргументом.

Приложение содержит тексты программ с краткими комментария ми.

Научная новизна исследования.

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

• На модельных примерах показано применение предлагаемого мето да для исследования вариационных задач для квадратичных функ ционалов, содержащих нелокальные оператры.

• Определены границы применимости предлагаемого подхода.

Апробация результатов исследования. Результаты диссертацион ного исследования были опубликованы в работах автора [28–33, 47] и в работах [21, 41], написанных в соавторстве. Из результатов совмест ных работ в диссертацию вошли только результаты, полученные автором лично.

Основные результаты диссертации докладывались и обсуждались на Пермском городском семинаре по функционально–дифференциальным уравнениям (1995–2006 гг., руководитель проф. Азбелев Н. В.), Ижев ском семинаре по дифференциальным уравнениям и задачам управле ния (1999 г., руководитель проф. Тонков Е. Л.), семинаре Лаборатории конструктивных методов исследования динамических моделей (в 2005 и 2008 гг., руководитель проф. Максимов В. П.), Пермском городском тео ретическом семинаре “Фундаментальные и прикладные модели управле ния экономикой” (2006 г., руководители проф. Аверин В. И. и проф. Пер ский Ю. К.), Междуродной конференции “Информационные технологии в инновационных проектах” (Ижевск, 1999 г.), III Всероссийском семина ре “Теория сеточных методов для нелинейных краевых задач” (Казань, 2000 г.), региональной конференции “Экономика и управление: актуаль ные проблемы и поиск путей решения” (Пермь, 2004 г.), IV Международ ном конгрессе по математическому моделированию (Нижний Новгород, 2004 г.), конференции “Современные проблемы прикладной математики и математического моделирования” (Воронеж, 2005 г.).

Глава Предварительные сведения Приведены примеры вариационных задач для квадратичных функционалов.

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

1.1 Примеры вариационных задач для квад ратичных функционалов 1.1.1 Классическая механика Рассмотрим замкнутую систему взаимодействующих материальных то чек [20, стр. 17]. Пусть в моменты t0 и t1 система занимает положения x(0) и x(1) соответственно. Согласно принципу стационарного действия (принцип Гамильтона) [7, стр. 30], траектория движения системы x(t) между этими положениями должна быть такова, чтобы интеграл t (1.1) L(x, x, t) dt t принимал экстремальное значение. Лагранжиан (интегрант) L в (1.1) имеет вид T U, где T кинетическая энергия системы, а U потен циальная энергия.

Если потенциальная энергия системы в определённый момент вре мени зависит только от расположения точек в этот момент времени, то можно сказать, что взаимодействия между точками распространяют ся мгновенно и mi x i (1.2) L= U (x1, x2,...).

i Здесь xi радиус–вектор i-й материальной точки, а mi её масса. Необ ходимость таких взаимодействий следует из основных предпосылок клас сической механики абсолютности времени и принципа относительно сти Галилея (во всех инерциальных системах отсчёта свойства простран ства и времени одинаковы и одинаковы все законы механики).

Зная лагранжиан (1.2) можно выписать уравнения Эйлера для функ ционала (1.1) уравнения движения:

d L L или mi xi = Fi, = dt xi xi где Fi = Ui /xi сила, действующая на i-ю точку.

Вид лагранжиана (1.2) показывает, что время не только однородно, но и изотропно, т.е. его свойства одинаковы в обоих направлениях, так как замена t на t оставляет неизменным лагранжиан, а значит и урав нения движения.

Если для описания движения используются не декартовы координа ты точек xi, а обобщённые координаты qi, то, сделав в (1.2) подстановку fi xi = fi (q1,..., qs ), xi = qk, qk k получим L= aik (q)qi qk U (q).

2 i,k Кинетическая энергия в обобщённых координатах является квадратич ной функцией от скорости, но может зависеть и от координат.

Пусть имеются две незамкнутые системы A и B, взаимодействующие друг с другом, причём система A + B, являющаяся объединением систем A и B, является замкнутой. Тогда L = TA (qA, qA ) + TB (qB, qB ) U (qA, qB ).

Если известно qB (t) движение системы B, то лагранжиан системы A имеет вид LA = TA (qA, qA ) U (qA, qB (t)) (лагранжиан определён с точностью до полной производной от любой функции координат и времени, поэтому член TB (qB (t), qB (t)) может быть опущен, так как он зависит только от времени и поэтому является пол ной производной от некоторой другой функции времени). Таким обра зом движение открытой системы во внешнем поле описывается обычным лагранжианом, но потенциальная энергия уже может зависеть от време ни.

Пример: колебания маятника Пусть требуется исследовать малые колебания двойного маятника, на ходящегося в однородном поле тяжести [20, стр. 21].

Точный лагранжиан (физический маятник) такой системы имеет вид m1 + m2 2 2 m2 2 L= l1 1 + l + m2 gl2 cos(2 )+ + m2 l1 l2 1 2 cos(1 2 ) + (m1 + m2 )gl1 cos(1 ) (g ускорение свободного падения). При малых x колебаниях значения косинусов близки к единице, поэтому исходный лагранжиан можно упростить l (математический маятник):

m 2 l2 m1 + m2 2 2 m2 2 2 L= l1 1 + l2 2 + m2 l1 l2 1 2 + 2 y m + (m1 + m2 )gl1 + m2 gl2.

В соответствии с принципом Гамильтона тра ектория движения системы должна быть экстремалью функционала t L dt.

t Пример: фигуры Лиссажу Пусть потенциальная энергия системы с двумя степенями свободы опре деляется следующим выражением [6, стр. 26]:

x2 + 2 x U= 1.

В этом случае траектория движения ( фигура Лиссажу ) будет экстре малью функционала t x2 + x2 x2 + 2 x 1 1 dt.

2 t Так как система состоит из двух подсистем, не связанных друг с дру гом, то движение по каждой координате можно рассматривать отдель но. Поэтому элементы траектории (x1, x2 ) должны быть экстремалями функционалов t1 t (x2 x2 ) dt (x2 2 x2 ) dt.

и 1 1 t0 t 1.1.2 Теория относительности Согласно современным опытным данным, максимальная скорость рас пространения взаимодействия конечна и ограничена скоростью света c.

Объединение принципа относительности (все законы природы одинако вы во всех инерциальных системах отсчёта) и конечной скорости рас пространения взаимодействий называется принципом относительности Эйнштейна.

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

Пример Рассмотрим упрощённую задачу взаимодействия двух объектов, распо ложенных на таком расстоянии l друг от друга, при котором становится заметно запаздывание передаваемого сигнала [48, p. 431].

Пусть имеется гармонический осциллятор t A B A, собственное движение которого описывает ся уравнением x(t) + 2 x(t) = (t), а на расстоянии l от осциллятора A находится другой осциллятор B. Колебание заряжен x ных частиц осциллятора A порождает электро магнитные волны, которые воздействуют на осциллятор B. Иницииро ванные колебания заряженных частиц осциллятора B в свою очередь порождают электромагнитные волны, которые действуют уже на A. Так как скорость распространения электромагнитных волн конечна и рав на скорости света c, между исходным движением частиц осциллятора A и индуцированным им взаимодействием от осциллятора B имеется вре менной лаг = l/c. Вследствие изотропности времени, замена t на t не должна менять уравнение движения, поэтому гармонический осцил лятор A будет описываться уравнением 1 x(t) + 2 x(t) = x(t ) + x(t + ) + (t).

2 Если внешнее воздействие на систему A + B равно нулю ( 0), то решение x должно быть экстремалью функционала действия [46] t x(t)2 2 x(t)2 + x(t /2)x(t + /2) dt.

t При этом предполагается, что заданы (известны) значения x(t) при t [t0 /2, t0 ] и t [t1, t1 + /2].

1.1.3 Экономика Рассмотрим простейшую односекторную модель монопольного производ ства1.

Обозначим в однопродуктовой модели количество произведённого то вара через x, а функцию стоимости производства через C(x). Пусть цена товара описывается непрерывно дифференцируемой функцией p. Пред положим, что рыночный спрос является функцией от цены и от скорости изменения цены D(p, p). Тогда объём прибыли на временном интервале Трофимов В. В. Геометрический анализ динамики больших экономических си стем. Статья в [23], стр. 174–198.

[t0, t1 ] определяется выражением t (1.3) D p, p p C D(p, p) dt, t причём предполагается, что рынок находится в состоянии равновесия, то есть объём производства равен объёму спроса.

Если D линейная функция от p и p, а C квадратичная функция от x, то (1.3) квадратичный функционал.

Пусть x = (x1,..., xn ), то есть монополист выпускает и продаёт n наименований товаров по ценам p = (p1,..., pn ). В этом случае задача максимизации прибыли на заданном промежутке времени будет задачей максимизации функционала n t Di (p, p)pi C D(p, p) dt max.

t0 i= Предположим, что в однопродуктовой модели изменение цены влияет на спрос с запаздыванием. Тогда t D p(t), p(t ) p C D(p(t), p(t )) dt.

t В многопродуктовой модели n t Di (p(t), p(t ))pi (t) C D(p(t), p(t )) dt max.

t0 i= Пусть Di зависит только от pi и pi (взаимонезависимые товары). Ко эффициенты в линейной функции Di могут быть константами:

Di (pi, pi ) = d0i + d1i pi (t) + d2i pi (t i ).

В более реалистичных моделях функция спроса часто должна иметь се зонную составляющую. Тогда Di (pi, pi, t) = d0i (t) + d1i (t)pi (t) + d2i (t)pi (t i ) и n t Di (pi (t), pi (t i ), t)pi (t) t0 i= C D1 (p1 (t), p1 (t i ), t),..., Dn (pn (t), pn (t n ), t) dt max.

1.1.4 Оптимальное управление Пусть поведение управляемой системы описывается системой дифферен циальных уравнений с отклоняющимся аргументом (1.4) x(t) = f t, x(t), x(h(t)) + u(t), где u управление. Требование минимизировать затраты (энергии, рез на управление приводит к задаче минимизации функционала t u(t)2 dt t или, более подробно, t1 (1.5) x(t) f t, x(t), x(h(t)) dt.

t Если f линейная функция (т.е. (1.4) система линейных дифферен циальных уравнений), то (1.5) квадратичный функционал.

1.2 Минимизация квадратичного функционала В классическом вариационном исчислении обычно исследуются задачи, в которых интегрантом является локальный оператор2. Пермский семи нар по функционально–дифференциальным уравнениям (ФДУ) разра ботал подход, позволяющий работать с функционалами более общего вида [1, 2, 12, 13], при этом решение ищется в пространстве, которое яв ляется в некотором смысле наиболее подходящим для рассматриваемой конкретной задачи.

Рассмотрим применение пермского подхода на следующем примере [2, стр. 182]: требуется минимизировать квадратичный функционал m b (1.6a) Ix = (T1i x)(s)(T2i x)(s) + (T0 x)(s) ds 2 a i= при наличии линейно независимой системы ограничений (1.6b) ix = i, i = 1,..., n.

Решение будем искать в банаховом пространстве X, изоморфном прямо му произведению L2 [a, b] Rn.

Полагаем, что операторы Tij, j = 1, 2, i = 1,..., m, и T0 являются эле ментами L (X, L2 [a, b]), а функционалы принадлежат L (X, R).

1,..., n Введём изоморфизм J = {, Y } : L2 [a, b] Rn X, построенный по данной системе функционалов 1,..., n на основе какой-нибудь одно значно разрешимой задачи x = z, ix = i, i = 1,..., n Оператор F : X Y, где X, Y функциональные пространства, называется локальным, если значения y(t) = (F x)(t) в любой окрестности точки t = t0 зависят только от значений x(t) в этой окрестности [2, стр. 9].

с линейным оператором : X L2 [a, b]. Известно, что такой оператор всегда существует (теорема 1.3.4 в [2, стр. 28]).

Подстановка x = z + Y в (1.6a) преобразует (1.6a)–(1.6b) в задачу о безусловной минимизации функционала def I1 (z) = I(z + Y ) = в гильбертовом пространстве L2 [a, b]. Обозначив Y = u, Qji = Tji, Q0 = T0, получим bm I1 (z) = (Q1i z)(s)(Q2i z)(s) ds+ 2 a i= bm + (Q1i z)(s)(T2i u)(s) + (Q2i z)(s)(T1i u)(s) ds+ 2 a i= m b + (T1i u)(s)(T2i u)(s) + (Q0 z)(s) + (T0 u)(s) ds.

2 a i= Это выражение можно переписать более компактно:

I1 (z) = Hz, z, z + q, L2 [a,b] L2 [a,b] где m (Q Q2i + Q Q1i ), (1.7) H= 1i 2i 2 i= m 1 (Q T2i + Q T1i )u Q (1), = 1i 2i 2 i= m b q= (T1i u)(s)(T2i u)(s) + (T0 u)(s) ds.

2 a i= Из определения H : L2 L2 видно, что это самосопряжённый опера тор.

Введём обозначение X = {x X | i x = i, i = 1,..., n}.

Точка x0 X (z0 L2 ) называется точкой локального минимума функционала I (I1 ), если существует такое 0, что I(x) I(x0 ) (I1 (z) I1 (z0 )) при всех x X (z L2 ), удовлетворяющих усло вию x x0 X ( z z0 L2 ). Если неравенство I(x) I(x0 ) (I1 (z) I1 (z0 )) выполняется при всех x X (z L2 ), то x0 (z0 ) называ ется точкой глобального минимума. Значение I(x0 ) (I1 (z0 )) называется локальным (соответственно глобальным) минимумом функционала [2, стр. 185].

Оператор H : L2 L2 называется положительно определённым, если 0 для всех z L2, и строго положительно определённым, если Hz, z он положительно определён и Hz, z = 0 только для z = 0.

Обозначим m 1 (Q T2,i + Q T1,i ), 0 = Q (1).

L= 1,i 2,i 2 i= Теорема 1.1 ( [2, стр. 185]). Всякий локальный минимум функционала (1.6a) является глобальным. Точка x0 X является точкой миниму ма функционала I тогда и только тогда, когда выполнены следующие условия 1. x0 является решением краевой задачи Lx = 0, ix = i, i = 1,..., n, 2. оператор H : L2 L2, определённый равенством (1.7), положи тельно определён.

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

Возможно, первые предложения по автоматическому анализу оши бок и верификации результатов компьютерных вычислений появились в конце 50-х начале 60-х годов прошлого века в работах Реймона Мура (см., например, [42–44]). Для этого было предложено использовать ин тервальное представление результатов вычислений. В частности, было показано, как можно применить такой подход при решении систем ли нейных алгебраических уравнений и обыкновенных дифференциальных уравнений.

В дальнейшем применение интервальных методов развивалось мно гими математиками. Можно отметить, например, труды Е. Каучера и др., У. Кулиша, В. Миранкера, Дж. Ньюмана [36,39,45]. На русском язы ке хороший обзор интервальных методов можно найти в монографии Г. Алефельда и Ю. Херцбергера [4].

В трудах С. П. Шарого рассматривается применение интервальных методов при решении систем уравнений (см., например, [27]).

Доказательные вычисления (вычисления с гарантированной точно стью) на ЭВМ рассматривались в трудах К. И. Бабенко (см., например, [9]) и С. К. Годунова (см., например, [10]).

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

Конструктивная математика коротко может быть охарактеризована следующими основными чертами [40, стр. 561]:

1. Предметом изучения служат конструктивные процессы и конструк тивные объекты, возникающие как результаты этих процессов.

2. Используется специальная конструктивная логика, которая учиты вает специфику конструктивных процессов и объектов.

3. Рассмотрение конструктивных процессов и объектов проводится в рамках абстракции потенциальной осуществимости с полным ис ключением идеи актуальной бесконечности.

4. Интуитивное понятие эффективности связывается с точным поня тием алгоритма.

В конструктивной математике начальными, базисными понятиями явля ются понятия конструктивного процесса и конструктивного объекта.

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

Абстракция потенциальной осуществимости даёт возможность не обращать внимания на реальные ограничения на ресурсы при рассуж дениях о выполнении операций. Если сегодня для решения задачи не хватает материалов, времени или пространства, то, возможно, эту зада чу можно будет решить в будущем. Также можно рассматривать такие объекты, которые никогда не смогут быть построены реально (напри мер, очень длинные формулы), но при построении которых всё же не используется понятие актуальной бесконечности.

Считается, что конструктивный объект с заданным свойством суще ствует, если можно указать потенциально выполнимый конструктивный процесс, результатом выполнения которого является построение этого объекта. Для описания конструктивных процессов используется понятие алгоритма.

Согласно [17, стр. 31–32], современное значение слова алгоритм во многом аналогично таким понятиям, как рецепт, метод, способ, процеду ра, программа, но при этом оно имеет дополнительный смысловой отте нок. Нельзя сказать, что алгоритм это просто набор конечного числа правил, задающих последовательность выполнения операций для реше ния задач определённого рода. Помимо этого, он имеет пять важных особенностей:

1. Конечность. Алгоритм всегда должен заканчиваться после выпол нения конечного числа шагов.

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

2. Определённость. Каждый шаг алгоритма должен быть точно опре делён. Действия, которые нужно выполнить, должны быть строго и недвусмысленно определены для каждого возможного случая.

3. Ввод. Алгоритм имеет некоторое (возможно, равное нулю) число входных данных, т.е. величин, которые задаются до начала его работы или определяются динамически во время его работы. Эти входные данные берутся из определённого набора объектов.

4. Вывод. У алгоритма есть одно или несколько выходных данных, т.е.

величин, имеющих вполне определённую связь с входными данны ми.

5. Эффективность. Обычно алгоритм считается эффективным, если все его операторы достаточно просты для того, чтобы их можно было точно выполнить в течение конечного промежутка времени с помощью карандаша и бумаги.

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

1.3.2 Машинно–вычислимые функции Рассматриваемые в данной работе конструктивные процессы реализуют ся в виде компьютерных программ.

Современная вычислительная техника обычно работает с числами, имеющими конечную точность. При вычислениях с плавающей точкой используется числовая система F R подмножество вещественных чисел, которые могут быть представлены в форме y = ±m et [35, стр. 40]. Система F характеризуется четыремя целыми параметрами:

1. основанием ;

2. точностью t;

3. диапазоном, в котором может изменяться значение порядка: emin e emax.

Мантисса m это целое число, удовлетворяющее неравенству 0 m t 1. Для единственности представления каждого y F предполагает ся, что m t1 для y = 0, то есть система является нормализованной.

Диапазон ненулевых чисел с плавающей точкой в F определяется выра жением emin 1 emax (1 t ).

|y| Большинство применяемых в настоящее время стандартных подпро грамм работает почти исключительно с нормализованными числами:

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

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

Для улучшения ситуации, был создан стандарт IEEE 754–1985, опре деляющий двоичные арифметические системы с плавающей точкой. При разработке данного стандарта учитывались принципы робастности3, эф фективности и переносимости программного обеспечения, возможность обработки арифметических исключений, обеспечение вычисления транс цендентных функций и арифметики высокой точности [35, стр. 45]. В дальнейшем этот стандарт был расширен возможностью работы с числа ми с плавающей точной, которые не зависели бы от основания (стандарт IEEE 854–1987).

В данном случае под робастностью понимается устойчивость результатов вычис лений при малых изменениях аргументов.

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

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

Одним из возможных решений данной проблемы является использо вание интервальных оценок [43, 44], когда значение истинного результа та неизвестно, но имеются границы интервала, в котором этот результат гарантированно находится. При этом можно использовать как стандарт ную интервальную арифметику [4], так и полную интервальную ариф метику Каучера [27, стр. 49], в которой допускаются интервалы вида [a, b], где b a.

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

1. Числа с плавающей точкой. Такой подход реализован, например, в пакете INTLIB [37, 38] на языке Fortran. При выполнении ариф метических операций используется аппаратно или программно ре ализованное направленное округление.

Поле над множеством таких чисел будем обозначать IF.

2. Рациональные числа, числители и знаменатели которых являются целыми числами неограниченной (а точнее, ограниченной только размером машинной памяти) числами. При реализации, например, в пакете символьных вычислений Maple фирмы Waterloo Maple Inc., могут быть использованы тип данных INTERVAL и функция evalr.

Поле над множеством таких чисел будем обозначать IQ.

Если несущественно, какое поле используется, будем использовать обо значение IR.

Другим вариантом вычислений, которые гарантируют верность по лученного результата, является работа с рациональными числами. Если x Q представить в виде отношения x1 /x2 целого числа x1 и нату рального числа x2, причём эти числа могут иметь произвольно высокую точность4 [18, стр. 304], то все арифметические операции можно будет совершать точно, без ошибок округления [18, стр. 373]. Обозначим через Q множество рациональных чисел, которые могут быть точно представ лены на используемой вычислимой платформе5.

Арифметика рациональных чисел реализована в математических па кетах компьютерной алгебры, таких как Maple, Mathematica, Maxima.

Заметим, что последний пакет (Maxima6 ) распространяется по лицензии GNU GPL7 и может свободно использоваться для любых целей.

При написании программ удобно использовать библиотеку GNU MP8, позволяющую работать с целыми и рациональными числами произволь ной точности, а также с числами с плавающей точкой, имеющими ман тиссу произовольной точности. Этот пакет распространяется по лицен зии GNU LGPL9.

Будем называть машинно–вычислимой или, более кратко, вычисли мой, функцию f, которая отображает IR в IR или Q в Q. Иногда тре буется, чтобы производные и интегралы от вычислимой функции также были вычислимыми функциями.

Точнее, точность (длина) таких чисел будет ограничена только объёмом вирту альной памяти.

Заметим, что мощность множества Q зависит от компьютера, на котором выпол няются вычисления (например, от того, какой объём оперативной памяти установлен на нём).

http://maxima.sourceforge.net/ru/ http://www.gnu.org/licenses/licenses.html#GPL http://gmplib.org/ http://www.gnu.org/copyleft/lesser.html Вычислимыми функциями являются многочлены с коэффициентами из IR или из Q.

Глава Абстрактная вариационная задача для квадратичного функционала Описывается постановка абстрактной вариационной задачи для квадратично го функционала. Приводятся необходимые и достаточные условия существо вания решения.

2.1 Постановка задачи 2.1.1 Задача минимизации в банаховом пространстве Будем называть абстрактной вариационной задачей для квадратичного функционала задачу поиска минимума функционала N (2.1a) I(x) = T1i x, T2i x + F0, x H X i= при ограничениях (2.1b) pi (x) = 0, i = 1,..., n1, (2.1c) qi (x) 0, i = 1,..., n2.

Здесь H вещественное сепарабельное гильбертово пространство;

X заданное вещественное банахово пространство (B-пространство). Линей ные ограниченные операторы T1i, T2i, i = 1,..., N, действуют из X в H.

значение линейного ограниченного функционала F0 X на F0, x X x;

pi, i = 1,..., n1, и qi, i = 1,..., n2 определённые на X аффинные функционалы:

i X, pi (x) = i, x i, i R, X i X, qi (x) = i, x i, i R.

X Предполагаем, что система ограничений совместна и невырождена.

Будем рассматривать задачи, где в качестве операторов Tki, k = 1, 2, i = 1,..., N, используются тождественные операторы, дифференциаль ные операторы любого порядка, интегральные операторы, а также опе раторы сосредоточенного и распределённого отклонения аргумента.

Решение задачи (2.1a)–(2.1c) существенно зависит от выбора про странства X. Подход к исследованию вариационных задач, разработан ный Пермским семинаром по функционально–дифференциальным урав нениям, основан на том, что для решение конкретной задачи ищется в пространстве, наиболее подходящем для данного случая.

Из бесконечномерных B-пространств наиболее близко к конечномер ным евклидовым пространствам, по крайней мере с точки зрения эле ментарных геометрических свойств, стоят гильбертовы пространства [14, стр. 269]. В дальнейшем будем предполагать, что пространство X, в ко тором ищется решение, изоморфно прямому произведению H Rm, где n1. В следующем разделе описывается метод редукции исходной m задачи к задаче минимизации в гильбертовом пространстве.

2.1.2 Редукция к задаче в гильбертовом простран стве Сведём задачу (2.1a)–(2.1c) к задаче минимизации в гильбертовом про странстве, воспользовавшись методом редукции, который был описан в разделе 1.2 (см. также [1, 2, 12, 13]).

Пусть некоторый линейный ограниченный оператор, отображает X на всё пространство H. Пусть при этом размерность ядра равна m, причём m n1. Выберем из (2.1b) первые m ограничений (при необ ходимости перенумеровав исходную систему ограничений), и составим оператор [, ], где = {1,..., m }. Пара [, ] отображает X на прямое произведение H Rm и ker[, ] = {0}, поэтому существует обратный оператор [, ]1 = {, Y } L (H Rm, X).

Таким образом на основе некоторой однозначно разрешимой краевой задачи (2.2) x = z, i, x = i, i = 1,..., m, X построен изоморфизм X H Rm. Подстановка в (2.1a)–(2.1c) решения задачи (2.2) (2.3) x = z + Y, = col{1,..., m }.

приводит к задаче минимизации в гильбертовом пространстве H:

(2.4a) I1 (z) = Qz, z H + f0, z H + 0, (2.4b) gi (z) = 0, i = 1,..., n1 m, (2.4c) hi (z) 0, i = 1,..., n2.

После подстановки (2.3) в (2.1a) i-й элемент суммы, с учётом веществен ности пространства H, будет иметь вид T1i x, T2i x = T1i (z + Y ), T2i (z + Y ) = H H = T1i z, T2i z + T1i z, T2i Y H + T1i Y, T2i z H + H 1 + T1i Y, T2i Y H = (T1i T2i + T2i T1i )z, z + 2 H + (T1i T2i + T2i T1i )Y, z H + T1i Y, T2i Y H.

Таким образом, в (2.4a) N (T1i T2i + T2i T1i ), Q= i= N f0 = F0 + (T1i + T2i )Y, i= N 0 = T1i Y, T2i Y + F0, Y X.

H i= Заметим, что оператор Q является самосопряжённым по построению.

После удаления первых m ограничений из (2.1b), получаем систему (2.5) pm+1 (x) = 0,..., pn1 (x) = 0.

Подстановка (2.3) даёт gi (z) pm+i (z + Y ) = = m+i, z + m+i, Y m+i = i, z ai, X X H i = 1,..., n1 m где i = m+i, ai = m+i m+i, Y X.

Аналогично, для ограничений–неравенств получаем hi (z) = i, z bi, i = 1,..., n2, H i = i, bi = i i, Y X.

2.2 Условия существования решения Одним из методов решения гладких оптимизационных задач с ограниче ниями в виде равенств и неравенств является метод множителей Лагран жа [3, стр. 252].

Функция Лагранжа задачи (2.4a)–(2.4c) имеет вид (2.6) L(z, 1, 2 ) = I1 (z) + 1, g(z) + 2, h(z) Rn2, Rn1 m где 1 (Rn1 m ) и 2 (Rn2 ).

Более подробно, n1 m L(z, 1, 2 ) = Qz, z + f0, z + 0 + 1i i, z ai + H H H 2 i= n + (), (2.7) + 2i i, z H bi = Qz, z f (), z H H i= где n1 m n f () = f0 1i i 2i i, i=1 i= n1 m n () = 0 1i ai 2i bi.

i=1 i= 2.2.1 Необходимые условия Известно (правило множителей Лагранжа для гладких задач с равен ствами и неравенствами, [3, стр. 252–253]), что если z доставляет ло кальный минимум в задаче (2.4a)–(2.4c), то найдутся такие не равные одновременно нулю множители Лагранжа 1 и 2, что будут выполне ны:

1. условие стационарности функции Лагранжа по z z Lz (, 1, 2 )[h] Q f (), z = 0, H;

H 2. условие согласования знаков: 2i 0, i = 1,..., n2 ;

3. условие дополняющей нежёсткости 2i hi () = 0, z i = 1,..., n2.

Из условия стационарности функции Лагранжа по z следует, что ре шение z задачи (2.4a)–(2.4c) должно быть корнем параметрического опе раторного уравнения (2.8) Qz = f ().

Запишем оператор Q в виде разности двух самосопряжённых опера торов A1 A2. Пусть оператор A1 имеет ограниченный обратный, а A является компактным оператором. В этом случае (2.8) можно переписать в виде (2.9) z Kz = y(), где K = A1 A2, y() = A1 f ().

1 Так как множество компактных операторов образует в кольце линейных ограниченных операторов двусторонний идеал, оператор K также будет компактным оператором [16, теорема 2, стр. 325]. Следовательно I K фредгольмов как сумма обратимого и компактного операторов [25, стр. 219–220].

Пусть в (2.1a) существует такой индекс i, что T1i = T2i. Если при ре дукции (2.1a)–(2.1c) к задаче минимизации в гильбертовом пространстве за оператор : X H принять T1i, то оператор Q сразу будет иметь вид I K, где (T1i T2i + T2i T1i ).

K= i=1,...,N i=i В дальнейшем будем рассматривать только те задачи (2.1a)–(2.1c), которые допускают представление соответствующего уравнения (2.8) в виде (2.9).

Для таких задач необходимым условием существования решения яв ляется разрешимость уравнения (2.9).

2.2.2 Достаточные условия Пусть z и (1, 2 ) удовлетворяют необходимым условиям существования решения. Выделим из набора ограничений–неравенств (2.4c) удержива ющие ограничения, т.е. такие, что hi () = 0, z i = 1,..., n2.

Введём функцию Лагранжа, содержащую только удерживающие огра ничения–неравенства, 1 L(z, 1, 2 ) = Qz, z f (), z + (), H H где n n1 m f () = f0 1i i 2i i, i=1 i= n n1 m = 0 1i ai 2i bi.

i=1 i= Известно (достаточные условия экстремума для задач с равенствами и неравенствами [3, стр. 293–294]), что если 2 0 и cуществует такое число 0, что z (2.10) Lzz (, 1, 2 )[, ] для любого, лежащего в подпространстве { H | gi (), = 0, i = 1,..., n1 m;

hi (), = 0, i = 1,..., n2 }, z z то z доставляет локальный минимум в задаче (2.4a)–(2.4c).

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

Вообще говоря, если в задаче присутствуют ограничения–неравенс тва, то множество допустимых решений может быть невыпуклым, и мо гут существовать несколько локальных минимумов. Однако, так как в рассматриваемом нами случае все функционалы в (2.1c) аффинные, то множество допустимых решений выпукло, как пересечение множе ства выпуклых множеств. Следовательно, локальный минимум задачи является и глобальным минимумом.

Если в пробной точке z, которая удовлетворяет необходмым услови ям существования решения, нет ни одного удерживающего ограничения, то в задаче (2.4a)–(2.4c) можно оставить только ограничения–равенства (2.4b).

Вторая производная функции Лагранжа является самосопряжённым оператором Q = I K L (H). Известно [16], что самосопряжённый оператор Q положительно определён:

Q, 0, H, H тогда и только тогда, когда его спектр (Q) не содержит отрицательных чисел: (Q) [0, ). Если оператор Q строго положительно определён и, кроме того, фредгольмов, то (Q) (0, ) [2, стр. 191]. Очевидно, что в случае (Q) (, ), 0, будет выполнено и условие (2.10).

Если K ограниченный самосопряжённый положительно определён ный компактный оператор, то для строгой положительной определённо сти Q = I K необходимо и достаточно выполнения условия K L (H) 1.

Пусть норма оператора K больше единицы. Тогда для проверки по ложительной определённости I K требуется оценить значение верхней границы спектра (K), так как, очевидно, выполнение max (K) влечёт (I K) (0, ).

2.3 Проекционный метод исследования Как было показано в предыдущем разделе, для проверки необходимого и достаточного условий существования решения требуется установить, является ли оператор I K обратимым и (строго) положительно опре делённым. К сожалению, даже в относительно простых задачах оператор K имеет достаточно сложный вид, что существенно затрудняет или даже делает невозможным исследование вручную.

Одним из методов приближённого решения сложных задач является замена исходной задачи некоторой близкой к ней задачей, имеющей более простой вид. В данной работе исходное уравнение z Kz = f заменяется некоторым приближённым z Kz = f, где оператор K, близкий по норме к K, является конечномерным.

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


2.3.1 Приближённое уравнение Выберем в гильбертовом пространстве H некоторую полную систему функций {i }, ортогональных с единичным весом:

i, i = ii ci, ci R, i = 1, 2,...

H Здесь ii символ Кронекера, ci = i 2, i = 1, 2,...

H Определим последовательность конечномерных пространств {Hn }, n = 1, 2,..., где Hn евклидово пространство, являющееся линейной оболочкой множества {1,..., n }. Обозначим через Pn ортопроектор H на Hn.

Элемент z пространства H может быть представлен в виде ряда z, i H z= i i, i H i= или, в случае ортонормированного базиса, z= z, i H i.

i= Сопоставив каждому z из H последовательность { z, i H }, получим изо морфизм между рассматриваемым пространством H и пространством l суммируемых с квадратом последовательностей.

Разложив Kz и f по базису {i }, получим Kz, i H z, j H Kj, i H Kz = i = i, i, i H j, j H i, i H i=1 i=1 j= f, i H f= i.

i, i H i= Запишем z Kz = f в виде равенства обобщённых многочленов z, i H z, j H Kj, i H f, i H i i = i.

i, i H j, j H i, i H i, i H i=1 i=1 j=1 i= Многочлены в правой и левой частях уравнения равны, если равны ко эффициенты при соответстующих функциях i :

z, i H z, i H Kj, i H f, i H (2.11) =, i = 1, 2,...

i, i H i, i H j, j H i, i H j= Обозначим z, i H Kj, i H f, i H zi =, Kij =, fi =.

i, i H j, j H i, i H Зададим последовательность конечномерных операторов {Kn }, где Kn отображает Hn в Hn, и рассмотрим систему n линейных алгебраических уравнений с n неизвестными n (n) (n) (2.12) zi Kij zj = fi, i = 1,..., n.

j= (n) Заметим, что, вообще говоря, zi = zi. Однако, z (n) z при n.

Таким образом мы заменили исходную бесконечную систему уравне ний (2.11) её конечномерным приближением (2.12).

Так как оператор K самосопряжённый, то Kij = Kji и матрица си стемы симметрична. Для обращения таких матриц существуют специ альные численные методы, позволяющие существенно сократить объём необходимых вычислений (см., например, [11, стр. 150]).

Все собственные числа симметричной матрицы вещественны. После вычисления характеристического многочлена det(E Kn ) их можно локализовать с любой степенью точности, используя, например, метода Штурма [22, стр. 42–43].

2.3.2 Проверка существования обратного оператора Пусть уравнение z (n) Kn z (n) = f (n) имеет решение. В общем случае, однако, это ещё не гарантирует обратимости исходного оператора I K.

Известно, что множество обратимых операторов открыто, т.е. свой ство обратимости является грубым, сохраняющимся при достаточно ма лых возмущениях. Если обратим оператор I Kn, то все операторы I K, для которых выполнено неравенство def (2.13) n K = K Kn =, L (H) L (H) (I Kn )1 L (H) также имеют обратные (теорема 4 [16, стр. 212]).

Если оператор I K обратим, то (теоретически) всегда можно найти такое приближение Kn, чтобы было выполнено (2.13).

Неравенство (2.13) не будет выполняется в двух случаях: когда опера тор I K не имеет обратного и когда полученное приближение Kn лежит слишком далеко от K в пространстве L (H). Если при попытке найти более точное приближение Kn Hn, n n неравенство (2.13) снова 1/ (I Kn )1 L (H), то вопрос, не будет выполнено, т.е. n K L (H) какой из двух причин это вызвано, останется открытым.

2.3.3 Проверка положительной определённости Для доказательства строгой положительной определённости оператора I K достаточно показать, что все собственные числа оператора K мень ше единицы.

Теорема 2.1. Пусть K компактный самосопряжённый оператор и µmin, µmax являются нижней и верхней границами его спектра соответ ственно. Тогда значения верхней и нижней границ µmin, µmax спектра возмущённого оператора K = K + K, где K компактный самосо пряжённый оператор, удовлетворяют неравенствам µmin K µmin µmin + K и µmax K µmax µmax + K.

Доказательство. Следует из теоремы Вейля [24, стр. 257] (соответству ющие собственные числа двух компактных самосопряжённых операто ров K и K отстоят друг от друга не более чем на K K ).

(n) Пусть µmax верхняя граница спектра оператора Kn (наибольшее собственное число). Тогда для µmax верхней границы спектра опера тора K выполнено неравенство (n) µmax µmax + n K L (H).

Таким образом, если µmax 1, то требуется подобрать такое приближе ние Kn, чтобы норма погрешности n K L (H) была меньше 1 µmax.

Если µmax 1 n K L (H), то оператор I K гарантированно не является положительно определённым.

Вариант µmax = 1 не рассматривается, так как в этом случае 1 явля ется собственным числом оператора K, и вычисления останавливаются ещё на этапе проверки необходимого условия существования решения.

2.3.4 Оценка точности решения Пусть оператор I K имеет обратный, и удалось найти такое прибли жение Kn, что выполнено (2.13). В этом случае справедлива следующая оценка нормы (I K)1 (теорема 4, [16, стр. 212]) (I Kn )1 L (H) (I K).

L (H) 1 (I Kn )1 L (H) · n K L (H) Пусть z (n) решение z Kn z = f (n), а z решение исходного урав нения z Kz = f. Тогда (I Kn )1 L (H) (n) zz H 1 (I Kn )1 L (H) · n K L (H) · z (n) + f f (n). (2.14) n K L (H) H H Заметим, что при проверке условий существования решения и при оценке точности требуется, чтобы оценки норм операторов (I Kn )1, n K и функции f fn были известны точно.

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

Глава Доказательный вычислительный эксперимент Исследование операторного уравнения z Kz = f существенно зависит от того, элементами какого пространства H являются z и f. В этой главе даётся описание наиболее часто встречающегося случая H = L2 [a, b]. Обобщение на случай H = L2 (), Rm, кратко рассматривается в последнем разделе.

3.1 Алгоритм доказательного вычислитель ного эксперимента 1. Для интегрального уравнения Эйлера строится модельное интегральное уравнение Фредгольма второго рода с вырожденным ядром.

2. Строится система линейных алгебраических уравнений (СЛАУ).

3. Вычисляется обратная матрица СЛАУ и оценивается норма обратного оператора I + R.

4. Если K K 1/ I + R, то возврат к первому шагу и вычисление более точного нового приближения.

5. Вычисление характеристического многочлена матрицы СЛАУ и оценка границ спектра конечномерного оператора.

6. Если max(i ) 1 + K K, то оператор не является положительно определённым.

7. Если 1 K K 1 + K K, то возврат к первому шагу max(i ) и вычисление более точного нового приближения.

8. Решение СЛАУ. Построение приближённого решения с гарантированной оценкой точности.

3.2 Приближённое интегральное уравнение В пространстве L2 [a, b] операторное уравнение z Kz = f () обычно является интегральным уравнением Фредгольма второго рода b (3.1) z(t) K(t, s)z(s) ds = f (t, ), a t b.

a Заменим его приближенным уравнением с вырожденным ядром b (3.2) z(t) K(t, s)z(s) ds = f (t, ), a t b, a где a, Q, а K : Q Q Q и f : Q Q вычислимые функции.

b Выберем a и так, чтобы a a b и a a a, b b.

b b b Для аппроксимации функции на заданной области можно использовать её приближённое представление в виде интерполяционных многочленов, сплай нов, рациональных функций. В данной работе приближение функции ищется в виде отрезка ряда Фурье по заданной системе ортогональных функций. Пре имущества такого подхода:

• если в процессе вычислений оказывается, что точность приближения недостаточна, то требуется вычислить только дополнительные члены ря да, не пересчитывая всё;

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

В качестве базиса аппроксимации наилучшим выбором в данном случае явля ются система ортогональных многочленов Лежандра или система ортогональ ных функций Радемахера–Уолша [19, стр. 395 и 404]:

• элементы этих систем попарно ортогональны с единичным весом;

• коэффициенты многочленов Лежандра и узлы перемены знака функций Радемахера–Уолша являются рациональными числами.

3.2.1 Аппроксимация правой части Параметризованная правая часть уравнения (3.1) имеет вид L f (t, ) = f0 (t) + l fl (t), l= где L = n1 m + n2 (см. (2.4b) и (2.4c)).

Введём на [, сетку с координатами i Q, i = 0, 1,..., n:

a b] a = 0 1 · · · n = b, и будем аппроксимировать fl, l = 0, 1,..., L, отдельно на каждом отрезке [i1, i ], i = 1,..., n.

Пусть на [, ],, Q, задана полная система компьютерно–вычислимых функций {j }, ортогональных с единичным весом:

c Q, j = j, j j (t)j (t) dt = 0, j=j.

Введём функцию i, аффинно отображающую отрезок [i1, i ] на [, ]:

t i i (t) = + ( ).

i+1 i Множество суперпозиций S = {j i } по построению будет полной орто j= гональной системой (базисом) на [i1, i ], i = 1,..., n.

Разложим fl в ряд Фурье на [i1, i ] по системе S :

(3.3) fli (t) = flij j i (t), j= где 1 fl i, j L2 [,] 1 (3.4) flij = = fl i (t) j (t) dt.

j, j L2 [,] cj В общем случае значение интеграла можно найти только приближённо, поэто му предположим, что нам известно только приближённое значение flij. Обо (m ) значим через fli i отрезок ряда Фурье с mi слагаемыми mi (m ) fli i (t) = fij j i (t).

j= Введём оценку погрешности аппроксимации на i-м отрезке (m ) fl (t) fli i (t), i fl max i fl IQ.


t[i1,i ] Вычисление минимального значения i fl (т.е. вычисление maxt[i1,i ] fl (t) (m ) f i (t) ) является, вообще говоря, сложной задачей. Для наших целей, однако, li достаточно относительно грубой оценки этого значения. Единственное условие эта оценка должна быть гарантирована.

Пусть 1 1 (например, {i } система ортогональных многочленов Ле жандра [19, стр. 397] или система ортогональных функций Радемахера–Уолша [19, стр. 405]). В этом случае f (t) f(t), где mi n fl (t) = flij j i (t) [i1,i ] (t), i=1 j= при этом коэффициенты flij имеют вид [f i f, f + i f ], j = 1, lij l lij l = flij f, j 1.

lij Обозначим через fl (t) середину интервала fl (t):

fl (t) = mid fl (t).

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

Разложение K в ряд Фурье на (i, j)-й ячейке сетки (область [i1 1, i1 ] [i2 1, i2 ]) имеет вид (3.5) Ki1 i2 (t) = Ki1 i2 j1 j2 j1 i1 (t) j2 i2 (t), j1 =1 j2 = где Ki1 i2 j1 j2 = c1 c1 1 (3.6) K i1 (t), i2 (s) j1 (t)j1 (s) dt ds.

j1 j Здесь также обычно можно получить только приближённые значения Ki1 i2 j1 j2.

(mi,mi ) Обозначим через K 1 2 отрезок ряда Фурье с mi mi слагаемыми i1 i2 1 mi1 mi (mi mi2 ) (t) Ki1 i2 1 = Ki1 i2 j1 j2 j1 i1 (t) j2 i2 (t).

j1 =1 j2 = Оценим погрешность на области [i1 1, i1 ] [i2 1, i2 ]:

(mi mi2 ) |K(t, s) Ki1 i2 i1 i2 K max (t, s)|, i1 i2 K IQ.

t[i1 1,i1 ] s[i2 1,i2 ] В случае 1 1 опять определим интервальную функцию mi1 mi n n K(t) = Ki1 i2 j1 j2 j1 i1 (t) j2 i2 (s) [i1 1,i1 ] (t)[i2 1,i2 ] (s), i1 =1 i2 =1 j1 =1 j1 = где [Ki i j j i i K, Ki i j j + i i K], j1 = 1 и j2 = 1, 1212 12 1212 = Ki1 i2 j1 j i i j j, j1 1 или j2 1.

K Через K(t, s) снова обозначим середину интервала K(t, s):

K(t, s) = mid K(t, s).

3.2.3 Оценка погрешности аппроксимации По построению, при t (i1 1, i1 ), i1 = 1,..., n, и s (i2 1, i2 ), i1 = 1,..., n, |K(t, s) K(t, s)| i1 i2 K, |f (t) fl (t)| i1 fl, l = 0, 1,..., L.

Пусть K(t, s) 0 при (t, s) [, [, и fl (t) 0 при t a или t Тогда a b] a b] b.

b b K(t, s) K(t, s)2 dt ds K K = L (L2 [a,b]) a a n n i1 i2 K 2 (i1 i1 1 )(i2 i2 1 ).

K(t, s) dt ds + [a,b][a,b]\[, a, a b][ b] i1 =1 i2 = Введя оценку сверху для значения |K(t, s)| на [a, b] [a, b] \ [, [, a b] a b] 0 K max |K(t, s)|, (t,s)[a,b][a,b] (t,s)[, a, a b][ b] получим 0 K 2 · (b a)2 ( a)2 + K K b L (L2 [a,b]) n n i1 i2 K 2 · (i2 i2 1 ).

+ (i1 i1 1 ) · i1 =1 i2 = По аналогии, оценка погрешности аппроксимации для f (t) имеет вид n · (b + (a a) + fl 2 2 [a,b] 0 fl2 i fl2 · (i i1 ).

fl b) L i= 3.2.4 Решение приближённого уравнения Для решения полученного интегрального уравнение Фредгольма второго рода b z(t) K(t, s)z(s) ds = f (t, ) a с вырожденным ядром mi1 mi n n K(t, s) = Ki1 i2 j1 j2 j1 i1 (t) j2 i2 (s) [i1 1,i1 ] (t)[i2 1,i2 ] (s) i1 =1 i2 =1 j1 =1 j2 = и правой частью L f (t, ) = f0 (t) + l fl (t), l= где mi n fl (t) = flij j i (t) [i1,i ] (t), l = 0, 1,..., L, i=1 j= применим стандартный метод преобразование в систему линейных алгебра ических уравнений (СЛАУ) L (3.7) AZ = B0 + l Bl.

l= Перепишем K(t, s) в виде m n K(t, s) = j1 i1 (t) [i1 1,i1 ] (t) i1 =1 j1 = m n Ki1 i2 j1 j2 j2 i2 (s) [i2 1,i2 ] (s).

i2 =1 j2 = n Очевидно, что СЛАУ будет иметь N = уравнений с N неизвестными.

i=1 mi При сведении интегрального уравнения к СЛАУ вырожденное ядро K(t, s) обычно записывается в виде N K(t, s) = vi (t)ui (s).

i= Для преобразования двойного суммирования по i1 и j1 в суммирование по i, введём функции : N N (определяет номер отрезка) : {1,..., N } {1,..., n} и : N N (определяет номер базисной функции) : {1,..., N } {1,..., max mi } i Пара [, ] : N N N задаёт взаимно однозначное соответствие между номе рами уравнений (неизвестных) и номерами отрезков (i1 ) и базисных функций (j1 ):

уравнение: i 1 2 3... m1 1 m 1 m 1 + 1 m 1 + 2... N отрезок: (i) 1 1 1... 1 1 2 2... n функция: (i) 1 2 3... m1 1 m 1 1 2... mn Если m1 = m2 = · · · = mn = m, то i (i) = 1 +, m (i) = 1 + (i 1) mod m, i = 1,..., N.

Перепишем K(t, s) и fl (t), l = 0, 1,..., L, с использованием функций и :

N K(t, s) = (i) (i) (t) [(i)1,(i) ] (t) i= ui (t) N K(i)(k)(i)(k) (k) (k) (s) [(k)1,(k) ] (s), k= vi (s) N fl (t) = fl(i)(i) (i) (i) (t) [(i)1,(i) ] (t).

i= Используя стандартную процедуру перехода к СЛАУ, получаем следующее выражение для элементов матрицы коэффициентов:

b Aij = ij uj (t)vi (t) dt = a N b = ij K(i)(k)(i)(k) (j) (j) (t) (k) (k) (t) a k= [(j)1,(j)] (t)[(k)1,(k) ] (t) dt Удаляя произведения характеристических функций по непересекающимся от резкам, получаем (j) Aij = ij K(i)(j)(i)(j) (j) (j) (t) dt (j) или, с учётом произведения ортогональных функций, (j) (j) (3.8) Aij = ij K(i)(j)(i)(j) · c(j).

Аналогично, выражение для элементов правой части СЛАУ имеет вид:

N N b Bli = fl (t)vi (t) dt = fl(k)(k) K(i)(k )(i)(k ) a k=1 k = b (k) (k) (t) (k ) (k ) (t) [(k)1,(k) ] (t)[(k )1,(k ) ] (t) dt = a N (k) = fl(k)(k) K(i)(k)(i)(k) (k) (k) (t) dt = (k) k= N (k) (k). (3.9) = fl(k)(k) K(i)(k)(i)(k) · c(k) k= Заметим, что при определении вида Aij и Bi существенно использовалась ортогональность базисных функций. Если бы базисные функции были не ор тогональны, то при вычислении, например, Aij требовалось бы рассматривать интегралы от всевозможных комбинаций (j) (j) (t) (k) (k) (t), k = 1,..., N, где (j) = (k). Кроме того, использование ортогональных функций позволяет в дальнейшем при необходимости повысить точность аппроксимации просто вычисляя дополнительные элементы ряда Фурье (переходя от пространства Hm к Hm с m m).

По построению все элементы матрицы A и вектора B принадлежат мно жеству Q. Используя системы компьютерной алгебры, например программы Maple или Maxima, можно (правда, иногда за довольно длительное время) вычислить обратную матрицу A1 и получить (точное) решение Z = A1 B.

Применение специализированного программного обеспечения, использующе го специфические возможности компьютерной архитектуры, позволяет суще ственно сократить время вычислений. Например, при программировании на C++ хорошие результаты даёт использование библиотеки процедур GNU MP.

Решение системы линейных алгебраических уравнений (3.7) имеет вид L l A1 Bl.

Z=A B0 + l= Z (0) Z (l) Решение интегрального уравнения (3.2) можно записать в виде N z(t, ) = f (t, ) + ui (t)Zi.

i= или, с учётом того, что f(t, ) также является отрезком ряда по ui, N L (0) (l) (3.10) z(t, ) = ui (t) f0i + Zi + l (fli + Zi ).

i=1 l= Подставим (3.10) в выражения для gi (z) (см. (2.4b)) и hi (z) ((2.4c)). Для опре деления значений множителей Лагранжа l, l = 1,..., n1 m + n2, требуется решить систему уравнений:

gi (z) = 0, i = 1,..., n1 m, (3.11) n1 m+i hi (z) = 0, i = 1,..., n2.

При этом n1 m+i, i = 1,..., n2, должны иметь неотрицательные значения.

После решения системы (3.11) и подстановки в (3.10), получим искомое решение приближённого интегрального уравнения:

N z(t) = ui (t)(fi + Zi ), i= где L L (0) (l) fi = f0i + l fli, Zi = Zi + l Zi.

l=1 l= После того, как оценки значений множителей Лагранжа будут найдены, решение z(t) можно записать в виде N N N A1 Bj = z(t) = f (t) + ui (t)Zi = f (t) + ui (t) ij i=1 i=1 j= N N b A1 ui (t)vj (s)f (s) ds.

= f(t) + ij a i=1 j= Следовательно, ядро резольвентного оператора R (I + R = (I K)1 ) имеет вид N N A1 ui (t)vj (s). (3.12) R(t, s) = ij i=1 j= Определим решение интервального уравнения b (3.13) z(t) K(t, s)z(t) dt = f(t), a как множество решений всевозможных уравнений b z(t) K(t, s)z(s) dt = f (t), a где K(t, s) K(t, s) и f (t) f(t, ). Границы K(t, s) и K(t, s) ядра K(t, s) яв ляются по построению вырожденными ядрами, поэтому в процессе решения уравнение (3.13) сводится к системе линейных уравнений с интервальными коэффициентами AZ = B, где (j) (j) Aij = ij K(i)(j)(i)(j) · c(j), N (k) (k) Bi = f(k)(k) K(i)(k)(i)(k) · c(k).

k= При решении системы линейных алгебраических уравнений с интервальными элементами можно использовать, например, интервальную модификацию ме тода исключения Гаусса [5]. Другой подход применение следующей теоремы:

Теорема 3.1 ( [11, стр. 84]). Предположим, что A Rnn, 0 = b Rn, Ax = b, A Rnn, b Rn (A + A)y = b + b, и что |A| |A|, |b| |b|. Если (A) = r 1, то A + A невырож денна, и yx |A1 | · |A|.

x 1r Здесь |X| матрица, состоящая из абсолютных значений матрицы X.

1 Введём обозначения A = wid A и B = wid B. Тогда 2 |Aij | |Bi | = max max, max |Aij | |Bi | i,j i и · A r=· A.

Пусть r 1. Тогда разница Z = Z Z между любым решением Z и най денным приближённым решением Z = A1 B может быть оценена по норме как 2 |A| · |A1 | z = Z Z.

1r Элементы интервального вектора Z имеют вид Zi = [Zi z, Zi + z ].

Таким образом решение интегрального уравнения zKz = f лежит внутри интервального решения N z(t) = f(t) + Zi (i) (i) (t) [(i)1 (t),(i) ] (t).

i= Середину интервальной оценки можно принять за приближённое решение z (t) = mid z(t).

3.2.5 Оценка погрешности приближённого решения Обозначим обратный к I K оператор через I + R. Оценку нормы разности между найденным приближением z и истинным решением z определим по формуле I +R L (L2 [a,b]) zz L2 [a,b] 1 I +R · K K L (L2 [a,b]) L (L2 [a,b]) K K ·z + f f.

L (L2 [a,b]) L2 [a,b] L2 [a,b] При вычислении оценок для z и I + R будем полагать, что z (t) 0 при t [, и R(t, s) 0 при (t, s) [, [, a b] a b] a b].

Оценим сверху норму резольвенты I + R:

b b 1/ R(t, s)2 dt ds I +R 1+.

L (L2 [a,b]) a a Резольвентное ядро R(t, s) имеет вид (см. (3.12)) N N N A1 K(j)(k)(j)(k) R(t, s) = ij i=1 j=1 k= (i) (i) (t) (k) (k) (s) [(i)1,(i) ] (t)[(k)1,(k) ] (s), а квадрат ядра (с учётом произведений характеристических функций):

N N N N A1 A1 K(j1 )(k)(j1 )(k) R(t, s) = ij1 ij j1 =1 j2 =1 i=1 k= 2 K(j2 )(k)(j2 )(k) (i) (i) (t) (k) (k) (s) [(i)1,(i) ] (t)[(k)1,(k) ] (s).

Следовательно, N N N N b b A1 A1 K(j1 )(k)(j1 )(k) R(t, s)2 dt ds = ij1 ij ( ) a a j1 =1 j2 =1 i=1 k= K(j2 )(k)(j2 )(k) c(i) c(k) ((i) (i)1 )((k) (k)1 ). (3.14) При больших размерностях матрицы A элементы обратной матрицы A рациональные дроби могут иметь очень большую длину числителя и знаме нателя. Поэтому, чтобы выполнение N 4 сложений происходило быстрее, можно заменить все Ki1 i2 j1 j2 интервальными оценками из IF и использовать интер вальную арифметику, программная реализация которой основана на аппарат ной поддержке вычислений с направленным округлением.

3.2.6 Исследование положительной определённости оператора Если верхняя граница спектра самосопряжённого оператора K не превосходит единицы, то I K положительно определён. Если верхняя граница спектра строго меньше едицины, то I K строго положительно определён. Для ком пактных симметричных операторов I K и I K значения их верхних и нижних границ отличаются друг от друга не более чем на K K L (L [a,b]), поэтому, найдя достаточно точную оценку верхней границы спектра оператора K, можно определить, будет ли I K положительно определённым. Данный метод не будет работать, если спектр K содержит единицу, однако, в этом случае оператор I K не имеет обратного, поэтому процесс исследования за дачи остановится ещё на этапе проверки необходимого условия существования решения.

Спектральное множество конечномерного оператора K совпадает с множе ством собственных чисел матрицы A = E A.

Спектральное множество самосопряжённого оператора состоит из веще ственных чисел.

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

Если + + K K 1, L (L2 [a,b]) то I K 0. При K K L (L2 [a,b]) оператор I K гарантированно не будет положительно определённым. Случай K K 1 + + K K L (L2 [a,b]) L (L2 [a,b]) оставляет ситуацию неопределённой и требует дополнительных исследований.

Существует много методов вычисления коэффициентов характеристиче ского многочлена [11, 26]. Надёжным, хотя и довольно трудоёмким является видоизменение Д. К. Фаддеева метода Леверье [26, стр. 311–316]. В нём для матрицы A вычисляется последовательность матриц A1,..., AN :

A1 = A Sp A1 = q1 B1 = A1 q1 E A2 = AB1 (Sp A2 )/2 = q2 B2 = A2 q2 E...

...

...

AN 1 = ABN 2 (Sp AN 1 )/(N 1) = qN 1 BN 1 = AN 1 qN 1 E AN = ABN 1 (Sp AN )/N = qN BN 1 = AN qN E и характеристический многочлен имеет вид (1)n ( N q1 N 1 q2 N 2 · · · qN ).

Если A неособенная матрица, то BN A1 =.

qN По построению элементы A рациональные числа, следовательно харак теристический многочлен det(A E) будет вычислимой функцией. Так как требуется оценить только верхнюю границу спектра, нет необходимости вы числять все корни характеристического многочлена.

Для оценки значения наибольшего корня будем использовать решение урав нения методом Ньютона.

Известно [22, стр. 10, Теорема 1.2], что все корни многочлена p() = N q1 N 1 q2 N 2 · · · qN лежат внутри круга || = 1 + max |qi | i или, в случае когда все корни вещественные числа, внутри интервала [1 max |qi |, 1 + max |qi |].

i i В качестве исходного приближения (0) возьмём правую границу этого интер вала. Следующие приближения определяются по формуле метода Ньютона p( i ) i+1 = (i), i = 0, 1, 2,...

p ( i ) Так как в методе Ньютона исходная функция p заменяется касательной в точке приближения, последовательность (0), (1), (2),...

сходится к наибольшему собственному числу сверху. Если на каком-то шаге будет выполнено условие (i) + K K 1, L (L2 [a,b]) то оператор I K положительно определённый.

Другой способ локализации наибольшего корня метод Штурма [22, стр. 42– 43]: строится последовательность Штурма и, например с помощью метода де ления отрезка пополам определяются интервалы, на каждом из которых на ходится в точности по одному корнью многочлена.

Однако при построении последовательности Штурма возникают проблемы.

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

3.3 Обобщение на случай функций многих переменных Рассмотрим изменения в описанной схеме исследования в случае H = L2 (), заданный прямоугольный гиперпараллелепипед в Rd :

где = [a1, b1 ] · · · [ad, bd ] Rd.

Интегральное уравнение Фредгольма второго рода в этом случае имеет вид b1 bd x(t1,..., td )... K(t1,..., td, s1,..., sd )x(s1,..., sd ) ds1,... dsd = a1 ad = f (t1,..., td ), (t1,..., td ).

3.3.1 Построение приближённого уравнения Заменим во всех отрезках ai и bi, i = 1,..., d на ai, i Q так, чтобы b ai i ai b bi.

d a Обозначим через i=1 [i, bi ] На каждом отрезке [i, i ], i = 1,..., d зададим сетку (i) :

ab (i) (i) (i) ai = 0 1 · · · ni = i, (i) b j Q, i = 1,..., d, j = 0,..., ni.

Выберем полную систему функций {k }, ортогональных на [, ] с ве (i) (i) (i) сом 1. Пусть аффинные функции j отображают [j1, j ] на [, ]. Тогда приближения для ядра уравнения и его правой части можно записать в виде K(t1,..., td, s1,..., sd ) K(t1,..., td, s1,..., sd ) = = Ki11...i1d i21...i2d j11...j1d j21...j2d i11 =1,...,n1 i21 =1,...,n1 j11 =1,...,mi11 j21 =1,...,mi............

i1d =1,...,nd i2d =1,...,nd j1d =1,...,mi1d j2d =1,...,mi2d d d (k) (l) (3.15) j1k i1k (tk ) [i,i1k ] (tk ) j1l i1l (sl ) [i,i1l ] (sl ) 1k 1 1l k=1 l= и f (t1,..., td ) f (t1,..., td ) = d (k) (3.16) = fi1...id j1...jd jk ik (tk ) [i,ik ] (tk ) k i1 =1,...,n1 j1 =1,...,mi1 k=......

id =1,...,nd jd =1,...,mid соответственно. Здесь Ki11...i1d i21...i2d j11...j1d j21...j2d = d c1 c1 1 1 1 = ·... K i11 (t1 ),..., i1d (td ), i21 (s1 ),..., i2d (sd ) j1i j2i i= 2d j11 (t1 )... j1d (td )j21 (s1 )... j2d (sd ) dt1... dtd ds1... dsd и d c1 · 1 fi1...id j1...jd =... f i1 (t1 ),..., id (td ) ji i= d j1 (t1 )... jd (td ) dt1... dtd.

3.3.2 Построение системы линейных алгебраических уравнений Заметим, что даже при отосительно небольших значениях d размерность мат рицы A может быть очень большой:

nd n N= ··· mi1... mid.

i1 =1 id = Введём функцию : N N2d, которая является биективным отображением множества {1,..., N } на множество всевозможных векторов (i1,..., id, j1,..., jd ), ik = 1,..., ni, jk = 1,..., mik, k = 1,..., d.

Для сокращения записи будем обозначать список (i)k, (i)k+1,..., (i)l че рез (i)k,...,l.

Перепишем (3.15) в виде N K(t1,..., td, s1,..., sd ) = ui (t1,..., td )vj (s1,..., sd ).

i= Здесь d (k) ui (t1,..., td ) = (i)d+k (i) (tk ) [(i),(i)k ] (tk ) k k k= и N vi (s1,..., sd ) = K(i)1,...,d (l)1,...,d (i)d+1,...,2d (l)d+1,...,2d l= d (k) (l)d+k (l)k (tk ) [(l),(l)k ] (tk ).

k k= С помощью стандартной процедуры перехода к СЛАУ, получим следующее выражение для матрицы коэффициентов:

Aij = ij uj (t1,..., td )vi (t1,..., td ) dt1... dtd = N = ij K(i)1,...,d (l)1,...,d (i)d+1,...,2d (l)d+1,...,2d l= d (k ) (l)d+k1 (l)k (tk1 ) [(l),(l)k ] (tk1 ) k1 1 k1 = d (k ) (j)d+k (j) (tk2 ) [(j),(j)k ] (tk2 ) dt1... dtd = k2 2 k2 k2 = (k) (k) d (j)k (j)k = ij K(i)1,...,d (j)1,...,d (i)d+1,...,2d (j)d+1,...,2d c(j)d+k.

k= Как и в случае = [a, b], выбор в качестве базиса ортогональных функций существенно упрощает формулу.

Записав (3.16), правую часть интегрального уравнения, в виде N d (k) f (t1,..., td ) = f(i) (j)d+k (j)k (tk ) [(j),(j)k ] (tk ) k j=1 k= и применив стандартную процедуру перехода, получим следующее выражение для элементов правой части СЛАУ Bi = f (t1,..., td )vi (t1,..., td ) dt1... dtd = N N = f(j) K(i)1,dots,d (l)1,...,d (i)d+1,...,2d (l)d+1,...,2d l=1 j= d (k ) (l)d+k1 (l)k (tk1 ) [(l),(l)k ] (tk1 ) k1 1 k1 = d (k ) (j)d+k2 (j)k (tk2 ) [(j),(j)k ] (tk2 ) dt1... dtd k2 2 k2 = или, после упрощения, (k) (k) N d (j) (j) k k Bi = f(j) K(i)1,...,d (j)1,...,d (i)d+1,...,2d (j)d+1,...,2d c(j)d+k.

j=1 k= Дальнейшие вычисления обращение матрицы, оценка границы спектра, решение СЛАУ и т.п. аналогичны случаю одной переменной (H = L2 [a, b]).

3.3.3 Оценка погрешности приближённого решения Пусть приближённое решение найдено. При вычислении оценки его погрешно сти основную сложность составляет оценка нормы резольвентного оператора R.

Для оценки сверху нормы резольвенты вычислим правую часть неравен ства 1/ I +R 1+ R(t1,..., td, s1,..., sd ) dt1... dtd ds1... dsd L (L2 ()) Резольвентное ядро R(t, s) имеет вид R(t1,..., td, s1,..., sd ) = N N N A1 K(j)1,...,d (k)1,...,d (j)d+1,...,2d (k)d+1,...,2d = ij i=1 j=1 k= d (l ) (i)d+l (i) (tl1 ) [(i),(i)l ] (tl1 ) l1 1 l1 l1 = d (l ) (k)d+l2 (k)l (sl2 ) [(k),(k)l ] (sl2 ).

l2 2 l1 = Его квадрат (с учётом произведения характеристических функций):

R(t1,..., td, s1,..., sd )2 = N N N N A1 A1s K(j1 )1,...,d (k)1,...,d (j1 )d+1,...,2d (k)d+1,...,2d = ij1 ij j1 =1 j2 =1 j=1 k= K(j1 )1,...,d (k)1,...,d (j1 )d+1,...,2d (k)d+1,...,2d d (l ) (i)d+l1 (i)l (tl1 ) [(i),(i)l ] (tl1 ) l1 1 l1 = d (l ) (k)d+l (k)l (sl2 ) [(k),(k)l ] (sl2 ).

l2 2 2 l1 = Проинтегрировав R(t1,..., td, s1,..., sd )2 по области, получим R(t1,..., td, s1,..., sd )2 dt1... dtd ds1... dsd = N N N N A1 A1 K(j1 )1,...,d (k)1,...,d (j1 )d+1,...,2d (k)d+1,...,2d = ij1 ij j1 =1 j2 =1 j=1 k= K(j1 )1,...,d (k)1,...,d (j1 )d+1,...,2d (k)d+1,...,2d d c(i)d+l c(k)d+l (l) (l) (l) (l) (i)l (i)l1 (k)l (k)l1.



Pages:   || 2 | 3 |
 





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

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