Презентация ИНТЕГРАЛ. Методы численного интегрирования в Mathcad
 Скачать 0.83 Mb.
|
|
Методы численного интегрирования в MathCAD Одномерный определённый интеграл Все рассмотренные ниже методы, в сущности, между собой похожи - если одномерный определённый интеграл есть площадь криволинейной трапеции под графиком: , то есть вопрос только в том, какой именно из простых зависимостей (прямая, парабола и т.п.) мы заменим подынтегральную функцию. Ясно, что можно заменить и вот так: Простейшее "интегрирование" - интеграл как площадь прямоугольника Площадь жирного прямоугольника приблизительно равна искомой площади под кривой, но это будет очень уж неточно, поэтому отрезок интегрирования по оси x всегда разбивают на небольшие интервалы (проще всего, с постоянным шагом h) и находят значение интеграла как сумму площадей простых фигур, например, прямоугольников, нижняя сторона которых равна h, а высота - значению f(x), взятому в некоторой точке интервала (на рисунке - в серединах): Определим тестовую функцию f(x), пределы интегрирования [a,b] и число интервалов n, на которое разбивается отрезок [a,b]. Величину шага h затем вычислим как (b-a)/n. Реализуем три основных метода прямоугольников. Разница между ними в том, в какой точке каждого отрезка на интервале интегрирования - левой, правой или в середине - берётся значение функции f(x). Численное интегрирование Пусть на отрезке [a,b] задана функция y=f(x). Разобьем отрезок на n элементарных отрезков с помощью точек x0, x1, . . . xn. x0=a; xn=b. На каждом отрезке выберем произвольную точку и найдем произведение значения функции в этой точке f(ζi) на длину отрезка Составим сумму этих произведений: Сумма Sn называется интегральной суммой. Определенным интегралом от функции f(x) на отрезке [a,b] называется предел интегральной суммы при неограниченном увеличении числа точек разбиения, при этом длина наибольшего из элементарных участков стремится к нулю xi-1< x i (1) (2) x0 x1 x2 xn x1 x2 a b Рис. 1. К вычислению интеграла S1 S2 Sn f(x) Выражение (1) при i=1,2,...,n описывают площади элементарных прямоугольников (рис.1). Интегральная сумма (2) - площадь ступенчатой фигуры, образуемой этими прямоугольниками. При неограниченном увеличении числа точек деления и стремлении к нулю всех элементов Dxi верхняя граница фигуры (ломаная) переходит в линию y=f(x). Площадь полученной фигуры - криволинейной трапеции - будет равна определенному интегралу. Во многих случаях, когда функция задана аналитически, определенный интеграл вычисляется с помощью неопределенного интеграла по формуле Ньютона - Лейбница. Она состоит в том, что определенный интеграл равен приращению первообразной F(x) на отрезке интегрирования: На практике часто нельзя воспользоваться формулой Ньютона – Лейбница по двум причинам: -вид функции f(x) не допускает непосредственного интегрирования, т.е. первообразную нельзя выразить в элементарных функциях; -значения функции f(x) заданы в виде таблицы. В этих случаях используются методы численного интегрирования. Они основаны на аппроксимации под интегральной функции некоторыми более простыми выражениями, например многочленами. Обычно используют локальную интерполяцию (Интерполя́ция, интерполи́рование — в вычислительной математике способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений). Формулы численного интегрирования называются квадратурами. Метод прямоугольников Этот метод использует замену определенного интеграла интегральной суммой. В качестве точек ζi могут выбираться левые ( ζi = xi-1 ) или правые ( ζi =xi ) границы участков. Для левых прямоугольников с постоянным шагом формула имеет вид: i=1,2, . . .n. (3) Для правых прямоугольников i=1,2, . . .n. (4) Более точной является формула прямоугольников, использующая значения функции в средних точках участков: i=1,2, . . .n. (5) и имеет порядок точности O(h3). Главный член погрешности формулы (5) на элементарном участке равен Рис. 5.2. Методы левых (а) и правых (б) прямоугольников. x x0 x1 xn f0 f1 x x0 x1 xn fn f1 б) а) Рис. 5.3. Метод средних прямоугольников. x x0 x1 xn y x y xi-1 xi fi-1 fi Рис. 5.4. Метод трапеций. hi f(x) i(x) x0 x1 x2 xn x1 x2 a b Рис. 1 S1 S2 Sn f(x) Интегрирование методами прямоугольников в MathCAD Метод трапеций Этот метод использует линейную интерполяцию, т.е. график функции y=f(x) представляется в виде ломаной, соединяющей точки (xi, yi). В этом случае площадь всей фигуры складывается из площадей элементарных трапеций (рис. 2). (xi-1,yi-1) (xi,yi) yi-1 yi xi-1 Рис. 2 Площадь каждой такой трапеции равна произведению полусуммы оснований на высоту: i=1,2, . . n. xi Рис. 7.3 Приведя в формуле (7.1) подобные члены, окончательно получим: (7.14) Формулу (7.14) называют формулой трапеций. Формулой трапеций часто пользуются для практических вычислений. Складывая все эти равенства, получим формулу трапеций: . Главный член погрешности этой формулы на элементарном участке равен .Точность интегрирования зависит от степени многочлена, количества участков и расположения точек. Во многих случаях формула центральных прямоугольников дает лучшую точность, чем формула трапеций. Программа вычисления интеграла методом трапеций Метод Симпсона Разобьем отрезок интегрирования [a,b] на четное число n равных частей с шагом h. На каждом отрезке подынтегральную функцию заменим многочленом второй степени: f(x) ji(x) = aix2 + bix +ci , xi-1< x < xi+1. Коэффициенты этих квадратных трехчленов могут быть найдены из условий равенства многочлена в точках xi соответствующим значениям функции yi. В качестве многочлена можно принять многочлен Лагранжа второй степени, проходящий через точки (xi-1,yi-1), (xi,yi), (xi+1,yi+1). Затем на каждом участке взять определенный интеграл и их просуммировать. В результате получим формулу: Главный член погрешности метода Симпсона на элементарном участке равен Программа вычисления интеграла методом Симпсона Если подынтегральная функция вычисляется по некоторому алгоритму, оформленному в виде функции с единственным параметром – аргументом функции, то можно предложить следующий программный блок, реализующий метод Симпсона: Параметры функции fsimpson: a, b – пределы интегрирования; n – число шагов интегрирования; f – имя функции, вычисляющей подынтегральное выражение. Интегрирование методами прямоугольников в MathCAD В методе трапеций мы для каждого отрезка интегрирования [xi,xi+1] соединяем отрезком прямой линии точки f(xi) и f(xi+1), считая интеграл как сумму площадей трапеций. Это всегда точнее, а сам метод ещё достаточно прост. По-моему, близок к оптимуму при массовых расчётах. В методе Симпсона (парабол) функцию f(x) на каждом отрезке интегрирования заменяют параболой, то есть, кривой второго порядка. Расчёт становится сложнее, но точность повышается в разы. Существует немало разновидностей формулы для метода Симпсона, вот 2 неплохих способа расчёта: Метод Симпсона (парабол) в MathCAD Увеличивая число интервалов n, можно оценить и порядок точности всех методов. Например, для метода первого порядка точности (методы левых и правых прямоугольников) при увеличении числа интервалов разбиения по оси x вдвое (n:=20 вместо n:=10 в начале документа) погрешность решения должна уменьшиться примерно в 2 раза. Для методов второго порядка точности (средних прямоугольников, трапеций) при уменьшении шага по x вдвое погрешность уменьшится примерно в 4 раза (второй по h порядок точности и означает, что погрешность уменьшается пропорционально величине h2). Метод Симпсона имеет четвёртый порядок точности, то есть, при уменьшении шага вдвое (увеличении вдвое числа интервалов n) погрешность решения уменьшится примерно в 24=16 раз. Решение обыкновенных дифференциальных уравнений Функция для решения дифференциальных уравнений имеет вид: Z: =rkfixed (y, x1, x2, npoints, D), где у – вектор начальных значений решений X1 – начало отрезка интегрирования X2 – конец отрезка интегрирования Npoints – число точек интегрирования D – вектор правых частей Число строк вектора D равно порядку уравнения. При этом дифференциальное уравнение записывается как бы в виде системы. На рис. 18 y1 - значение первой производной. Следовательно, первая строка вектора D – это первая производная, вторая строка – вторая производная. Результаты представляются в виде таблицы, где первый столбец число значений независимой переменной, которое равно npoints, второй – значения искомой функции, третий – значения первой производной. Для построения графика найденной функции необходимо по оси абсцисс отложить первый столбец таблицы, по оси ординат – второй. Столбец указывается верхним индексом (горячие клавиши Ctrl + 6) (см. рис. 18). Решение обыкновенного дифференциального уравнения. |