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

Размещено на http://www.allbest.ru//

Размещено на http://www.allbest.ru//

МИНИСТЕРСТВО ОБРАЗОВАНИЯ НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

Брянский государственный университет имени академика И. Г. Петровского

Физико-математический факультет

Кафедра Информатики и прикладной математики

Курсовая работа

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

Выполнила:

Студентка 3 курса ФМФ

Камсулина М. И.

Научный руководитель:

кандидат физико-математических наук,

Трубников С.В.

Брянск 2016



Введение

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

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

Создание обзора литературы по численному интегрированию;

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

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



Глава 1. Методы численного интегрирования

Данная глава посвящена обзору методов численного интегрирования. Они находят обширное применение для решения многих математических задач [2]. В ней последовательно рассмотрены формулы Ньютона-Котеса, схемы Мотне-Карло, формулы Гаусса-Кристоффеля, кубатурные формулы.

1.1 Квадратурные формулы Ньютона-Котеса

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

Для написания параграфа использовалась книга [1].

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

. (1.1.1)

Один из наиболее широко используемых способов приближенного вычисления интегралов связан с применением квадратурных формул. Идея использования квадратурных формул заложена в самом определении интеграла: отрезок разбивается на n отрезков точками . На всех отрезках разбиения произвольным образом выбираются точки . Следует ввести обозначения: , . По определению интеграл Римана представляет собой предел интегральных сумм:

. (1.1.2)

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

, (1.1.3)

называемая квадратурной формулой. Точки называются узлами, а числа - коэффициентами квадратурной формулы. Разность вида

(1.1.4)

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

. (1.1.5)



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

. (1.1.6)

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

(1.1.7)

Если обозначить , то формула (1.1.7) превратится в квадратурную формулу. Квадратурные формулы, полученные таким способом, называют формулами Ньютона - Котеса.

Если записать формулу (1.1.7) при , , то ,



(1.1.8)

Формула (1.1.8) называется формулой левых прямоугольников. Геометрический смысл формулы левых прямоугольников: точное значение интеграла I равно площади криволинейной трапеции, ограниченной сверху графиком подынтегральной функции , справа и слева - вертикальными прямыми и , а снизу - осью x. А приближенное значение Q численно равно площади прямоугольника, ограниченного сверху горизонтальной прямой , справа и слева - вертикальными прямыми и , а снизу - осью x.

Если записать теперь формулу (1.1.7) при , , то

(1.1.9)

Формула (1.1.9) называется формулой правых прямоугольников. Геометрический смысл формулы правых прямоугольников: приближенное значение интеграла Q численно равно площади прямоугольника, ограниченного сверху горизонтальной прямой , справа и слева - вертикальными прямыми и , а снизу - осью x.

Если записать (1.1.7) при , , то ,

(1.1.10)



Формула (1.1.10) называется формулой средних прямоугольников. Геометрический смысл формулы средних прямоугольников: приближенное значение интеграла Q численно равно площади прямоугольника, ограниченного сверху горизонтальной прямой , справа и слева - вертикальными прямыми и , а снизу - осью x.

Если записать формулу (1.1.7) при , , , то

,

(1.1.11)

Формула (1.1.11) называется формулой трапеций. Геометрический смысл формулы трапеций: приближенное значение интеграла Q численно равно площади трапеции, ограниченной сверху прямой , справа и слева - вертикальными прямыми и , а снизу - осью x.

Если записать формулу (1.1.7) при , , , , то

,

. (1.1.12)

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

1.2 Формулы Гаусса

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

Для написания данного параграфа использовалась книга [1].

Пусть существует квадратурная формула

, (1.2.1)

которая является точной для любого многочлена, степень которого не превышает . Здесь - узлы, а - коэффициенты квадратурной формулы. Необходимо найти эти узлы и коэффициенты. Для этого нужно составить многочлен , степень которого равна . Здесь - многочлены Лежандра. При любом значении m от 0 до степень составленного многочлена не превышает . Поэтому для него должна быть точной квадратурная формула (1.2.1):

, (1.2.2)



причем эта формула справедлива при любом .

Разложив многочлен в ряд Фурье по системе многочленов Лежандра, получится следующая формула:

, (1.2.3)

где при . Кроме того, и при , поскольку разлагается многочлен n-ой степени. Поэтому разложение (1.2.3) превращается в равенство

. (1.2.4)

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

Пусть имеется многочлен -ой степени

(1.2.5)

Этот многочлен обладает очевидными свойствами:

(1.2.6)

Для многочлена (1.2.5) также должна быть точной квадратурная формула (1.2.1):

.

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

, . (1.2.7)

В итоге получилась квадратурная формула (1.2.1), узлы которой совпадают с корнями многочлена Лежандра , а коэффициенты вычисляются по формуле (1.2.7). Эта формула называется квадратурной формулой Гаусса для . Формулу (1.2.1) можно обобщить на случай произвольного отрезка интегрирования путем линейной замены переменной, отображающей на :

. (1.2.8)

Квадратурная формула Гаусса для произвольного отрезка примет вид

. (1.2.9)

Узлы общей формулы определяются из узлов с помощью преобразования (1.2.8):

, (1.2.10)

а коэффициенты получаются путем замены переменной (1.2.8) под знаком интеграла:

, ,

Где

В результате замены переменных получим:

, . (1.2.11)

Квадратурная формула (1.2.9) с узлами и коэффициентами, вычисляемыми по формулам (1.2.10) и (1.2.11), называется формулой Гаусса для произвольного отрезка интегрирования. При формула Гаусса совпадает с формулой средних прямоугольников.

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

Отрезок разбивается на n равных отрезков точками

, где , . (1.1.12)

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

. (1.1.13)

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

(1.1.14)

Полученная квадратурная формула (1.1.14) называется обобщенной формулой левых прямоугольников.

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

(1.1.15)

Полученная квадратурная формула (1.1.15) называется обобщенной формулой правых прямоугольников.

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

(1.1.16)

Полученная квадратурная формула (1.1.16) называется обобщенной формулой средних прямоугольников.

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

(1.1.17)

Полученная квадратурная формула (1.1.17) называется обобщенной формулой трапеций.

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

, (1.1.18)

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

(1.1.19)

Заменив точные значения интегралов, стоящих в сумме (1.1.19), на приближенные значения, получаемые по формуле Симпсона, получится формула:

(1.1.20)

Полученная квадратурная формула (1.1.20) называется обобщенной формулой Симпсона.

1.3 Метод Монте-Карло

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

Для написания данного параграфа использовалась книга [1].

Первая схема метода Монте-Карло.

Пусть - плотность распределения вероятности случайной величины . Случайная величина . Тогда математическое ожидание величины будет равно

,

а дисперсия

.

Если величина равномерно распределена на , то есть

то математическое ожидание

, (1.3.1)

а дисперсия

. (1.3.2)

Статистическими оценками (приближенными значениями) математического ожидания и дисперсии являются выборочное среднее

(1.3.3)

и выборочная дисперсия

(1.3.4)

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

Согласно центральной предельной теореме , по вероятности при . Кроме того, также представляет собой случайную величину, распределение вероятности которой стремится к нормальному (распределению Гаусса) при с дисперсией

(1.3.5)

Согласно правилу 3-х сигм, с вероятностью 0,997 будет выполняться неравенство

. (1.3.6)

Подставив сюда выражения (1.3.1), (1.3.3) ,(1.3.5) и умножив обе части неравенства на , получим:

(1.3.7)

Последнее неравенство выполняется при достаточно больших значениях n с вероятностью 0,997. Таким образом, если в качестве приближенного значения интеграла выбирается величина , то оценка погрешности этого приближенного значения с вероятностью 0,997 и при достаточно большом значении n будет равна

(1.2.8)

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

Вторая схема метода Монте-Карло.

Пусть на . Ищется приближенное значение Q интеграла . Точное значение интеграла равно площади криволинейной трапеции, ограниченной слева и справа прямыми , , а сверху и снизу - осью Ох и графиком функции . Пусть - площадь прямоугольника , содержащего указанную криволинейную трапецию, и - значения случайных величин и , равномерно распределенных на и соответственно. Эти значения можно получить с помощью датчика случайных чисел. Случайная точка на плоскости будет равномерно распределена в прямоугольнике . Поэтому вероятность попадания очередной случайной точки в криволинейную трапецию равна . Проведя n статистических испытаний, можно получить случайные значения и при . Теперь необходимо подсчитать количество попаданий точек в криволинейную трапецию. Условие попадания очередной точки в трапецию очевидно: . Тогда, согласно центральной предельной теореме, относительная частота попаданий будет стремиться к вероятности попадания при . Отсюда следует, что последовательность будет стремиться к точному значению интеграла I по вероятности при .

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

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

1.4 Кубатурные формулы

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

Для написания данного параграфа использовалась книга [1].

Требуется вычислить приближенное значение двойного интеграла

(1.4.1)

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

Тогда искомый двойной интеграл можно записать в виде повторного:

, (1.4.2)

Где

. (1.4.3)

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

(1.4.4)

где

, , (1.4.5)

Согласно (1.4.3),

, . (1.4.6)

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

(1.4.7)

где

,

, (1.4.8)

, .

Заменив в формуле (1.4.4) точные значения интегралов приближенными значениями , получится еще одно приближенное значение двойного интеграла:

(1.4.9)

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

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

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

(1.4.10)

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

, , , , …

Очевидно, , следовательно, .

Далее вычисляются последовательно приближенные значения при , и при каждом значении n проверяется выполнение неравенства

, . (1.4.11)

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

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



Глава 2. Программная реализация кубатурной формулы на основе формулы Симпсона

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

2.1 Кубатурная формула на основе формулы Симпсона

интеграл симпсон функция программный

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

Для написания данного параграфа использовалась книга [1].

Построение кубатурной формулы на основе формулы Симпсона будет осуществляться по аналогии с построением кубатурной формулы на основе формулы средних прямоугольников.

Пусть дана формула (1.4.1). Вначале необходимо задать сетку равноотстоящих узлов и шаг разбиения на оси OX. Для этого отрезок [a, b] разбивается на равных узлов точками

, (2.1.1)

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

(2.1.2)

Тогда, согласно (1.1.20), приближённое значение интеграла можно высчитать по формуле:

(2.1.3)

Теперь необходимо проделать то же самое для оси OY. Тогда

(2.1.4)

(2.1.5)

(2.1.6)

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

(2.1.7)

Подставив (2.1.7) в (2.1.3), получим искомую функцию.

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

(2.1.8)

где k - порядок точности формулы, r(0,1).



2.2Программная реализация кубатурной формулы на основе формулы Симпсона

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

В разделе описаний описаны 3 типа массивов: mas0 - одномерный массив с элементами вещественного типа, mas1 - одномерный массив с элементами целочисленного типа, mas2 - двумерный массив с элементами вещественного типа.

Далее описываются переменные, которые будут использованы в процессе реализации программы. Переменные a,b,hx,sum11,sum22,e,I11,I12,s1,s2,P являются переменными вещественного типа, i,j,m - целочисленного типа, sum1,sum2,x,hy,Q - типа mas0, n - типа mas1, y - типа mas2. Переменные a и b соответствуют значениям нижнего и верхнего пределов интегрирования, hx - отрезок разбиения (2.1.2), sum11 - значению суммы в (2.1.3), sum22 - значению суммы в (2.1.3), - погрешность приближённого значения (2.1.8), I11, I12 - формулы в (2.1.8), s2 и s1 - коэффициенты в (2.1.8) при , P - левая часть неравенства (2.1.8), i, j - индексы строк и столбцов в массивах, m - число в формуле (2.1.1), sum1 и sum2 - значения сумм и в формуле (2.1.7), х - точки разбиения отрезка [a, b] (2.1.1), hy - шаги разбиения по оси OY (2.1.6), Q - приближённое значение интеграла (2.1.7), n - числа, высчитываемые по формуле (2.1.5), у - точки разбиения, высчитываемые по формуле (2.1.4).

За разделом описаний следует раздел команд. Оператор writeln('Введите а'); выводит на экране сообщение «Введите а». Затем следует оператор readln(a); который считывает введённое пользователем значение и присваивает его переменной а. Те же самые операторы выполняется для переменной b. Далее следует строка hx:=(b-a)/m. В ней вычисляется шаг разбиения отрезка [a, b] (2.1.2).

Затем выполняется условный оператор, состоящий из 2m+1 итераций. В каждой итерации по формуле x[i]:=a+i*hx; генерируются значения x[i] элементов массива и при помощи оператора read(x[i]); эти значения заносятся в массив на i позицию. Данный оператор выполняет разбиение отрезка [a, b] на 2m+2 отрезка длиной hx (2.1.1).

Следующие условные операторы по своему принципу работы аналогичны. Формируются значения n[i] и считывает их на i позицию в массив n. Действие данного оператора аналогично выполнению (2.1.5).

После формируются значения шагов разбиения оси OY - hy[i], которые считываются в массив hy на i позицию. Работа данного оператора аналогична выполнению (2.1.6).

Далее выполняются два условных оператора, которые позволяют сформировать двумерный массив (2.1.4) в виде матрицы размера (2m+1)x(2n[2m]). Элементы y[i, j] считываются в i строку, j столбец матрицы у.

Затем элементам массивов sum1 и sum2, имеющим нулевой индекс, присваиваются нулевые значения. Это делается для того, чтобы впоследствии можно было воспользоваться рекуррентными формулами sum1[i]:=sum1[i-1]+x[i]*y[i,2*j]; и sum2[i]:=sum2[i-1]+x[i]*y[i,2*j-1];. Элементы массивов sum1 и sum2 представляют собой и соответственно.

После вычисления сумм, используя условный оператор, вычисляются приближённые значения интегралов по формуле Q[i]:=(hy[i]/3)/(x[i]*y[i,0]+x[i]*y[i,2*n[i]]+2*sum1[i]+4*sum2[i]); (2.1.7). Все вычисленные значения заносятся в массив Q на i позицию.

Далее, по формулам sum11:=sum11+Q[2*i]; и sum22:=sum22+Q[2*i-1]; вычисляются и соответственно из формулы (2.1.3).

Затем следует вычислить приближённое значение интеграла с заданной погрешностью. Для этого с помощью оператора writeln('Введите допустимую погрешность: '); на экран выводится указание о вводе допустимой погрешности. Оператор readln(e); считывает значение, указанное пользователем.

Следующий шаг - задание . Для этого переменным s1 и s2 задаются значения 1 и Ѕ соответственно, после чего переменной I11 присваивается выражение I11:=(hx*s1/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22); а I12 - выражение I12:=(hx*s2/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22);. Таким образом, I11 соответствует , а I12 - . После этого переменной Q присваивается выражение P:= abs(I12-I11)/7;, соответствующее .

Далее выполняется цикл с предусловием while P <= e do begin, в котором условием окончания итераций является выполнение неравенства Строка s1:=s2;s2:=s1/2; эквивалентна заданию последовательности , I11:=(hx*s1/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22); - заданию , I12:=(hx*s2/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22); - заданию , а P:=abs (I12-I11)/7; соответствует выражению (2.1.8).

В конце программы при помощи оператора writeln('Искомое значение интеграла равно ',I11,' с погрешностью, не превышающей ',e); пользователю выводится приближённое значение интеграла с погрешностью не превышающей e.

При использовании программа выдаёт приближённое значение интеграла, равное 0,167. Точное значение интеграла равно 1/6 (приложение 2).

Заключение

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

Был создан обзора литературы по численному интегрированию;

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

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



Литература

Акулич И. Л. Математическое программирование в примерах и задачах. Учебное пособие. - М.: Высшая школа, 1993.

Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. - Численные методы. - М., Наука, 1987.

Епанешников А. М., Епанешникова В. А. Программирование в среде Turbo Pascal 7.0. - M.: «ДИАЛОГ-МИФИ», 1996.

Краснов М. Л. Вся высшая математика: учебник / А. В. Краснов, А. И. Киселёв, Г. И. Макаренко, Е. В. Шишкин, В. И. Заляпин, С. К. Соболев; М: Эдиториал УРСС, 2001. - Т. 5. - 240с.

Соболь И. М. Численные методы Монте-Карло. М.: Наука, 1973.

Трубников С В. Численные методы. Часть 2: Аппроксимация функций. Численное дифференцирование и интегрирование: Учебное пособие для студентов вузов / С. В. Трубников - Брянск: Издательство БГУ, 2005. - 121с.

Трубников С. В. Компьютерное моделирование: Учебное пособие для студентов вузов. - Брянск: Изд-во БГУ, 2004.



Приложение 1

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

Program Sympson;

type mas0 = array [0..100] of real;

mas1 = array [0..100] of integer;

mas2 = array [0..100] of mas0;

var a,b,hx,sum11,sum22,e,I11,I12,s1,s2,P: real;

i,j,m: integer;

sum1,sum2,x,hy,Q: mas0;

n: mas1;

y: mas2;

begin

writeln('Введите а');

readln(a);

writeln('Введите b');

readln(b);

writeln('Введите число m: ');

readln(m);

hx:=(b-a)/2*m;

for i:=0 to 2*m do begin

x[i]:=a+i*hx;

read(x[i]);

end;

for i:=0 to 2*m do begin

n[i]:=trunc(((sqr(x[i]-1)+1)-(sqr(x[i]-1)-1))/2*hx)+1;

read(n[i]);

end;

for i:=0 to 2*m do begin

hy[i]:=((sqr(x[i]-1)+1)-(sqr(x[i]-1)-1))/2*n[i];

read(hy[i]);

end;

for i:=0 to 2*m do begin

for j:=0 to 2*n[i] do begin

y[i,j]:=sqr(x[i]-1)-1+j*hy[i];

read(y[i,j]);

end;

end;

sum1[0]:=0;

sum2[0]:=0;

for i:=0 to 2*m do begin

for j:=1 to n[i]-1 do begin

sum1[i]:=sum1[i-1]+x[i]*y[i,2*j];

read(sum1[i]);

end;

for j:=1 to n[i] do begin

sum2[i]:=sum2[i-1]+x[i]*y[i,2*j-1];

read(sum2[i]);

end;

Q[i]:=(hy[i]/3)/(x[i]*y[i,0]+x[i]*y[i,2*n[i]]+2*sum1[i]+4*sum2[i]);

read(Q[i]);

end;

sum11:=0;

for i:=1 to m-1 do begin

sum11:=sum11+Q[2*i];

end;

sum22:=0;

for i:=1 to m do begin

sum22:=sum22+Q[2*i-1];

end;

writeln('Введите допустимую погрешность: ');

readln(e);

s1:=1;s2:=1/2;

I11:=(hx*s1/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22);

I12:=(hx*s2/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22);

P:= abs(I12-I11)/7;

while P <= e do begin

s1:=s2;s2:=s1/2;

I11:=(hx*s1/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22);

I12:=(hx*s2/3)*(Q[1]+Q[2*(m+1)]+2*sum11+4*sum22);

P:=abs (I12-I11)/7;

end;

writeln('Искомое значение интеграла равно ',I11,' с погрешностью, не превышающей ',e);

end.



Приложение 2

Точное значение кратного интеграла при ,

Размещено на Allbest.ru

...