Вычисление второй производной по одной переменной

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
Строка 57: Строка 57:
'''Замечание.''' Число узлов предполагалось произвольным; очевидно, симметричное расположение узлов относительно точки <tex>x_k^{(n)}</tex> означает, что при нечетном числе узлов точка <tex>x_k^{(p)}</tex> совпадает с центральным узлом, а при четном - лежит между средними узлами.
'''Замечание.''' Число узлов предполагалось произвольным; очевидно, симметричное расположение узлов относительно точки <tex>x_k^{(n)}</tex> означает, что при нечетном числе узлов точка <tex>x_k^{(p)}</tex> совпадает с центральным узлом, а при четном - лежит между средними узлами.
 +
 +
На произвольной сетке условие симметрии реализуется только в исключительных случаях. Но если сетка равномерна, то каждый ее узел симметрично окружен соседними узлами. Это позволяется составить несложные формулы хорошей точности для вычисления производных в узлах сетки.
 +
 +
Например, возьмем три соседний узла <tex>x_0</tex>, <tex>x_1</tex> и <tex>x_2</tex>. и вычислим первую и вторую производную в среднем узле. Выражая в одночленных формулах {{eqref|2}} разделенные разности через узловые значения функции, легко получим:
 +
 +
::<tex>y'(x_1)=(y_2-y_0)/2h+O(h_2)</tex>, <tex>h=x_{i+1}-x_i=const</tex>, <tex>y_i=y(x_i)</tex>,
 +
 +
{{eqno|6}}
 +
::<tex>y''(x_1)=(y_2-2y_1+y_0)/h^2+O(h^2)</tex>.
 +
 +
Аналогично можно вывести формулы более высокого порядка или для более выслких производных. Например, формула для второй производной в центральной точке по пяти узлам:
 +
 +
{{eqno|7}}
 +
::<tex>y''(x_2)=(-y_4+16y_3-30y_2+16y_1-y_0)/(12h^2)+O(h^4)</tex>
 +
 +
Заметим, что формулы {{eqref|6}} и {{eqref|7}} написаны для случая равномерной сетки; применение их на неравномерной сетке приведет к грубой ошибке.
 +
 +
На равномерной сетке для априоной оценки точности формул часто применяют способ разложения по формуле Тейлора-Макрорена. Предположим, например, что функция <tex>y(x)</tex> имеет непрерыную четвертую производную, и выразим значения функции в узлах <tex>x_{i\pm 1}</tex> через значиния функции и ее производных в центре симметрии узлов (в данном случае этим центром является узел <tex>x_i</tex>):
 +
 +
::<tex>y(x_{i\pm 1})=y(x_i \pm h)= y_i \pm h y_i'+\frac{1}{2}h^2 y_i'' \pm \frac{1}{6}h^3 y_i''' + \frac{1}{24}h^4 y^{(4)}(\eta_\pm)</tex>,
 +
 +
где <tex>\eta_{+}</tex> - некая точка интервала <tex>(x_i, x_{i+1})</tex>, а <tex>\eta_{-}</tex> - некоторая точка итервала <tex>(x_{i-1},x_i)</tex>. Подставляя эти разложения во вторую разность, стоящую в правой части формулы {{eqref|6}} для второй производной имеем:
 +
 +
::<tex>\frac{1}{h^2}(y_{i+1}-2y_i+y_{i-1})=y_i''+\frac{h^2}{24}[y^{(4)}(\eta_{+}+y^{(4)})(\eta_{-})]= y_i'' +O(h^2) </tex>.
 +
 +
Это подтверждает ранее сделанную оценку и уточняет величину остаточного члена, который оказался равен <tex>h^2y^{(4)}(\eta)/12</tex>. Такой способ получения остаточного члена проще, чем непосредственное вычисление по формуле {{eqref|1}}.
 +
 +

Версия 19:59, 15 октября 2008

Содержание

Введение

Постановка математической задачи

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

Изложение метода

Полиномиальные формулы

При численном дифференцировании функцию y(x) аппроксимируют легко вычисляемой функцией \varphi(x) и приближенно полагают y^{(k)}(x)\approx\varphi^{(k)}(x). При этом можно использовать различные способы аппроксимации. Рассмотрим простейший случай - аппроксимацию интерполяционным многочленом Ньютона. Вводя обозначение \xi_i=x-x_i, запишем это многочлен и продифференцируем его почленно:

\varphi(x)=y(x_0)+\xi_0 y(x_0,x_1)+\xi_0\xi_1 y(x_0,x_1,x_2) + \xi_0\xi_1\xi_2  y(x_0,x_1,x_2,x_3)+\dots
\varphi'(x)=y(x_0,x_1)+(\xi_0+\xi_1) y(x_0,x_1,x_2) + (\xi_0\xi_1+\xi_0\xi_2 +\xi_1\xi_2) y(x_0,x_1,x_2,x_3)+\dots
\varphi''(x)=2 y(x_0,x_1,x_2) + 2(\xi_0+\xi_1 +\xi_2) y(x_0,x_1,x_2,x_3)+\dots

Общая формула примет следующий вид:

( 1 )
 \varphi^{(k)}(x)=k!\left[ y(x_0,x_1,\dots,x_k) + \left( \sum_{i=0}^k \xi_i \right) y(x_0,x_1,\dots,x_{k+1}) + \left( \sum_{i>j\geq 0}^{i=k+1}\xi_i\xi_j\right)y(x_0,x_1,\dots,x_{k+2})   + \left( \sum_{i>j>l\geq 0}^{i=k+2}\xi_i\xi_j\xi_l\right)y(x_0,x_1,\dots,x_{k+3}) +\dots\right]

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

y'(x)\approx y(x_0,x_1) = \frac{y(x_0)-y(x_1)}{x_0-x_1},
( 2 )
\frac{1}{2}y''(x)\approx y(x_0,x_1,x_2) = \frac{1}{x_0-x_2}\left( \frac{y_0-y_1}{x_0-x_1}- \frac{y_1-y_2}{x_1-x_2}\right),
\frac{1}{k!} y^{(k)}(x) \approx y(x_0,x_1,\dots,x_k) = \sum_{p=0}^{k}y_p \prod_{i=0, i\neq p}^k {(x_p-x_i)}^{-1}

Исследование точности полученных выражений при численных расчетах удобно делать при помощи апостериорной оценки, по скорости убывания членов ряда (1). Если шаг сетки достаточно мал, то погрешность близка к первому отброшенному члену. Пусть мы используем узлы x_i, i=1\dots n. Тогда первый отброшенный член содержит разделенную разность y(x_0,x_1,\dots,x_{n+1}), которая согласно (2) примерно равна y^{(n+1)}(x)/(n+1)!. Перед ней стоит сумма произведений различных множителей \xi_i ; каждое произведение содержит n+1-k множителей, а вся сумма состоит из C_{n+1}^k слагаемых. Отсюда следует оценка погрешности формулы (1) с  n+1 узлами:

(3)
R_n^{(k)}\leq\frac{M_{n+1}}{(n+1-k)!}\max_i {|\xi_i|}^{n+1-k}, M_{n+1}=\max{|y^{(n+1)}|}

В частности, если сетка равномерная, то M_{n+1}=\max{|\xi_i|}\leq nh, откуда

(4)
R_n^{(k)}<M_{n+1} {\left( \frac{en}{n+1-k}h \right)}^{n+1-k}=O(h^{n+1-k}) .

Стоит заметить, что строгое априорное исследование погрешности формулы (1), аналогичное выводу остаточного члена многочлена Ньютона в форме Коши, для произвольного расположения узлов приводит к той же оценке (3).

Таким образом, порядок точности формулы (1) по отношению к шагу сетки равен числу оставленных в ней членов, или, что то же самое, он равен числу узлов интерполяции минус порядок производной. Поэтому минимальное число узлов, необходимое для вычисления k-й производной, равно k+1; оно приводит к формулам (2) и обеспечивает первый порядок точности. Эти выводы соответсвуют общему принципу: при почленном дифференцировании рдяа скорость его сходимости уменьшается.

Простейшие формулы

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

Рассмотрим сначала причину повышения точности. Остаточный член общей формулы (1) - это чногочлен \sum\prod(x-x_i) степени n+1-k относительно x. Если x равен корню этого многочлена, то главный остаточный член обращается в нуль, то есть в этой точке формула имеет порядок точности на еденицу больше, чем согласно оценке (4). Эти точки повышенной точности будем обозначать x_k^{(p)}, где k - порядок производной, а p=n+1-k - число оставленных в формуле (1) членов. Очевидно, p-членная формула имеет p точек повышенноу точности.

У одночленной формулы (2) для k-й производной точка повышенной точности на произвольной сетке определяется учловием \sum \xi_i=\sum(x-x
_i)=0, что дает

(5)
x_k^{(1)}=(x_0+x_1+\dots+x_k)/(k+1);

в этой точке одночленная формула имеет погрешность O(h^2) вместо обычной O(h). Для двучленной формулы задача нахождения точек повышенной точности приводит к квадратному уравнению, корни которого действительны, но формула для их нахождения громоздка. Если p>2, то найти точки повышенной точности очень сложно, за исключением одного частного случая, который мы сейчас и рассмотрим.

Утверждение. Пусть p нечетно, а узлы в формуле (2) выбраны так, что они расположены симметрично относительно точки x; тогда x является одной из точек повышенной точности x_k^{(p)}.

Доказательство. В самом деле, при этом величины \xi_i=x-x_i имеюи попарно равные абсолютные величины, но противоположные знаки. В остаточном члене множитель w=\sum\prod\xi_i имеет нечетную степень, и при одномерном изменении знаков всех \xi_iон должен изменить знак. Но поскольку одновременное изменение знаков \xi_i сводится при таком расположении узлов лишь к перемене их нумерации, то величина w должна сохраниться, что возможно только при w=0. Утверждение доказано.

Замечание. Число узлов предполагалось произвольным; очевидно, симметричное расположение узлов относительно точки x_k^{(n)} означает, что при нечетном числе узлов точка x_k^{(p)} совпадает с центральным узлом, а при четном - лежит между средними узлами.

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

Например, возьмем три соседний узла x_0, x_1 и x_2. и вычислим первую и вторую производную в среднем узле. Выражая в одночленных формулах (2) разделенные разности через узловые значения функции, легко получим:

y'(x_1)=(y_2-y_0)/2h+O(h_2), h=x_{i+1}-x_i=const, y_i=y(x_i),
(6)
y''(x_1)=(y_2-2y_1+y_0)/h^2+O(h^2).

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

(7)
y''(x_2)=(-y_4+16y_3-30y_2+16y_1-y_0)/(12h^2)+O(h^4)

Заметим, что формулы (6) и (7) написаны для случая равномерной сетки; применение их на неравномерной сетке приведет к грубой ошибке.

На равномерной сетке для априоной оценки точности формул часто применяют способ разложения по формуле Тейлора-Макрорена. Предположим, например, что функция y(x) имеет непрерыную четвертую производную, и выразим значения функции в узлах x_{i\pm 1} через значиния функции и ее производных в центре симметрии узлов (в данном случае этим центром является узел x_i):

y(x_{i\pm 1})=y(x_i \pm h)= y_i \pm h y_i'+\frac{1}{2}h^2 y_i'' \pm \frac{1}{6}h^3 y_i''' + \frac{1}{24}h^4 y^{(4)}(\eta_\pm),

где \eta_{+} - некая точка интервала (x_i, x_{i+1}), а \eta_{-} - некоторая точка итервала (x_{i-1},x_i). Подставляя эти разложения во вторую разность, стоящую в правой части формулы (6) для второй производной имеем:

\frac{1}{h^2}(y_{i+1}-2y_i+y_{i-1})=y_i''+\frac{h^2}{24}[y^{(4)}(\eta_{+}+y^{(4)})(\eta_{-})]= y_i'' +O(h^2) .

Это подтверждает ранее сделанную оценку и уточняет величину остаточного члена, который оказался равен h^2y^{(4)}(\eta)/12. Такой способ получения остаточного члена проще, чем непосредственное вычисление по формуле (1).



Числовой пример

Рекомендации программисту

Заключение

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

  • А.А.Самарский, А.В.Гулин.  Численные методы. Москва «Наука», 1989.
  • Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков.  Численные методы. Лаборатория Базовых Знаний, 2003.
  • Н.Н.Калиткин.  Численные методы. Москва «Наука», 1978.


Личные инструменты