|
|
| Строка 1: |
Строка 1: |
| | == Введение == | | == Введение == |
| - | === Постановка математической задачи ===
| |
| - | Рассмотрим задачу вычисления второй производной <tex>\frac{d^2}{d x^2}</tex> функции <tex>y(x)</tex>.
| |
| - | Если задачу точно решить либо не удаётся, либо слишком сложно, то приходится использовать методы численного дифференцирования.
| |
| - |
| |
| - | В данном отчёте теоретически обосновываются приближенные формулы для вычисления производной:
| |
| - | {{ eqno | 1 }}
| |
| - | ::<tex>y''(x)=(y(x+h)-2y(x)+y(x-h))/h^2+O(h^2)</tex>.
| |
| - |
| |
| - | {{ eqno | 2 }}
| |
| - | ::<tex>y''(x)=(-y(x+2h)+16y(x+h)-30y(x)+16y(x-h)-y(x-2h))/(12h^2)+O(h^4)</tex>
| |
| - |
| |
| - | Эти формулы можно использовать и для вычисления частной производной функции многих переменных.
| |
| - |
| |
| - | После этого приводятся оценки ошибок формул {{eqref|1}} и {{eqref|2}}, а также рассматривается вопрос о выборе оптимального шага <tex>h</tex>.
| |
| - |
| |
| | == Изложение метода == | | == Изложение метода == |
| - | === Определения ===
| |
| - | Пусть в некоторой окрестности точки <tex>x_0 \in \mathbb{R}</tex> определена функция <tex>y: U(x_0) \subset \mathbb{R} \to \mathbb{R}</tex>. Производной функции <tex>y(x)</tex> в точке <tex>x_0</tex> называется предел, если он существует,
| |
| - |
| |
| - | <tex>y'(x_0)=\lim_{x\to x_0} \frac{y(x)-y(x_0)}{x-x_0}</tex>.
| |
| - |
| |
| - | Второй производной функции <tex>y(x)</tex> в точке <tex>x_0</tex> называется производная функции <tex>y'(x)</tex> в точке <tex>x_0</tex>. Аналогично определяются производные и более высоких порядков.
| |
| - |
| |
| - | === Полиномиальные формулы ===
| |
| - | При численном дифференцировании функцию <tex>y(x)</tex> аппроксимируют легко вычисляемой функцией <tex>\varphi(x)</tex> и приближенно полагают <tex>y^{(k)}(x)\approx\varphi^{(k)}(x)</tex>. При этом можно использовать различные способы аппроксимации. Рассмотрим простейший случай - аппроксимацию интерполяционным многочленом Ньютона.
| |
| - |
| |
| - | Пусть задана сетка <tex>x_0<x_1<\dots<x_N. y(x)</tex> - исследуемая функция. Введем обозначение <tex>\xi_i=x-x_i</tex> и введем понятие разделенные разности <tex>y(x_0,x_1,\dots,x_k)</tex>, определяемые следующим образом:
| |
| - |
| |
| - | ::<tex>y(x_i,x_j)=\frac{y(x_i) - y(x_j)}{x_i - x_j}</tex>,
| |
| - |
| |
| - | ::<tex>y(x_i,x_j,x_k) = \frac{y(x_i,x_j)-y(x_j,x_k)}{x_i-x_k}</tex> и т.д.
| |
| - |
| |
| - | Можно доказать, что
| |
| - |
| |
| - | ::<tex>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}</tex>.
| |
| - |
| |
| - | Запишем интерполяционный многочлен Ньютона и продифференцируем его почленно:
| |
| - |
| |
| - | ::<tex>\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</tex>
| |
| - |
| |
| - | ::<tex>\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</tex>
| |
| - |
| |
| - | ::<tex>\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</tex>
| |
| - |
| |
| - | Общая формула примет следующий вид:
| |
| - | {{ eqno | 3 }}
| |
| - | ::<tex> \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] </tex>
| |
| - |
| |
| - | Обрывая ряд на некотором числе членов, получим приближенное выражение для соответствующей производной. Наиболее простые выражения получим, оставляя в формуле {{eqref|3}} только первый член:
| |
| - |
| |
| - | ::<tex>y'(x)\approx y(x_0,x_1) = \frac{y(x_0)-y(x_1)}{x_0-x_1}</tex>,
| |
| - |
| |
| - | {{ eqno | 4 }}
| |
| - | ::<tex>\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)</tex>,
| |
| - |
| |
| - | ::<tex>\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}</tex>
| |
| - |
| |
| - | Исследование точности полученных выражений при численных расчётах удобно делать при помощи апостериорной оценки, по скорости убывания членов ряда {{eqref|3}}. Если шаг сетки достаточно мал, то погрешность близка к первому отброшенному члену. Пусть мы используем узлы <tex>x_i, i=1\dots n</tex>. Тогда первый отброшенный член содержит разделенную разность <tex>y(x_0,x_1,\dots,x_{n+1})</tex>, которая согласно {{eqref|4}} примерно равна <tex>y^{(n+1)}(x)/(n+1)!</tex>. Перед ней стоит сумма произведений различных множителей <tex>\xi_i </tex>; каждое произведение содержит <tex>n+1-k</tex> множителей, а вся сумма состоит из <tex>C_{n+1}^k</tex> слагаемых. Отсюда следует оценка погрешности формулы {{eqref|3}} с <tex> n+1 </tex> узлами:
| |
| - |
| |
| - | {{eqno|5}}
| |
| - | ::<tex>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)}|}</tex>
| |
| - |
| |
| - | В частности, если сетка равномерная, то <tex>M_{n+1}=\max{|\xi_i|}\leq nh</tex>, откуда
| |
| - | {{eqno|6}}
| |
| - | ::<tex>R_n^{(k)}<M_{n+1} {\left( \frac{en}{n+1-k}h \right)}^{n+1-k}=O(h^{n+1-k}) </tex>.
| |
| - |
| |
| - | Стоит заметить, что строгое априорное исследование погрешности формулы {{eqref|3}}, аналогичное выводу остаточного члена многочлена Ньютона в форме Коши, для произвольного расположения узлов приводит к той же оценке {{eqref|5}}.
| |
| - |
| |
| - | Таким образом, порядок точности формулы {{eqref|3}} по отношению к шагу сетки равен числу оставленных в ней членов, или, что то же самое, он равен числу узлов интерполяции минус порядок производной. Поэтому минимальное число узлов, необходимое для вычисления <tex>k</tex>-й производной, равно <tex>k+1</tex>; оно приводит к формулам {{eqref|4}} и обеспечивает первый порядок точности. Эти выводы соответствуют общему принципу: при почленном дифференцировании ряда скорость его сходимости уменьшается.
| |
| - |
| |
| - | === Простейшие формулы (случай равномерной сетки)===
| |
| - |
| |
| - | Рассмотрим равномерную сетку <tex>x_0<x_1<\dots<x_N,\; x_{i+1}-x_i=h=const,\; y_i=y(x_i). </tex>
| |
| - | При этом вид формул {{eqref|3}} заметно упрощается, а точность нередко повышается.
| |
| - |
| |
| - | Рассмотрим сначала причину повышения точности. Остаточный член общей формулы {{eqref|3}} - это многочлен <tex>\sum\prod(x-x_i)</tex> степени <tex>n+1-k</tex> относительно <tex>x</tex>. Если <tex>x</tex> равен корню этого многочлена, то главный остаточный член обращается в нуль, то есть в этой точке формула имеет порядок точности на еденицу больше, чем согласно оценке {{eqref|6}}. Эти точки повышенной точности будем обозначать <tex>x_k^{(p)}</tex>, где <tex>k</tex> - порядок производной, а <tex>p=n+1-k</tex> - число оставленных в формуле {{eqref|3}} членов. Очевидно, <tex>p</tex>-членная формула имеет p точек повышенной точности.
| |
| - |
| |
| - | У одночленной формулы {{eqref|4}} для <tex>k</tex>-й производной точка повышенной точности на произвольной сетке определяется учловием <tex>\sum \xi_i=\sum(x-x_i)=0</tex>, что дает
| |
| - |
| |
| - | {{eqno|7}}
| |
| - | ::<tex>x_k^{(1)}=(x_0+x_1+\dots+x_k)/(k+1)</tex>;
| |
| - |
| |
| - | в этой точке одночленная формула имеет погрешность <tex>O(h^2)</tex> вместо обычной <tex>O(h)</tex>. Для двучленной формулы задача нахождения точек повышенной точности приводит к квадратному уравнению, корни которого действительны, но формула для их нахождения громоздка. Если <tex>p>2</tex>, то найти точки повышенной точности очень сложно, за исключением одного частного случая, который мы сейчас и рассмотрим.
| |
| - |
| |
| - | '''Утверждение.''' Пусть <tex>p</tex> нечетно, а узлы в формуле {{eqref|4}} выбраны так, что они расположены симметрично относительно точки <tex>x</tex>; тогда <tex>x</tex> является одной из точек повышенной точности <tex>x_k^{(p)}</tex>.
| |
| - |
| |
| - | На произвольной сетке условие симметрии реализуется только в исключительных случаях. Но если сетка равномерна, то каждый ее узел симметрично окружен соседними узлами. Это позволяется составить несложные формулы хорошей точности для вычисления производных в узлах сетки.
| |
| - |
| |
| - | Например, возьмем три соседний узла <tex>x_0</tex>, <tex>x_1</tex> и <tex>x_2</tex>. и вычислим первую и вторую производную в среднем узле. Выражая в одночленных формулах {{eqref|4}} разделенные разности через узловые значения функции, легко получим:
| |
| - |
| |
| - | ::<tex>y'(x_1)=(y_2-y_0)/2h+O(h_2),\; h=x_{i+1}-x_i=const,\; y_i=y(x_i)</tex>,
| |
| - |
| |
| - | {{eqno|8}}
| |
| - | ::<tex>y''(x_1)=(y_2-2y_1+y_0)/h^2+O(h^2)</tex>.
| |
| - |
| |
| - | Аналогично можно вывести формулы более высокого порядка или для более выслких производных. Например, формула для второй производной в центральной точке по пяти узлам:
| |
| - |
| |
| - | {{eqno|9}}
| |
| - | ::<tex>y''(x_2)=(-y_4+16y_3-30y_2+16y_1-y_0)/(12h^2)+O(h^4)</tex>
| |
| - |
| |
| - | Заметим, что формулы {{eqref|8}} и {{eqref|9}} написаны для случая равномерной сетки; применение их на неравномерной сетке приведет к грубой ошибке. Также формулы {{eqref|8}}и {{eqref|9}} совпадают с формулами {{eqref|1}} и {{eqref|2}}.
| |
| - |
| |
| - | На равномерной сетке для априоной оценки точности формул часто применяют способ разложения по формуле Тейлора-Макрорена. Предположим, например, что функция <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|8}} для второй производной имеем:
| |
| - |
| |
| - | ::<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|3}}.
| |
| - |
| |
| - | === Выбор оптимального шага (Регуляризация метода) ===
| |
| - | [[Изображение:Derives.PNG|thumb|200px|Рис.1]]
| |
| - | При численном дифференцированиии приходится вычитать друг из друга близкие значения функций. Это приводит к уничтожению первых значащих цифр, то есть к потере части достоверных знаков чила. Если значения функции известны с малой точностью, то встает вопрос - останется ли в ответе хотя один достоверный знак?
| |
| - |
| |
| - | Для ответа на этот вопрос исследуем ошибки при численном дифференцировании. При интерполировании обобщенным многочленом производная <tex>k</tex>-го порядка определяется согласно {{eqref|4}}, {{eqref|5}} формулой типа
| |
| - |
| |
| - | ::<tex>y^{(k)}(x)=h^{-k}\sum_q C_q(x)y(x_q) +R_k(x).\; C_q(x)=O(1)</tex>.
| |
| - |
| |
| - | Если формула имеет порядок точности <tex>p</tex>, то, значит, ее остаточный член равен <tex>R_k(x)\approx C(x)h^p</tex>. Этот остаточный член определяет погрешность метода, и он неограниченно убывает при <tex>h\to 0</tex>. Его зависимость от шага изображена на рис.1 жирной линией.
| |
| - |
| |
| - | Но есть еще неустарнимая погрешность, связанная с погрешностью функции <tex>\delta y(x)</tex>. Поскольку точный вид этой погрешности неизвестен, можно оценить только мажоранту неустранимой погрешности <tex>r_k=\delta*h^{-k}\sum_q|C_q|</tex>; она неограниченно возрастает при <tex>h \to 0</tex> (тонкая линия на рис.1). Фактически же неустранимая погрешность будет нерегулярно зависеть от величины шага, беспорядочно осцилируя в границах, определяемых мажорантой (точки на рис.1).
| |
| - |
| |
| - | Пока шаг достаточно велик, при его убывании неустранимая погрешность мала по сравнению с погрешностью метода; поэтому полная погрешность убывает. При дальнейшем уменьшении шага неустранимая погрешность становится заметной, что проявляется в не вполне регулярной зависимости результатов вычислений от величины шага. Наконец, при достаточно малом шаге неустранимая погрешность становится преобладающей, и при дальнейшем уменьшении шага результат вычислений становится все меньше достоверным.
| |
| - |
| |
| - | Полная погрешность мажорируется суммой <tex>R_k+r_k</tex> (штриховая кривая на рис.1). Оптимальным будет шаг, соответсвующий этой кривой. Нетрудно подсчитать, что
| |
| - |
| |
| - | {{eqno|10}}
| |
| - | ::<tex>h_0(\delta)={\left( \frac{k\delta \sum|C_q|}{pC} \right)}^{\frac{1}{p+k}}=O(\delta^{\frac{1}{p+k}})</tex>.
| |
| - |
| |
| - | ::<tex>\min(R_k+r_k)=Ch_0^p \left( 1+\frac{p}{k} \right) =O\left(\delta^{\frac{p}{p+k}} \right) </tex>.
| |
| - |
| |
| - | Меньший шаг невыгоден, а меньшая погрешность, вообще говоря, недостижима (хотя отдельные вычисления случайно окажутся более точными, но мы не сможем этого узнать). Эта минимальная ошибка тем менее, чем меньше погрешность входных данных и порядок вычисляемой производной и чем выше порядок точности формулы.
| |
| - |
| |
| - | Очевидно, при <tex>\delta y(x)\to 0</tex> можно получить сколь угодно высокую точность результата, если шаг стремится к нулю, будучи всегда не менее <tex>h_0(\delta)</tex>. Но если допустить <tex>h<h_0(\delta)</tex>, то результат предельного перехода может быть неправильным.
| |
| - |
| |
| - | Эта тонкость связана с некорректностью задачи дифференцирования. Рассмотрим погрешность входных данных вида <tex>\delta y(x)= m^{-1}\sin (m^2 x)</tex>. Она приводит к погрешности первой производной <tex>\delta y'(x) = m \cos (m^2 x)</tex>. При <tex>m \to \infty</tex> погрешность функции неограниченно убывает, а погрешность производной в той же норме неограниченно растет. Значит, нет непрерывной зависимости производной от функции, то есть дифференцирование некорректно. Особенно сильно это сказывается при нахождении производных высокого порядка.
| |
| - |
| |
| - | Изложенный выше способ определения оптимального шага и запрещение вести расчет шагом меньше оптимального есть некоторый способ регуляризации дифференцирования, так называемая регуляризация по шагу.
| |
| - |
| |
| - | Стоит заметить, что с ростом <tex>p</tex>, количества узлов в сетке для вычисления производных, константы <tex>C_q</tex> и <tex>C</tex> в формулах {{eqref|10}} повышаются, поэтому увеличение <tex>p</tex> разумно лишь в определенных пределах.
| |
| - |
| |
| | == Числовые примеры == | | == Числовые примеры == |
| - | Далее приведены графики зависимости ошибки вычисления от шага:
| |
| - | {| class="standard"
| |
| - | !
| |
| - | !<tex>\frac{d^2}{dx^2}e_x|_{x=0} </tex>
| |
| - | !<tex>\frac{d^2}{dx^2}\sin x|_{x=\frac{\pi}{2}} </tex>
| |
| - | !<tex>\frac{d^2}{dx^2} x^4|_{x=1} </tex>
| |
| - | |-
| |
| - | ! Метод {{eqref|1}}
| |
| - | |[[Изображение:Er_exp1.png|thumb|300px|Рис.1]]
| |
| - | |[[Изображение:Er_sin1.png|thumb|300px|Рис.1]]
| |
| - | |[[Изображение:Er_pol1.png|thumb|300px|Рис.1]]
| |
| - | |-
| |
| - | ! Метод {{eqref|2}}
| |
| - | |[[Изображение:Er_exp2.png|thumb|300px|Рис.1]]
| |
| - | |[[Изображение:Er_sin2.png|thumb|300px|Рис.1]]
| |
| - | |[[Изображение:Er_pol2.png|thumb|300px|Рис.1]]
| |
| - | |}
| |
| - | Вычисления проводились в стандартном типе [http://en.wikipedia.org/wiki/Double_precision double] (позволяет хранить 15 значащих десятичных цифр) языка [http://ru.wikipedia.org/wiki/C%2B%2B C++].
| |
| - |
| |
| - | Видно, что при слишком малом шаге полученные значения неадекватны.
| |
| - |
| |
| - | Оба метода дают хорошую точность, не превосходящую теоретические оценки.
| |
| - |
| |
| - | Сравнивая работу двух метотов вычисления производной, можно заметить, что второй из них получает ответ с большей точностью при меньшем шаге. Это соответсвует теоретическим выкладкам.
| |
| - |
| |
| - |
| |
| - |
| |
| | == Рекомендации программисту == | | == Рекомендации программисту == |
| - |
| |
| - | При программировании не стоит забывать, что вопрос выбора шага <tex>h</tex> очень важен.
| |
| - |
| |
| - | Исходя и экспериментов лучше всего брать шаг 0.0001 для метода {{eqref|1}} и 0.001 - 0.01 для метода {{eqref|2}}.
| |
| - |
| |
| | == Заключение == | | == Заключение == |
| - | В данном отчете получены теоретически и проверены практически два метода вычисления второй производной по одной переменной ({{eqref|1}} и {{eqref|2}} ).
| |
| - | И теоретические выкладки, и эксперименты показывают, что вопрос выбора шага - очень важный. Для его разрешения можно взять оптимальный шаг {{eqref|10}}, совпадающий с шагами показавшими наилучшие результаты при экспериментах.
| |
| - |
| |
| - | Экспериментально подтверждено, что оба метода дают неплохие результаты и могут быть использованы при практическом вычислении второй производной по одной переменной.
| |
| - |
| |
| | == Ссылки == | | == Ссылки == |
| | * [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]] | | * [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]] |
| - | * [[Media:DeriveExample1.zip| Код функций вычисления второй производной, ZIP [3Кб]]]
| |
| | | | |
| | == Список литературы == | | == Список литературы == |
| Строка 193: |
Строка 13: |
| | | | |
| | {{stub}} | | {{stub}} |
| - | [[Категория:Численное дифференцирование]]
| |
| | [[Категория:Учебные задачи]] | | [[Категория:Учебные задачи]] |
