Дискуссионный математический форумМатематический форум

Математический форум Math Help Planet

Обсуждение и решение задач по математике, физике, химии, экономике

Теоретический раздел
Часовой пояс: UTC + 4 часа [ Летнее время ]
MathHelpPlanet.com RSS-лента Математического форума

Часовой пояс: UTC + 4 часа [ Летнее время ]


Методы численного дифференцирования и интегрирования

Методы численного дифференцирования и интегрирования


Общие понятия


Если функция [math]f(x)[/math] непрерывна на отрезке [math][a,b][/math] и известна ее первообразная [math]F(x)[/math], то определенный интеграл от этой функции может быть вычислен по формуле Ньютона-Лейбница


[math]\int\limits_{a}^{b} f(x)\,dx= \Bigl.{F(x)}\Bigr|_{a}^{b}= F(b)-F(a),[/math]

где [math]F'(x)=f(x)[/math]. Однако во многих случаях возникают большие трудности, связанные с нахождением первообразной, или эта задача не может быть решена элементарными способами. Например, в элементарных функциях не выражается интеграл [math]\textstyle{\int\limits_{1}^{2} \dfrac{dx}{\ln x}}[/math].


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


[math]\Omega_n= \bigl\{x_0,x_1,\ldots,x_n\bigr\},\quad x_{i+1}= x_i+h_{i+1},\quad i=\overline{0,n-1},\quad h_{i+1}= x_{i+1}-x_i.[/math]

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


Традиционно интегралы вычисляются с помощью квадратурных формул (явно), выражающихся линейными комбинациями значений функции [math]y_i=f(x_i)[/math] в узлах сетки [math]\Omega_n[/math], функционально-дифференциальной квадратурной формулы Эйлера-Маклорена, формул Гаусса-Кристоффеля, Маркова и различных нестандартных формул. В данной главе кратко описаны как классические, так и новые методы численного дифференцирования и интегрирования, разработанные на основе аппарата интегрально-дифференциальных сплайнов. Приведен ряд обобщенных явных квадратурных формул функционально-дифференциального типа, а также описан новый способ вычисления интегралов на основе решения систем линейных алгебраических уравнений. Последний подход в отличие от способа вычисления интегралов по квадратурным формулам носит неявный характер и связан с построением неявных алгоритмов.


Особенность аппарата численного дифференцирования и интегрирования, изложенного в данной главе, состоит в том, что в некоторых аппроксимацион-ных операторах дифференцирования могут использоваться как значения сеточной функции в узлах, так и значения определенных интегралов на частичных отрезках [math][x_i,x_{i+1}][/math], а в аппроксимационных операторах интегрирования — значения функции и ее производных в узлах. В соответствии с этим далее выделяются классические постановки задач численного дифференцирования и интегрирования, формулируемые обычно на равномерной сетке [math](h_{i+1}=\text{const})[/math], и обобщенные постановки, учитывающие связи производных с интегралами в задаче численного дифференцирования и интегралов с производными в задаче численного интегрирования, формулируемые на неравномерной сетке. Классические постановки ниже нумеруются цифрой 1, а обобщенные — цифрой 2. Данные постановки относятся только к явным аппроксимационным формулам численного дифференцирования и интегрирования. Соответствующие неявные алгоритмы рассматриваются отдельно.




Постановка задачи численного дифференцирования


Задача 1. Пусть на отрезке [math][a,b][/math] на равномерной сетке [math]\Omega_n~(h_{i+1}= h=\text{const})[/math] заданы:
а) сеточная функция [math]y_i= f(x_i),~ i=\overline{0,n}[/math] своими значениями [math]f_i=f(x_i)[/math];
б) точки [math]x_j[/math] сетки [math]\Omega_n[/math], в которых требуется найти значения производных;
в) желаемый порядок [math]t[/math] точности (аппроксимации) относительно [math]h[/math].

Требуется с заданным порядком точности (аппроксимации) вычислить значения производных [math]\Bigl.{\widehat{f}^{(p)}(x)}\Bigr|_{x=x_j}[/math] в точках [math]x_j[/math] сетки, где [math]p[/math] — порядок производной.


Иначе требуется получить аппроксимационный оператор [math]\widehat{f}^{(p)}(x_j)[/math], удовлетворяющий условию [math]\bigl|f^{(p)}(x_j)-\widehat{f}^{(p)}(x_j)\bigr|\leqslant C\,h^t[/math], где [math]C=\text{const}[/math], не зависящая от величины шага [math]h[/math].


Задача 2. Пусть на отрезке [math][a,b][/math] в общем случае на неравномерной сетке


[math]\Omega_n= \bigl\{x_0,x_1,\ldots,x_n\bigr\},\quad x_{i+1}=x_i+h_{i+1},\quad i=\overline{0,n-1},\quad h_{i+1}= x_{i+1}-x_i~ (h_{i+1}=\text{var})[/math]

заданы:


а) сеточная функция [math]y_i= f(x_i),~ i=\overline{0,n}[/math], своими значениями [math]f_i= f(x_i)[/math] и (или) значениями интегралов [math]\textstyle{I_{i}^{i+1}= \int\limits_{x_i}^{x_{i+1}} f(x)dx}[/math] на частичных отрезках [math][x_i,x_{i+1}][/math];


б) точки [math]x_j[/math] сетки [math]\Omega_n[/math], в которых требуется найти значения производных;


в) желаемый порядок [math]t[/math] точности (аппроксимации) относительно величины шага.


Требуется с заданным порядком точности (аппроксимации) вычислить значения производных [math]\Bigl.{\widehat{f}^{(p)}(x)}\Bigr|_{x=x_j}[/math] в точках [math]x_j[/math] сетки (получить аппроксимационный оператор), где [math]p[/math] — порядок производной.




Постановка задачи численного интегрирования


Задача 1. Пусть на отрезке [math][a,b][/math] на равномерной сетке [math]\Omega_n~ (h_{i+1}=h=\text{const})[/math] заданы:


а) сеточная функция [math]y_i= f(x_i),~ i=\overline{0,n}[/math], своими значениями [math]f_i= f(x_i)[/math] или сеточное представление формульной функции [math]y=f(x)[/math];


б) желаемый порядок [math]t[/math] точности (аппроксимации) относительно величины шага [math]h[/math].


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


[math]\widehat{I}_a^b\cong I_a^b= \int\limits_{a}^{b} f(x)\,dx\,.[/math]

Иначе требуется получить аппроксимационный оператор интегрирования [math]\widehat{I}_a^b[/math], удовлетворяющий условию [math]\bigl|\widehat{I}_a^b-I_a^b\bigr| \leqslant C\,h^t[/math], где [math]C=\text{const}[/math], не зависящая от [math]h[/math].


Задача 2. Пусть на отрезке [math][a,b][/math] в общем случае на неравномерной сетке


[math]\Omega_n= \bigl\{x_0,x_1,\ldots,x_n\bigr\}\quad (x_{i+1}= x_i+h_{i+1},~ i=\overline{0,n-1},~ h_{i+1}=x_{i+1}-x_i=\text{var})[/math]

заданы:


а) сеточная функция [math]y_i=f(x_i)~ i=\overline{0,n}[/math], своими значениями [math]f_i= f(x_i)[/math] и возможно значениями производных [math]f^{(p)}(x)[/math];


б) желаемый порядок [math]t[/math] точности (аппроксимации) относительно шага.


Требуется с заданным порядком точности вычислить значение интеграла (получить аппроксимационный оператор)


[math]\widehat{I}_a^b\cong I_a^b= \int\limits_{a}^{b} f(x)\,dx\,.[/math]

Отметим, что символом "^" здесь и далее обозначаются операторы дифференцирования и интегрирования.


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

1) разложения функции [math]f(x)[/math] по формуле Тейлора;

2) разложения первообразной функции [math]F(x)[/math] по формуле Тейлора;

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


Общая интегрально-функционально-дифференциальная формула для аппрок-симационного оператора дифференцирования [math]\widehat{h}^{(p)}(x_j)[/math] может быть записана в форме


[math]\widehat{h}^{(p)}(x_j)= \sum\limits_{k_1} \Bigl(\sum\limits_{l_1} a_{k_1,l_1}(h) \Bigr)I_{k_1}^{k_1+1}+ \sum\limits_{k_2} \Bigl(\sum\limits_{l_2} b_{k_2,l_2}(h)\Bigr)f_{k_2}+ \sum\limits_{k_3}\Bigl(\sum\limits_{l_3} c_{k_3,l_3}(h)\Bigr)f'_{k_3}+ \ldots,[/math]
(5.1)

где [math]a_{k_1,l_1}(h),\, b_{k_2,l_2}(h),\, c_{k_3,l_3}(h),\,\ldots[/math] — коэффициенты, индексы при [math]h[/math] для упрощения записи не указаны.


Общая функционально-дифференциальная формула для аппроксимационного оператора интегрирования [math]\widehat{I}_{i}^{i+1}[/math] может быть записана в виде


[math]\widehat{I}_{i}^{i+1}= h \sum\limits_{k_1}a_{k_1}(h) f_{k_1}+ h^2 \sum\limits_{k_2} b_{k_2}(h) f'_{k_2}+ h^3\sum\limits_{k_3}c_{k_3}(h) f''_{k_3}+\ldots,[/math]
(5.2)

где [math]a_{k_1}(h),\, b_{k_2}(h),\, c_{k_3}(h),\,\ldots[/math] - коэффициенты.


Возможность такого комбинированного представления следует из разложения первообразной [math]F(x)[/math] по формуле Тейлора относительно точки [math]x_i[/math] и последующего выражения из него производной некоторого порядка или определенного интеграла.


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


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


Аналогичная классификация справедлива и для формулы (5.2) численного интегрирования.




Метод Рунге уточнения результатов численного дифференцирования и интегрирования


Этот метод является апостериорным, и он основан на получении поправки, выражаемой через два приближенных результата [math]z_{hk}[/math] и [math]z_h[/math], получающихся на двух равномерных сетках с шагами [math]h[/math] и [math]kh[/math], где [math]k[/math] — коэффициент, характеризующий уменьшение шага [math](k<1)[/math] или его увеличение [math](k>1)[/math] при переходе с одной сетки на другую. Наиболее распространен способ генерации поправки в итерационном процессе с [math]k=1\!\!\not{\phantom{|}}\,2[/math] при необходимости снижения погрешности и [math]k=2[/math] при возможном ее увеличении.


Общая формула уточненного результата, получающегося на двух сетках с [math]h[/math] и [math]kh[/math], записывается в виде


[math]z= z_{hk}+\frac{z_{hk}-z_h}{(1\!\!\not{\phantom{|}}\,k)^p-1}\,,[/math]
(5.79)

где [math]p[/math] - порядок точности метода, по которому рассчитаны [math]z_{hk}[/math] и [math]z_h[/math]. Второе слагаемое в правой части (5.79) представляет собой поправку или апостериорную оценку точности результата, получаемого на сетке с параметром [math]kh[/math]. Если заданная точность не достигается, то при [math]k<1[/math] процесс продолжается путем дробления шага и перехода на новую, более мелкую сетку. Если же указанное.слагаемое намного меньше заданной точности, то в алгоритме необходимо перейти на укрупненную сетку, полагая, например, [math]k=2[/math].


Данную методику рассмотрим применительно к расчету производных на основе сгущающихся сеток [math](k=1\!\!\not{\phantom{|}}\,2)[/math], когда для первоначального вычисления производных принимается одна из формул (5.10). Схема уточнения приведена в табл. 5.5, где численные значения производных обозначены символом [math]Z[/math] с двумя нижними индексами, первый из которых (т>\) указывает на ступени уточнения ([math]m=1[/math] — первая, [math]m=2[/math] — вторая и т.д.), а второй — [math]k[/math] — на порядковые номера значений производных, полученных при различных шагах на каждой ступени уточнения. Начиная с третьей колонки таблицы помещаются значения [math]Z_{0,0},Z_{0,1},\ldots[/math], полученные по некоторой исходной аппроксимационной формуле при различных шагах разбиения сетки:


[math]h_0,\quad h_1=\frac{h_0}{2},\quad h_2=\frac{h_1}{2},\quad \ldots,\quad h_{k+1}=\frac{h_k}{2},\quad \ldots[/math]

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


[math]Z_{m,k}= Z_{m-1,k+1}+ \frac{Z_{m-1,k+1}-Z_{m-1,k}}{2^{2m}-1}\,,[/math]
(5.80)

где значение [math]2m[/math] соответствует порядку аппроксимации [math]p[/math].


Порядок уточнения численных результатов


Пусть по заданной сеточной функции [math]y_i= f(x_i),~ i=\overline{0,n}[/math], в некоторой фиксированной точке [math]x_j=x_i[/math] производная вначале рассчитывается по формуле (5.10) второго порядка аппроксимации


[math]Z_{0,k}= \frac{f(x_i+h_k)-f(x_i-h_k)}{2h_k}\quad (O(h^2)).[/math]
(5.81)

Тогда уточнение до порядка [math]O(h^6)[/math] производится следующим образом.


1. По формуле (5.81) рассчитываются значения [math]Z_{0,0},\, Z_{0,1},\, Z_{0,2}[/math] производных в заданной точке с использованием таблицы значений исходной сеточной функции. При этом для вычисления [math]Z_{0,k},~ k=0,1,2[/math], берутся значения [math]f(x_i)[/math] в тех точках сетки, которые соответствуют шагам [math]h_0,\, h_1[/math] и [math]h_2[/math]. Все величины, входящие в третью колонку таблицы, имеют порядок аппроксимации [math]O(h^2)[/math], что указано в нижней строке табл. 5.5.


2. Реализуется первая ступень процесса уточнения по формуле (5.80), в которой следует принять [math]m=1[/math] (знаменатель второго слагаемого формулы (5.80) равен 3), а [math]p=2m=2[/math]. При этом вычисляются [math]Z_{1,0}[/math] и [math]Z_{1,1}[/math]. Порядок аппроксимации результатов, получаемых в первой ступени уточнения, равен четырем.


3. Реализуется вторая ступень процесса уточнения по формуле (5.80), в которой следует принять [math]m=2[/math] (знаменатель второго слагаемого формулы (5.80) равен 15), а [math]p=2m=4[/math]. При этом вычисляется окончательное значение [math]Z_{2,0}[/math], порядок аппроксимации которого равен шести. Процесс уточнения до заданного порядка [math]p=6[/math] закончен.


Замечание. В случае, когда необходимо обеспечить порядок аппроксимации, превышающий полученный, например восьмой, производятся дополнительные расчеты значений [math]Z_{0,k}[/math] при более мелких шагах и расчеты [math]Z_{1,2},\, Z_{2,1}[/math] по первой и второй ступени уточнения. После этого рассчитывается [math]Z_{3,0}[/math] с порядком [math]O(h^8)[/math] (см. табл.5.5). Данный процесс может быть продолжен и далее.


Пример 5.13. Для сеточного представления функции [math]y=\sin x,~ x\in[0; 2\pi\!\!\not{\phantom{|}}\,5][/math], заданной на сетке с шагом [math]h= \pi\!\!\not{\phantom{|}}\,10[/math], на основе применения формулы (5.10) второго порядка и метода Рунге найти первую производную [math]\Bigl.{f'(x)}\Bigr|_{x=\frac{\pi}{5}}[/math] с порядком [math]O(h^4)[/math].


▼ Решение

1. По формуле (5.81) рассчитываются значения [math]Z_{0,0}[/math] и [math]Z_{0,1}[/math] производных в точке [math]x_2=\frac{\pi}{5}[/math]. При этом [math]Z_{0,0}[/math] получается при [math]h_0= \frac{2\pi}{10}[/math], а значение [math]Z_{0,1}[/math] при [math]h_1= \frac{\pi}{10}[/math]. Записывая формулу (5.81) для указанных шагов, получаем


[math]\begin{gathered}f'_{2,h}= Z_{0,0}= \frac{f_4-f_0}{2h_0}= \frac{0,\!951057-0}{2\cdot \pi\!\!\not{\phantom{|}}\,5}= 0,\!756827\,,\\[2pt] f'_{2,\frac{h}{2}}= Z_{0,1}= \frac{f_3-f_1}{2h_1}= \frac{0,\!809017-0,\!309017}{2\cdot \pi\!\!\not{\phantom{|}}\,10}= 0,\!795775\,. \end{gathered}[/math]

Исходная сеточная функция и результаты расчетов производной в точке [math]x_2[/math] помещены в табл. 5.6.


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


2. Реализуется первая (в данном случае и последняя) ступень процесса уточнения по формуле (5.80), в которой [math]m=1[/math]. При этом принимается [math]f'_{2,\frac{h}{2}}= Z_{0,1},~ f'_{2,h}= Z_{0,0}[/math]. Тогда


[math]Z_{1,0}\equiv f'_{2,\text{utochn}}= Z_{0,1}+ \frac{Z_{0,1}-Z_{0,0}}{2^2-1}= 0,\!795775+ \frac{0,\!795775-0,\!756827}{3}= 0,\!808758.[/math]

Полученный результат соответствует порядку аппроксимации [math]O(h^4)[/math].


Замечание. Для уточнения результатов численного интегрирования метод Рунге реализуется аналогично.


Часовой пояс: UTC + 4 часа [ Летнее время ]


Яндекс.Метрика

Copyright © 2010-2016 MathHelpPlanet.com. All rights reserved