Федоренко Р.П. Введение в вычислительную физику (Федоренко Р.П. Введение в вычислительную физику.djvu), страница 12

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

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

Он состоит в том, что с помощью специальных быстро работающих алгоритмов генерируется последовательность «случайных» точек единичного л-мерного куба х', хз, ..., хм е я", в каждой точке вычисляется значение У = У(х'"), и в качестве при- ближенного значения интеграла принимают величину М ' ~ У„,. Что касается алгоритмов, генерирующих случайные точки (соответствующие программы называют «датчиками случайных чисел»), то их разработка — отдельная достаточно тонкая наука. Эти точки должны быть «равномерно распределены» в л-мерном кубе, т,е.

они не должны «сбиваться в кучу» и оставлять в кубе «пустоты», в которые долго не попадают генерируемые точки. Можно еще иначе пояснить требование равномерности распределения случайных точек. Выделим в кубе некоторую часть О' не слишком причудливой формы и не слишком малой меры. Тогда число точек последовательности х', хз, ..., х"', попавших в о, должно быть близко к М шез о и не только асимптотически (при М- «ь), но и при конечных, не слишком больших М.

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

Добавим еще несколько слов о вычислении интеграла от У по слож- ВЫЧИСЛЕНИЕ ОПРЕДЕЛЕННЫХ ИНТЕГРАЛОЕ ной области й. Пусть область й есть часть единичного куба, выделенная неравенствами й)(х) ж О, Р' = 1, 2, ..., у. Техника вычисления интеграла по й выглядит примерно так. Генерируется та же последовательность псевдослучайных чисел, равномерно распределенных в единичном кубе.

Для каждой точки х проверяются неравенства а.(х ) жО, 1=1, 2, ..., у. Если все они выполнены, / т.е. х Е й, вычисляется,((хи), которая прибавляется к накапливающейся сумме. Кроме того, ведется подсчет числа попавших в й точек. Пусть вычислено М случайных точек, К из которых попало в й, и накоплена соответствующая сумма "', у(х"'«). Среднее значение у в области й можно определить величиной 1 ~ ~ х 1 И Вместе с тем это среднее вычисляется методом Монте-Карло: К 'х' У(х «). Что касается меры шез й, то она приближенно равна К/М.

Итак, имеем формулу к ) у(х) «1х= — Х. у(х"'«) + е . с« »=! О погрешности е» известно, что это случайная величина, математическое ожидание ее модуля имеет оценку СРЧ'К; константа С зависит от гладкости ~. Кващзатуры высокой точности, В качестве 7 можно взять интерполяцнонный полипом 2,„(1), построенный на сетке (1„)„'" . За счет того или иного выбора узлов сетки можно получить соответствующие преимущества. Широко используются чебышевские сетки, обеспечивающие устойчивость чебышевских квадратур при больших л«', Отметим еще гауссовы квадратуры.

Они основаны на следующей идее. При любой сетке квадратура, в которой используется интерполяционный полипом 1.„, точна в классе полиномов степени не выше ЛР, Можно так подобрать узлы сетки, что квадратура станет точной в классе поли номов степени не выше 2ЛР+ 1. Узлы сеток и коэффициенты квадратурных формул вычислены для разных ПР на стандартном интервале ( — 1, 1). Их можно найти в справочниках по методам приближенных вычислений (гауссовы узлы и коэффициенты, чебышевские и некоторые другие). 58 основы вьгчислитвльиой млтвмлтики й 5. Численное интегрирование задачи Коши дпя систем обыкновенных днффвренцивпьных уравнвний Рассмотрим задачу Коши для системы дифференциальных уравнений.

требуется найти функцию х(г), 0 ч г ч т, удовлетворяющую уравнению —,", =У(х, г), х(0) = х, Ок с<Т. (1) Здесь х — гримерный вектор, г(х, г) — заданная вектор-функция той же размерности, хо — заданная точка (данные Коши), Как известно, при весьма простых и общих предположениях о гладкости У решение задачи (1) существует н единственно. Хорошо также известно, что найти явное аналитическое выражение для х(г) удается в крайне редких случаях, н только для очень частных классов задач (например, линейных уравнений) имеются способы явного решения задачи Коши. Сама же задача (!) — одна нз наиболее часто встречающихся в различных приложениях (в физике, механике, астрономии, биологии, экономике и тд.). К таким задачам приходят, изучая движения планет и ракет, эволюцию биологических и экономических систем.

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

Начнем с основных обьектов метода сеток. На интервале (О, Т(, на котором ищется решение, введем покрывающую его дискретную сетку точек О=Г <Г,<Г «...Гн=т, (Гь) и 0' Ради простоты в дальнейшем будем в основном использовать равномерную сетку с шагом т = Т)И: Г„= нт. Приближенное решение будем искать в виде сенгочной функции, т.е. в виде функции дискретного аргумента и; обозначим ее как (х„)„'" .

Напомним, что каждый х„есть р-мерный вектор. С содержательной точки зрения х„будет представлять (приближенно) значение искомой функции х(г) в узле г„: х„х(г„). Сеточная функция (х„) не может удовлетворять никакому дифференциальному уравнению, н нужно построить какие-то другие уравнения, из которых можно было бы найти функцию (х„) и притом так, чтобы она была приближенным решением исходной задачи. й 51 ИИТЕГРИУОЕАИИЕ ЗАДАЧИ КОШИ ДЛЯ СИСТЕМ ОДУ Такие уравнения (так называемые разностные уравнения) строятся очень просто: входящие в дифференциальные уравнения производные заменяются соответствующими разностями. Это можно сделать разными способами, и они приводят к разным уравнениям. Например, "=/(х„, 1„) (явная схема Эйлера), (2) = /(х„+„ 1„Ч,) (неявная схема Эйлера), (3) х„+1- х„„1 = /(х„, 1„) (схема с центральной разностью). (4) Ограничимся пока этими простыми примерами, ниже мы рассмотрим и более сложные схемы.

Первая проблема, с которой мы сталкиваемся при решении разностных уравнений, — зто фактическое определение х„. Во всех случаях значения х„определяются последовательно, слева направо. Особенно просто вычисляются х„ при использовании явной схемы Эйлера: л = О, 1, ..., Ф вЂ” 1. х„ „ = х„ + т /(х„, 1„), Здесь х„известно (данные Коши). В правой части этой формулы используются уже найденные значения х„. В случае неявной схемы Эйлера ситуация несколько сложнее. Пусть значение х„уже найдено. Тогда хв м находится из уравнения х„„= х„+ т /(х„+,, 1„„), (5) которое являегся нелинейным относительно неизвестного х„+,, Правда, это слабая нелинейность, так как перед / стоит малый множитель т, что делает уравнение (5) не таким уж сложным.

Интуитивно ясно, что х„+, мало Отличается от х„, т.е. х„есть очень хорошее начальное приближение для какого-либо итерационного метода определения х„, (см. з 1). В схеме (4) мы сталкиваемся с характерным явлением: формальные порядки дифференциального и разносгного уравнений не совпадают. Под формальным порядком мы понимаем число произвольных постоянных в общем решении, или, если угодно, число дополнительных данных, полностью определяющих решение. В схемах Эйлера (2) и (3), как и в дифференциальном уравнении, достаточно задать х, чтобы все остальные значения х„хз, ... Определялись однозначно. В схеме (4) ситуация иная: только задав хв и х„мы опРеделим хз = хв + т /(х„1,) и т.д.

бо ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ Значение хо определяется постановкой задачи, х, формально можно задать каким угодно. Фактически, конечно, если мы стремимся получить хорошие результаты, х, нужно задавать достаточно аккуратно. Например, приемлемым (ио не лучшим) является х, = хо. Лучше (и правильнее) вычислять х, по какой-либо более простой разносгной схеме, например по явной схеме Эйлера: х, = хо + т У(х«). Обсудим некоторые вопросы, естественно возникающие при построении и использовании разностных схем.

Первый вопрос: чем следует руководствоваться при построении разностных уравнений? Интуитивно ясно, что разностные уравнения (2) — (4) имеют прямое отношение к исходному дифференциальному уравнению хотя бы потому, что в пределе (при т- 0) они переходят в (1). Но это соображение нужно четко оформить и дать ему количественное выражение. Таковым будет фундаментальное в вычислительной математике понятие аппроксимации. Второй вопрос. Пусть разносгное уравнение явно «соответствует» дифференциальному.