Численное решение интегральных уравнений методом квадратур

20. Численное интегрирование методом квадратур

Численное интегрирование методом квадратур

Приведенные ниже функции осуществляют интегрирование и двойное интегрирование, используя квадратурную формулу Симпсона или метод Гаусса-Лобатто. Квадратура — численный метод нахождения площади под графиком функции/(т), т. е. вычисление определенного интеграла вида

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

Функции quad и quadl используют два различных алгоритма квадратуры для вычисления определенного интеграла. Функция quad выполняет интегрирование по методу низкого порядка, используя рекурсивное правило Симпсона [4, 40]. Но она может быть более эффективной при негладких подынтегральных функциях или при низкой требуемой точности вычислений. Алгоритм quad в MATLAB 6 изменен по сравнению с предшествовавшими версиями, точность по умолчанию по сравнению с версиями 5.3х повышена в 1000 раз (с 10- 3 до 10- 6 ). Новая функция quadl (квадратура Лобатто) использует адаптивное правило квадратуры Гаусса— Лобатто очень высокого порядка. Устаревшая функция quads выполняла интегрирование, используя квадратурные формулы Ньютона—Котеса 8-го порядка [40]. Достижимая точность интегрирования гладких функций в MATLAB 6 поэтому также значительно выше, чем в предшествующих версиях.

quad(@fun,a.b) — возвращает численное значение определенного интеграла от заданной функции @fun на отрезке [а Ь]. Используется значительно усовершенствованный в MATLAB 6 адаптивный метод Симпсона;

quad(@fun,a.b.tol) — возвращает численное значение определенного интеграла с заданной относительной погрешностью tol. По умолчанию to1=l.e-6. Можно также использовать вектор, состоящий из двух элементов tol =[rel_tol abs_tol], чтобы точно определить комбинацию относительной и абсолютной погрешности;

quad(@fun,a.b,tol .trace) — возвращает численное значение определенного интеграла и при значении trace, не равном нулю, строит график, показывающий ход вычисления интеграла;

quad(@fun,a,b.tol,trace,PI,P2. ) — возвращает численное значение определенного интеграла по хот подынтегральной функции fun, использует дополнительные аргументы P1, P2. которые напрямую передаются в подынтегральную функцию: G=fun(X.Pl,P2. ). Примеры:

dblquad(@fun,inmin,inmax.outmin,outmax) — вычисляет и возвращает значение двойного интеграла для подынтегральной функции fun (Inner, outer), по умолчанию используя квадратурную функцию quad. Inner — внутренняя переменная, изменяющаяся на закрытом интервале от inmin до inmax, a outer — внешняя переменная, изменяющаяся на закрытом интервале от outmin до outmax. Первый аргумент @fun — строка, описывающая подынтегральную функцию. Этс может быть либо дескриптор функции, либо объект inline (в последнем случае символ «@» в ее записи отсутствует). Обычная запись в апострофах тепер недопустима. Эта функция должна быть функцией двух переменных вид. fout=fun(inner.outer). Функция должна брать вектор inner и скаляр outer возвращать вектор fout, который является функцией, вычисленной в outer и каждом значении inner [ Функция inime(‘expr’, ‘argl’. ‘argn’) так же создает объект, но без дескриптора, ‘ехрг’ — выражеши Строки ‘argx’ —входные аргументы. При их отсутствии по умолчанию подставляется х. Если вмест ‘arg’ — скаляр, то он означает количество дополнительных переменных Р. Примеры запис; g = inline(exp); g — inline(‘t A 2’); g = inline(‘sin(2*pi*f + theta)’); g = inline(‘sin(2*pi*f + theta)’, : ‘theta’); g — inline(‘x A Pl+P2’, 2). — Примеч. ред. ];

dblquad(@fun,inmin.inmax.outmin,oiitmax,tol .trace) — передает в функцию dblquad параметры tol и trace. Смотрите справку по функции quad для получения информации о параметрах to! и trace;

dblquad(@fun,inmin,inmax,outmin,outmax.tol .trace,order) — передает параметры tol и trace для функции quad или quadl в зависимости от значения строки order. Допустимые значения для параметра order — @quad , @quadl или имя любого определенного пользователем квадратурного метода с таким же вызовом и такими же возвращаемыми параметрами, как у функций quad и quad!. (Например, при проверке старых программ можно использовать @quad8 для большей совместимости с прежними версиями MATLAB). По умолчанию (без параметра order) вызывается @quad. поскольку подинтегральные функции могут быть негладкими.

Пример: пусть m-файл integl.m описывает функцию 2*y*sin(x)+x/2*cos(y), тогда вычислить двойной интеграл от той функции можно следующим образом:

Вычисление определённых интегралов: базовые алгоритмы

В этой публикации описаны простейшие методы вычисления интегралов функций от одной переменной на отрезке, также называемые квадратурными формулами. Обычно эти методы реализованы в стандартных математических библиотеках, таких как GNU Scientific Library для C, SciPy для Python и других. Публикация имеет целью продемонстрировать, как эти методы работают «под капотом», и обратить внимание на некоторые вопросы точности и производительности алгоритмов. Также хотелось бы отметить связь квадратурных формул и методов численного интегрирования обыкновенных дифференциальных уравнений, о которых хочу написать ещё одну публикацию.

Определение интеграла

Интегралом (по Риману) от функции на отрезке называется следующий предел:

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

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

Более наглядно геометрическое определение — интеграл равен площади криволинейной трапеции, ограниченной осью 0x, графиком функции и прямыми x = a и x = b (закрашенная область на рисунке).

Квадратурные формулы

Определение интеграла (1) можно переписать в виде

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

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

Вычислительная задача

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

Требуется найти приближённое значение интеграла .
Решения будут реализованы на языке Python 3.6.

Для проверки методов используется интеграл .

Кусочно-постоянная аппроксимация

Идейно простейшие квадратурные формулы возникают из применения выражения (1) «в лоб»:

Т.к. от метода разбиения отрезка точками и выбора точек значение предела не зависит, то выберем их так, чтобы они удобно вычислялись — например, разбиение возьмём равномерным, а для точек вычисления функции рассмотрим варианты: 1) ; 2) ; 3) .

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

Для анализа производительности квадратурных формул построим график погрешности в координатах «число точек — отличие численного результата от точного».

Что можно заметить:

1. Формула со средней точкой гораздо точнее, чем с правой или левой точками
2. Погрешность формулы со средней точкой падает быстрее, чем у двух остальных
3. При очень мелком разбиении погрешность формулы со средней точкой начинает возрастать
   Первые два пункта связаны с тем, что формула прямоугольников со средней точкой имеет второй порядок аппроксимации, т.е. , а формулы правых и левых прямоугольников — первый порядок, т.е. .
   Возрастание погрешности при измельчении шага интегрирования связано с нарастанием погрешности округления при суммировании большого числа слагаемых. Эта ошибка растёт как , что не даёт при интегрировании достигнуть машинной точности.
   Вывод: методы прямоугольников с правой и левой точками имеют низкую точность, которая к тому же медленно растёт с измельчением разбиения. Поэтому они имеют смысл разве что в демонстрационных целях. Метод прямоугольников со средней точкой имеет более высокий порядок аппроксимации, что даёт ему шансы на использование в реальных приложениях (об этом чуть ниже).

Кусочно-линейная аппроксимация

Следующий логический шаг — аппроксимировать интегрируемую функцию на каждом из подотрезков линейной функцией, что даёт квадратурную формулу трапеций:

Иллюстрация метода трапеций для n=1 и n=2.

В случае равномерной сетки длины всех отрезков разбиения равны, и формула имеет вид

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

Контроль точности вычисления

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

Как это реализовать? Один из простых методов оценки погрешности — правило Рунге — разность значений интегралов, рассчитанных по n и 2n точкам, даёт оценку погрешности: . Метод трапеций удобнее для удвоения мелкости разбиения, чем метод прямоугольников с центральной точкой. При расчёте методом трапеций для удвоения числа точек нужны новые значения функции только в серединах отрезков предыдущего разбиения, т.е. предыдущее приближение интеграла можно использовать для вычисления следующего.

Метод прямоугольников не требует вычислять значения функции на концах отрезка. Это означает, что его можно использовать для функций, имеющих на краях отрезка интегрируемые особенности (например, sinx/x или x -1/2 от 0 до 1). Поэтому показанный далее метод экстраполяции будет работать точно так же и для метода прямоугольников. Отличие от метода трапеций лишь в том, что при уменьшении шага вдвое отбрасывается результат предыдущих вычислений, однако можно утроить число точек, и тогда предыдущее значение интеграла также можно использовать для вычисления нового. Формулы для экстраполяции в этом случае необходимо скорректировать на другое соотношение шагов интегрирования.

Отсюда получаем следующий код для метода трапеций с контролем точности:

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

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

Кусочно-параболическая аппроксимация

Следующим шагом аппроксимируем функцию элементами парабол. Для этого требуется, чтобы число отрезков разбиения было чётным, тогда параболы могут быть проведены через тройки точек с абсциссами <(x₀=a, x₁, x₂), (x₂, x₃, x₄), . (x_n-2, x_n-1, x_n=b)>.

Иллюстрация кусочно-параболического приближения на 3 и 5 точках (n=2 и n=3).

Приближая интеграл от функции на каждом из отрезков [x_k;x_k+2] интегралом от параболической аппроксимации на этом отрезке и считая точки равномерно распределенными (x_k+1=x_k+h), получаем формулу Симпсона:

Из формулы (4) напрямую получается «наивная» реализация метода Симпсона:

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

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

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

Сравним эффективность метода трапеций и парабол:

Как видим, обоими методами ответ можно получть с достаточно высокой точностью, но количество вызовов подынтегральной функции разительно отличается — метод более высокого порядка эффективнее в 32 раза!

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

Отсюда и возникает такой рост эффективности по сравнению с простой формулой трапеций.

Что дальше?

Дальнейшая логика повышения точности квадратурных формул, в целом, понятна — если функцию продолжать приближать многочленами всё более высокой степени, то и интеграл от этих многочленов будет всё точнее приближать интеграл от исходной функции. Этот подход называется построением квадратурных формул Ньютона-Котеса. Известны формулы вплоть до 8 порядка аппроксимации, но выше среди весовых коэффициентов w_i в (2) появляются знакопеременные члены, и формулы при вычислениях теряют устойчивость.

Попробуем пойти другим путём. Ошибка квадратурной формулы представляется в виде ряда по степеням шага интегрирования h. Замечательное свойство метода трапеций (и прямоугольников со средней точкой!) в том, что для неё этот ряд состоит только из чётных степеней:

На нахождении последовательных приближений к этому разложению основана экстраполяция Ричардсона: вместо того, чтобы приближать подынтегральную функцию многочленом, по рассчитанным приближениям интеграла строится полиномиальная аппроксимация, которая при h=0 должна давать наилучшее приближение к истинному значению интеграла.

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

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

Легко показать, что исключение старшего члена погрешности формулы трапеций в точности даст формулу Симпсона:

Повторяя аналогичную процедуру для формулы Симпсона, получаем:

Если продолжить, вырисовывается такая таблица:

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

Элементы таблицы, как можно вывести из разложения (5), связаны рекуррентным соотношением:

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

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

Дополнительный метод добавляется в класс Quadrature

Проверим, как работает аппроксимация высокого порядка:

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

Некоторые замечания

Замечание 1. Количество вызовов функции в этих задачах характеризует число суммирований при вычислении интеграла. Уменьшение числа вычислений подынтегрального выражения не только экономит вычислительные ресурсы (хотя при более оптимизированной реализации и это тоже), но и уменьшает влияние погрешностей округления на результат. Так, при попытке вычислить интеграл тестовой функции метод трапеций зависает при попытке достигнуть относительной точности 5×10 -15 , метод парабол — при желаемой точности 2×10 -16 (что является пределом для чисел в двойной точности), а метод Ромберга справляется с вычислением тестового интеграла вплоть до машинной точности (с ошибкой в младшем бите). То есть, повышается не только точность интегрирования при заданном числе вызовов функции, но и предельно достижимая точность вычисления интеграла.

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

Замечание 3. Хотя метод Ромберга для ряда функций работает почти магическим образом, он предполагает наличие у подынтегральной функции ограниченных производных высоких порядков. Это значит, что для функций с изломами или разрывами он может оказаться хуже простых методов. Например, проинтегрируем f(x)=|x|:

Замечание 4. Может показаться, что чем выше порядок аппроксимации, тем лучше. На самом деле, лучше ограничить число столбцов таблицы Ромберга на уровне 4-6. Чтобы понять это, посмотрим на формулу (6). Второе слагаемое представляет собой разность двух последовательных элементов j-1-го столбца, поделенную на примерно 4 j . Т.к. в j-1-м столбце находятся аппроксимации интеграла порядка 2j, то сама разность имеет порядок (1/n_i) 2j

4 —ij . C учётом деления получается

7 второе слагаемое в (6) теряет точность после приведения порядков при сложении чисел с плавающей точкой, и повышение порядка аппроксимации может вести к накоплению ошибки округления.

Замечание 5. Желающие могут ради интереса применить описанные методы для нахождения интеграла и эквивалентного ему . Как говорится, почувствуйте разницу.

Заключение

Представлено описание и реализация базовых методов численного интегрирования функций на равномерной сетке. Продемонстрировано, как с помощью несложной модификации получить на базе метода трапеций класс квадратурных формул по методу Ромберга, что значительно ускоряет сходимость численного интегрирования. Метод хорошо работает для интегрирования «обычных» функций, т.е. слабо меняющихся на отрезке интегрирования, не имеющих особенностей на краях отрезка (см. Замечание 5), быстрых осцилляций и т.д.

Продвинутые методы численного интегрирования для более сложных случаев можно найти в книгах из списка литературы (в [3] — с примерами реализации на C++).

Содержание

Глава I. Общие сведения об интегральных уравнениях. 7

Глава II. Вычисление определенных интегралов на Mathcad. 11

2.1. Метод Ромберга. 11

2.2. Использование пакетов MathCAD для решения дифференциальных уравнений. 14

2.3. Метод Эйлера для дифференциальных уравнений первого порядка 16

2.4. Решение дифференциальных уравнений второго порядка. 17

Глава III Численные методы решения интегральных уравнений. 20

3.1. Квадратурный метод решения интегральных уравнений Фредгольма. 24

3.2. Квадратурный метод решения интегральных уравнений Вольтерры. 27

Глава IV. Прикладные задачи, использующие решение интегральных уравнений. 29

4.1. Расчет теплоизоляции. 29

4.2. Фильтр Калмана. 33

Листинг№1 Численное интегрирование. 39

1. Функция, возвращающая значение интеграла функции помощью метода Симпсона 39

2. Функция, возвращающая значение интеграла с помощью формулы трапеции 39

3. Функция, возвращающая значение интеграла, найденного по формулам треугольников. 40

Листинг № 2 Решение обыкновенных дифференциальных уравнений. 41

1.Функция, возвращающая численное решение ДУ методом Адамса. 41

2. Задание функции возвращающей решение ДУ методом Пикара. 42

3.Метод Эйлера. 42

Листинг №3. Решение линейного интегрального уравнения Вольтерра I-го рода. 43

Листинг №4. Решение линейного интегрального уравнения Вольтера II-го рода. 44

Листинг №6. Фильтр Калмана. 45

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

Введение.

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

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

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

cчитается задача Абеля, имеющая следующую формулировку:

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

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

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

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

В данной теме важным оказывается выбор базового программного средства.

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

Пакет Maple также весьма популярен. Кроме аналитических преобразований, пакет в состоянии решать задачи численно. Характерной особенностью пакета является то, что он позволяет конвертировать документы в формат LaTeX — стандартный формат подавляющего большинства научных издательств мирового класса. Кроме того, ряд других программных продуктов используют интегрированный символьный процессор Maple. Например, пакет подготовки научных публикаций Scientific WorkPlace позволяет обращаться к символьному процессору Maple, производить аналитические преобразования и встраивать полученные результаты в создаваемый документ[1].

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

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

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

Глава I. Общие сведения об интегральных уравнениях.

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

где — ядро интегрального преобразования, правая часть и являются заданными функциями, — параметр уравнения. Область интегрирования V может быть фиксированной (интегральные уравнения типа фредгольмовых) или переменной (интегральные уравнения типа вольтерровых).

Линейное интегральное уравнение первого рода получается при , и имеет вид:

Однородное линейное интегральное уравнение второго рода получается при и имеет вид:

Неоднородное интегральное уравнение второго рода получается при g(x) = 1 и имеет вид

Уравнения вида являются неоднородными.

Линейное интегральное уравнение Вольтерра первого рода имеет вид:

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

Представленная процедура решает уравнение методом квадратурных формул. Вычисление интеграла производится по формуле трапеций с постоянным шагом h:

,

где

Aj = 1 при j > 1 и Aj = 0.5 при j = 1

Линейное интегральное уравнение Вольтера второго рода имеет вид:

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

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

· Aj = 1 при j > 1 и Aj = 0.5 при j = 1.

Линейное интегральное неоднородное уравнение Фредгольма второго рода имеет вид:

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

· Aj = 1 при j, не равном 1 или n

· Aj = 0.5 при j, равном 1 или n.

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

Для сравнения с нолем в алгоритм передается малое число epsilon, и любое число, по модулю меньшее epsilon, считается нолем.

Глава II. Вычисление определенных интегралов на Mathcad

2.1. Метод Ромберга

Пусть требуется вычислить определенный интеграл на интервале [a;b].

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

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

Погрешность этой формулы равна .

Обозначим , где . Смысл введенного обозначения станет, ясен несколько позже.

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

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

В данном случае суммирование включает только один элемент. Обратите внимание, в новую оценку вошла старая оценка. Нам потребовалось определять значение функции только в новых узлах[3].

Если имеется 2n подинтервалов, то

Вообще, справедливо рекуррентное соотношение

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

Метод Ромберга заключается в том, что полученные оценки значения интеграла экстраполируют на случай бесконечного числа разбиений (величины шага равной нулю) по рекуррентной формуле

(1)

То есть строится следующий треугольник

R(5,1) R(5,2) R(5,3) R(5,4) R(5,5) ,

в котором первый столбец состоит из значений интеграла, полученных при последовательном удвоении числа интервалов. Второй столбец – результат уточнения значений первого столбца по рекуррентной формуле (1). Третий столбец – уточненные значения интеграла на основе второго столбца и т. д[4].

Формула (1) может быть получена различными способами. Можно, например, воспользоваться методом Невиля. Пусть имеется набор точек . Обозначим полином нулевой степени, проходящий через i-ю точку. Обозначим полином первой степени, проходящий через точки i и i+1. Совершенно аналогично будет означать полином n–1 степени, проходящий через все n точек. Легко убедиться, что

В нашем случае . В качестве выступают . Мы хотим получить значение интеграла в пределе , поэтому .

2.2. Использование пакетов MathCAD для решения дифференциальных уравнений.

Пусть необходимо найти решение уравнения

(2)

с начальным условием . Такая задача называется задачей Коши. Разложим искомую функцию в ряд вблизи точки и ограничимся первыми двумя членами разложения Учтя уравнение (2) и обозначив , получаем Эту формулу можно применять многократно, находя значения функции во все новых и новых точках.

(3)

Такой метод решения обыкновенных дифференциальных уравнений называется методом Эйлера. Геометрически метод Эйлера означает, что на каждом шаге мы аппроксимируем решение (интегральную кривую) отрезком касательной, проведенной к графику решения в начале интервала. Точность метода невелика и имеет порядок h. Говорят, что метод Эйлера – метод первого порядка, то есть его точность растет линейно с уменьшением шага h.

Существуют различные модификации метода Эйлера, позволяющие увеличить его точность. Все они основаны на том, что производную, вычисленную в начале интервала, заменяют на среднее значение производной на данном интервале[5]. Среднее значение производной можно получить (конечно же, только приближенно) различными способами. Можно, например, оценить значение производной в середине интервала и использовать его для аппроксимации решения на всем интервале

Можно также оценить среднее значение производной на интервале

Такие модификации метода Эйлера имеет уже точность второго порядка.

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

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

Еще один тип задач, часто встречающихся на практике, – краевые задачи. Пусть имеется дифференциальное уравнение второго порядка . Решение уравнения требуется найти на интервале , причем известно, что . Понятно, что произвольный интервал заменой переменных может быть сведен к единичному. Для решения краевой задачи обычно применяют метод стрельб. Пусть где k – некоторый параметр. Для некоторого пробного значения k может быть решена задача Коши, например, методом Рунге-Кутты. Полученное решение будет зависеть от значения параметра . Мы хотим найти такое значение параметра, чтобы выполнялось условие . Фактически мы свели исходную задачу к задаче решения трансцендентного уравнения с таблично заданной функцией. Если найдены такие значения параметра k1 и k2, что , то дальнейшее уточнение значения параметра можно проводить методом деления отрезка пополам[6].

2.3. Метод Эйлера для дифференциальных уравнений первого порядка

Решим задачу Коши для дифференциального уравнения первого порядка методом Эйлера.

Пусть правая часть уравнения равна

Зададим границы изменения x:

Зададим число точек и величину шага:

Зададим начальные условия:

Вычислим x и y по формулам Эйлера

Представим результат графически и сравним его с аналитическим решением

Точное аналитическое решение и решение, полученное численно, отличаются в точке x=1 на

То есть относительная ошибка составляет

2.4. Решение дифференциальных уравнений второго порядка

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

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

Пусть декремент затухания

Пусть циклическая частота

Зададим начальные условия

y0 соответствует начальной координате, а – начальной скорости. Зададим теперь матрицу D. С учетом того, что искомая величина соответствует нулевому элементу массива , ее первая производная – первому, а вторая – второму, имеем

Представим результаты расчета на графике и сравним их с аналитическим решением

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

Примечание: Mathcad имеет еще две функции для решения задачи Коши. Это функции Rkadapt и Bulstoer. Эти функции имеют те же самые аргументы и возвращают решения в такой же форме, что и функция rkfixed. Первая из этих функций использует метод Рунге–Кутты с переменным шагом, что позволяет повысить точность вычислений и сократить их объем, если искомое решение имеет области, где ее значения меняются быстро, и области плавного изменения. Функция Rkadapt будет варьировать величину шага в зависимости от скорости изменения решения[7].

Функция Bulstoer реализует иной численный метод – метод Булирша–Штёра. Ее следует применять, если известно, что решение является гладкой функцией.

Глава III Численные методы решения интегральных уравнений.

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

,

где D — некоторая область n-мерного пространства;

x — неизвестная функция, зависящая от времени;

K — функция относительно x(линейная или нелинейная).

Далее мы ограничим рассмотрение одномерным линейными интегральными уравнениями, в которой функция x(t) является функцией, зависящей от одной переменной, а область D – отрезком конечной длины, в каждой точке которого подъинтегральная функция K(t, s,x(s)) представима в виде Q(t, s)x(s).

Классификация типов линейных интегральных уравнений приводится по виду верхней границы интеграла в : если верхняя граница интегрирования является постоянной, то уравнение называется уравнением Фредгольма, если переменной — уравнением Вольтерры, которые, в свою очередь, подразделяются на уравнения первого и второго рода[8]. На практике наиболее широко применяются линейные интегральные уравнения второго рода:

,

,

где f(t) – неизвестная функция;

x(t)- решение уравнения;

Q(t, s)- ядро интегрального уравнения.

Ядро интегрального уравнения Фредгольма определяется на множестве точек квадрата [a, b]x[a, b],уравнения Вольтерры – в треугольнике .

Отметим, что доопределив ядро Q(t, s) уравнения Вольтерры нулем, в треугольнике уравнение Вольтера можно считать уравнением Фредгольма и применять для его решения методы уравнения Фредгольма. Однако при этом могут быть упущены некоторые специфические особенности уравнения Вольтерры, что определяет необходимость их раздельного рассмотрения[9].

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

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

,

,

относящиеся к классу некорректных задач.

Уравнения первого и второго рода можно записать в общем виде, используя функцию h(t), тождественно равную нулю для уравнений первого рода и единице — для уравнений второго рода:

.

Когда функция h(t) обращается в ноль в некоторых точках прямоугольника интегрирования, уравнение относится к интегральным уравнениям третьего рода. Уравнения данного типа встречаются в приложениях значительно реже, чем уравнения первых двух типов, значительно менее изучены.

Многие используемые на практике интегральные уравнения имеют ядро, зависящее только от разности t-s. Интегральные уравнения с данным типом ядра называются уравнениями с разностным ядром. Примером данного типа является уравнение, полученное в задаче Абеля[10].

Если Q(t, s) и f(t) – непрерывные функции, то при любых значениях параметра непрерывное решение уравнения Вольтерры второго рода существует, и оно единственное. Для уравнения Фредгольма второго рода при тех же требованиях единственное непрерывное решение существует, например, при условии, что

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

Известны формулы (или совокупность формул), позволяющие найти точное решение x(t). Например, решение уравнения Вольтерры, с с мультипликативным ядром

вычисляется по формуле :

Решение уравнения Фредгольма с вырожденным ядром

,

где числа — решения системы линейных алгебраических уравнений

(1)

;

.

Условие существования и единственности решения уравнения Фредгольма с вырожденным ядром, очевидно. Зависит от значения определителя D() системы линейных алгебраических уравнений (1), называемого определителем Фредгольма. Если D()≠0, то решение существует и единственно.

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

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

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

3.1. Квадратурный метод решения интегральных уравнений Фредгольма.

Заменим определенный интеграл

,

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

,

где j=1,2,…,n – номера узлов временной сетки; — весовые коэффициенты квадратурной формулы.

Подставив правую часть приближенного равенства с вместо интеграла в уравнение Фредгольма второго рода, получим

Данное выражение задает функцию, описывающую приближенное решение интегрального уравнения

,

Введем на отрезке [a, b] дискретную временную сетку узлы которой совпадают с узлами сетки . Для каждого момента времени выполняется равенство

,

И запишем равенство в виде системы n — линейных алгебраических уравнений с n неизвестными:

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

Таким образом, нахождение решения уравнения Фредгольма второго рода осуществляется в соответствии со следующим алгоритмом.

1. Задать временную сетку

2. Вычислить значение функции f(x) в узлах временной сетки.

3. Вычислить элементы матрицы, составленной из коэффициентов системы линейных алгебраических уравнений.

4. Решить систему линейных уравнений.

Точность численного решения интегрального уравнения зависит от нескольких факторов: применяемой квадратурной формулы, числа узлов временной сетки, свойств функции Q(t, s). В ряде книг приводятся аналитические выражения, позволяющие оценить максимальную погрешность численного решения при использовании различных вычислительных схем. Однако эти оценки оказываются малопригодными из-за их громоздкости, поэтому на практике используют менее строгий метод контроля точности численного решения — принцип Рунге.

Данный принцип состоит в сравнении численных решений, полученных на временных сетках с шагом 2h и h, в одних и тех же узлах временной сетки. Абсолютное значение разности данных решений характеризует величину погрешности численного решения. Недостаток данного подхода состоит в том, что при данном способе контроля приходится ограничиваться квадратурными формулами, пригодными только для сеток с равномерным шагом[11].

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

3.2. Квадратурный метод решения интегральных уравнений Вольтерры.

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

,где

Учитывая что уравнение Вольтерры формально можно считать уравнением Фредгольма вида:

,

K(t, s)=

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

Введем в рассмотрение временную сетку из [a, b], сотоящую из n узлов, и выберем конкретную квадратурную с весами , тогда приближенное решение интегрального уравнения принимает вид

,

Составим систему линейных алгебраических уравнений, аналогичную системе (1), которая в силу свойств ядра интегрального уравнения вырождается в треугольную:

Из данной системы видно, что искомые значения находятся последовательными вычислениями по следующим формулам:

где i=2,…,n.

Глава IV. Прикладные задачи, использующие решение интегральных уравнений.

4.1. Расчет теплоизоляции.

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

Определить наружный диаметр изоляции , при котором на внешней поверхности изоляции устанавливается температура . Определить: линейный коэффициент теплопередачи от воды к воздуху ; потери теплоты с одного метра длины трубопровода ; температуру наружной поверхности стального трубопровода . Наружный диаметр изоляции должен быть рассчитан с такой точность, чтобы температура отличалась от заданной не более чем на 0.1 K.

Упрощающие предположения: течение воды в трубе является термически стабилизированным; между сталью и асбестом существует идеальный тепловой контакт; теплопроводности стали и асбеста не зависят от температуры[12].

Вывод расчетных соотношений

Расчетная модель — бесконечная цилиндрическая труба, режим стационарный, объемных источников тепла в трубе нет. Уравнение теплопроводности имеет вид:

В цилиндрической системе координат получаем:

Граничные условия (закон Ньютона — Рихмана):

(положительным считается тепловой поток, идущий от центра )

Произведем замену переменных: ( ), (), тогда

()

()

где , ,.

(т. к. ). Тогда:

,

C другой стороны:

Приравнивая выражения для , находим постоянную K:

Тепловой поток, проходящий через стенку трубы на единице длины:

где

т. е. , где — линейный коэффициент теплопередачи

Окончательные расчетные формулы для определения температур (с учетом того, что ) имеют вид:

Для определения коэффициентов теплопередачи используются эмпирические формулы для переходного режима течения ():

, где .

Поскольку , а теплофизические свойства воды меняются с температурой не очень сильно, то и можно взять при температуре и считать, что . Тогда имеем:

Аналогично для естественной конвекции:

, где ,

Теплофизические свойства воды и воздуха берутся из книги [“Задачник по технической термодинамике и теории тепломассообмена” Под. ред. и Петражицкого , “Высшая школа” , 1986].

Задача решается методом последовательных приближений, первое приближение , последующие приближения находятся из соотношения:

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

4.2. Фильтр Калмана.

Представим себе некоторую систему, состояние которой в любой момент времени однозначно характеризуется определенным набором величин (например, координаты, скорости, уровни напряжения и т. д.), как правило, недоступных для непосредственного определения. Говоря терминами векторной алгебры, эти величины являются элементами вектора состояния системы, отнесенного к заданному моменту времени. Кроме того, имеется ряд переменных, некоторым образом связанных с состоянием системы, которые можно измерить с заданной точностью; такие величины составляют вектор измерений, относящихся к определенному моменту времени. Алгоритм фильтра Калмана позволяет в реальном времени построить оптимальную оценку состояния системы, основываясь на измерениях, неизбежно содержащих погрешности; при этом вектор измерений рассматривается в качестве многомерного выходного сигнала системы, отягощенного шумом, а вектор состояния — неизвестный многомерный сигнал, подлежащий определению. Условием оптимальности построенной оценки состояния является минимум ее средней квадратической ошибки[13].

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

Начальными условиями на каждом новом цикле алгоритма служат оценка состояния системы и величина, характеризующая ее погрешность. В случае скалярной переменной такой характеристикой является дисперсия, которая тем больше, чем сильнее разброс индивидуальных значений относительно истинного. Распространенная оценка дисперсии — среднеквадратическое отклонение, то есть квадрат стандартного отклонения, — выражает степень разброса величины относительно среднего. Обобщением дисперсии для вектора, то есть совокупности скалярных величин, служит ковариационная матрица. Ее диагональные элементы являются дисперсиями соответствующих составляющих вектора, а недиагональные — ковариациями, характеризующими взаимосвязь между парой составляющих. Совокупность измерений, отнесенных к каждому из моментов времени, обобщает вектор измерений. Алгоритм последовательно обрабатывает вновь поступающие векторы измерений, учитывая при этом значения, вычисленные на предшествующем цикле. Эта особенность отличает алгоритм фильтра Калмана от нерекуррентных алгоритмов, которым для работы требуется хранить весь массив обрабатываемых данных[14].

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

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

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

Если состояние рассматриваемой системы неизменно, то алгоритм фильтра Калмана сводится к последовательной форме классического метода наименьших квадратов, в котором матрица, обратная ковариационной, выступает в качестве весовой. Другими словами, фильтр Калмана является, по существу, рекуррентным способом решения задачи уравнивания по методу наименьших квадратов. Данная задача впервые решена в 1795 году, результаты были опубликованы в работе 1809 года под названием “Теория движения небесных тел”, в которой он применил метод наименьших квадратов к определению элементов орбит небесных тел (см. раздел “Замечания Гаусса”). Все изложенные в этой работе положения, касающиеся эффективности применения данного метода при обработке результатов измерений в равной степени относятся и к фильтру Калмана.

Применение фильтра Калмана в спутниковой навигационной аппаратуре.

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

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

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

Листинг№1 Численное интегрирование

1. Функция, возвращающая значение интеграла функции помощью метода Симпсона

2. Функция, возвращающая значение интеграла с помощью формулы трапеции.