Участник:Александр Двойнев/песочница

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

(Различия между версиями)
Перейти к: навигация, поиск
(Заключение)
Строка 1: Строка 1:
== Введение ==
== Введение ==
-
Ставится задача вычислить интеграл вида
 
-
{{eqno |1}}
 
-
::<tex>J=\int_a^bf(t)dt,</tex>
 
-
где <tex>a</tex> и <tex>b</tex> - нижний и верхний пределы интегрирования; <tex>f(t)</tex> - непрерывная функция на отрезке <tex>[a,b]</tex>.
 
-
Пусть на отрезке интегрирования дана сетка, т.е. имеется совокупность узлов
+
Пусть задана функция <tex>f(x)</tex> действительного переменного. Требуется найти корни уравнения
-
<tex>\left\{{t_i}\right\}_{i = 0}^{N},\; i\in \left[{a, b}\right].</tex> Тогда интервал <tex>[a,b]</tex> разобьется на участки <tex>\tau_i = t_i-t_{i-1},\; i=1,\dots,N.</tex>.
+
{{eqno|1}}
-
Пусть также известны значения функции в узлах этой сетки, т.е. задана таблица
+
::<tex>f(x)=0.</tex>
-
<tex>f_i = \left\{{f(t_i)}\right\}_{i = 0}^{N}.</tex>
+
-
Представим интеграл {{eqref|1}} в виде суммы интегралов по частичным отрезкам:
+
-
{{eqno|2}}
+
-
::<tex>\int_a^bf(t)dt=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}f(t)dt.</tex>
+
-
Сущность большинства методов вычисления определённых интегралов состоит в замене подынтегральной функции <tex>f(t)</tex> на отрезке <tex>[t_{i-1},\;t_i]</tex> аппроксимирующей функцией <tex>\varphi(t)</tex>, для которой можно легко записать первообразную в элементарных функциях, т.е.
+
Задача нахождения корней уравнения {{eqref|1}} обычно решается в 2 этапа. На первом этапе проводится отделение корней, т.е. выделение отрезков, содержащих только один корень. На втором этапе, используя начальное приближение, строится итерационный процесс, позволяющий уточнить значение отыскиваемого корня.
-
{{eqno |3}}
+
-
::<tex>\int_{t_{i-1}}^{t_i}f(t)dt=\int_{t_{i-1}}^{t_i}\varphi(t)dt+R_i=S_i+R_i,</tex>
+
-
{{eqno|1*}}
+
-
::<tex>J=\sum_{i=1}^N \left( \int_{t_{i-1}}^{t_i}\varphi(t)dt+R_i \right)=S+R,</tex>
+
-
где <tex>S_i</tex> - приближённое значение интеграла на i-м отрезке, <tex>S=\sum_{i=1}^NS_i,</tex> а <tex>R_i</tex> - погрешность вычисления интеграла, <tex>R=\sum_{i=1}^NR_i,</tex>. Лучше всего изучена замена <tex>f(t)</tex> алгебраическим многочленом.
+
== Изложение метода ==
== Изложение метода ==
-
=== Общий случай ===
 
-
Возьмём в {{eqref|3}} в качестве аппроксимирующей функции [[Интерполяция кубическими сплайнами|кубический сплайн]]:
 
-
{{eqno |4}}
 
-
::<tex>\int_{t_{i-1}}^{t_i}f(t)dt = \int_{t_{i-1}}^{t_i}\varphi_i(t)dt+R_i,</tex> где
 
-
{{eqno|φ}}
+
===Метод касательных (Ньютона-Рафсона)===
-
::<tex>\varphi_i(t) = a_i+b_i(t-t_{i-1})+c_i(t - t_{i-1})^2 + d_i(t - t_{i-1})^3, \; t \in [t_{i-1},\;t_i]. </tex>
+
-
Известно, что коэффициенты <tex> a_i, b_i</tex> и <tex>d_i</tex> вычисляются по следующим формулам:
+
Пусть на отрезке <tex>[a,b]</tex> существует единственный корень уравнения {{eqref|1}}: <tex>x*</tex>
-
{{eqno|5a}}
+
{{eqno|2}}
-
::<tex> a_i=f_{i-1}</tex>
+
::<tex>f(x*)=0</tex>,
 +
а <tex>f'(x)</tex> существует, непрерывна и отлична от нуля на <tex>[a,b]</tex>. Перепишем {{eqref|2}} следующим образом:
 +
::<tex>f(x^k+(x*-x^k))=0</tex>
 +
и применим к этому выражению [[формула Лагранжа|формулу Лагранжа]]:
 +
::<tex>f(x^k)+f'(\bar{x})(x*-x^k)=0, \;\bar{x} \in [a,b].</tex>
 +
Заменим <tex> \bar x</tex> на <tex>x^k</tex>, а <tex>x*</tex> - на <tex>x^{k+1}</tex> и получим формулу итерационного процесса:
 +
::<tex>f(x^k)+f'(x^k)(x^{k+1}-x^k)=0.</tex>
 +
Выразим отсюда <tex>x^{k+1}</tex>:
 +
{{eqno|3}}
 +
::<tex>x^{k+1}=x^k-\frac{f(x^k)}{f'(x^k)}.</tex>
-
{{eqno|5b}}
+
Метод касательных имеет (когда сходится) квадратичную скорость сходимости: <tex>|x^{k+1}-x*|=O((x^k-x*)^2).</tex>
-
::<tex>b_i=\frac{f_i-f_{i-1}}{\tau_i}\;-\;\frac{\tau_i(2c_i+c_{i+1})}{3}</tex>
+
-
{{eqno|5d}}
+
===Модифицированный метод касательных===
-
::<tex>d_i=\frac{c_{i+1}-c_i}{3\tau_i}</tex>
+
-
А коэффициенты <tex>c_i</tex> являются решением [[СЛАУ]]:
+
Если мы хотим избежать вычисления производной на каждом шаге, то можно взять <tex> f'(x^0)</tex> вместо <tex>f'(x),</tex> где <tex>x^0</tex> - начальное приближение:
-
{{eqno|5c}}
+
::<tex>x^{k+1}=x^k-\frac{f(x^k)}{f'(x^0)}.</tex>
-
::<tex>\tau_{i-1}c_{i-1}+2(\tau_{i-1}+\tau_i)c_i+\tau_i c_{i+1}=3\left(\frac{f_i-f_{i-1}}{\tau_i} \; - \; \frac{f_{i-1}-f_{i-2}}{\tau_{i-1}}\right),\;\; 2 \le i \le N</tex>
+
В отличие от обычного метода касательных, в модифицированном методе предоставляется меньше требований к выбору начального приближения, а так же гарантировано отсутствие деления на ноль, если <tex>f'(x^0) \ne 0.</tex>
-
Нетрудно видеть, что матрица для решения [[СЛАУ]] {{eqref|5c}} есть трёхдиагональная матрица с диагональным преобладанием. Поэтому коэффициенты <tex>c_i</tex> можно вычислить с помощью метода прогонки.
+
-
Коэффициенты <tex>c_1</tex> и <tex>c_{N+1}</tex> получаются из условий свободных концов сплайна. Обычно требуют нулевую кривизну на концах сплайна и берут <tex>c_1=c_{N+1}=0.</tex>
+
Однако модифицированный метод имеет один существенный недостаток: линейную скорость сходимости: <tex>|x^{k+1}-x*|=O(x^k-x*).</tex>
-
В итоге интеграл {{eqref|1*}} запишется как сумма интегралов от сплайнов:
+
=== Геометрическая интерпретация===
-
::<tex>J=\sum_{i=1}^N\;\int_{t_{i-1}}^{t_i}\varphi_i(t)dt+R=\sum_{i=1}^N\left(a_i\tau_i+{b_i\over2}\tau_i^2+{c_i\over3}\tau_i^3+{d_i\over4}\tau_i^4\right)+R.</tex>
+
[[Изображение:kasat.jpg|thumb|200px|Метод касательных]]
 +
Метод Ньютона-Рафсона называют также методом касательных, т.к. новое приближение <tex>x^{k+1}</tex> является абсциссой точки пересечения касательной, проведенной в точке <tex>(x^k,f(x^k))</tex> к графику функции <tex>f(x)</tex>, с осью ОХ.
-
Последняя формула упрощается при подстановке в неё выражений {{eqref|5a}}, {{eqref|5b}} и {{eqref|5d}} для коэффициентов <tex>a_i,\;b_i</tex> и <tex>d_i:</tex>
+
== Сходимость метода ==
-
{{eqno|6}}
+
Заметим, что метод касательных является частным случаем [[Метод простых итераций|метода простых итераций]]
-
::<tex>J = \sum_{i=1}^N\;\frac{f_i+f_{i-1}}{2}\tau_i\;-\;\sum_{i=1}^N\;\frac{\tau_i^3(c_{i+1}+c_i)}{12}+R.</tex>
+
::<tex>x^{k+1}=g(x^k),\; k=0,1,\ldots,</tex>
 +
для которого
 +
{{eqno|4}}
 +
::<tex>g(x)=x-\frac{f(x)}{f'(x)}.</tex>
-
=== Случай с равномерной сеткой ===
+
[[Метод простых итераций]] сходится тогда и только тогда, когда
 +
::<tex>|g'(x)|\le q < 1,</tex>
-
Пусть на отрезке <tex>[a,b]</tex> задана равномерная сетка, т.е. <tex>\tau_1=\tau_2=\dots=\tau_N=\tau.</tex> Тогда выражение {{eqref|6}} перепишется в виде:
+
Подставим в последнее условие выражение для g(x) {{eqref|4}} и получим условие сходимости метода касательных:
 +
::<tex>\left|\frac{f(x)f''(x)}{(f'(x))^2}\right| \le q < 1.</tex>
-
{{eqno|7}}
+
Сформулируем теорему о квадратичной скорости сходимости метода касательных:
-
::<tex>J\;=\;\frac{\tau}{2}(f_0+f_N)+\tau\sum_{i=1}^{N-1}f_i-\frac{\tau^3}{6}\sum_{i=2}^Nc_i+R.</tex>
+
-
Просуммируем уравнения {{eqref|5c}} от i=2 до N. Получим:
+
'''Теорема.''' Пусть <tex>x*</tex> - простой вещественный корень уравнения <tex>f(x) = 0</tex>, а функция <tex>f(x)</tex> - дважды дифференцируема в некоторой окрестности <tex>U_r(x*)</tex>, причем первая произодная нигде не обращается в нуль.
-
{{eqno|8}}
+
Тогда, следуя обозначениям
-
::<tex> 6\sum_{i=2}^Nc_i=(c_2+c_N)+\frac{3}{\tau^2}(f_0-f_1+f_N-f_{N-1}).</tex>
+
::<tex>0 < m_1 = \;\inf_{x\in U_r(x*)}|f'(x)|,\;\; M_2 = \;\sup_{x\in U_r(x*)}|f''(x)|</tex>,
-
 
+
при выборе начального приближения <tex>x^0</tex> из той же окрестности <tex>U_r(x^*)</tex> такого, что
-
Подставим {{eqref|8}} в {{eqref|7}} и получим окончательное выражение для <tex>J</tex>:
+
::<tex>\frac{M_2|x^0 - x*|}{2m_1} = q < 1</tex>,
-
 
+
итерационная последовательность
-
{{eqno|9}}
+
::<tex>x^{k+1} = x^k - \frac{f(x^k)}{f'(x^k)},\; k = 0,1,\ldots</tex>
-
::<tex> J\;=\;\tau\sum_{i=2}^{N-2}f_i+\frac{\tau}{12}(5f_0+13f_1+13f_{N-1}+5f_N)-\frac{\tau^3}{36}(c_2+c_N)+R</tex>
+
будет сходиться к <tex>x*</tex>, причем для погрешности на k-м шаге буддет справедлива оценка:
-
 
+
::<tex>|x^k - x*| \le q^{2^k - 1}|x^0 - x*|.</tex>
-
Несмотря на то, что <tex>c_2</tex> и <tex>c_N</tex> все равно придется вычислять методом прогонки, точность и скорость вычисления приближенного значения интеграла будут увеличены за счет сокращения числа слагаемых.
+
-
 
+
-
== Анализ метода и ошибок ==
+
-
 
+
-
Анализ формулы {{eqref|6}} показывает, что первый член в правой части совпадает с [[Методы прямоугольников и трапеций|формулой трапеций]]. Следовательно, второй член характеризует поправку к методу трапеций, которую дает использование сплайнов.
+
-
 
+
-
Как следует из формулы {{eqref|φ}}, коэффициенты <tex>c_i</tex> выражаются через вторые производные <tex>\varphi_i''(t):</tex>
+
-
 
+
-
::<tex>c_i=\frac{1}{2}\varphi_i''(t_{i-1})\approx\frac{1}{2}f_{i-1}''.</tex>
+
-
 
+
-
Это позволяет оценить второй член правой части формулы {{eqref|6}}:
+
-
 
+
-
::<tex>\frac{\tau_i^3}{12}(c_{i-1}+c_i)\approx\frac{\tau_i^3}{12}f_*'',</tex>
+
-
 
+
-
где <tex>f_*''</tex> - вторая производная в некоторой внутренней точке. Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций.
+
== Числовой пример ==
== Числовой пример ==
-
[[Изображение:Splineint.jpg|thumb|200px| рис.1. График функции]]
+
Рассмотрим функцию <tex>f(x)=
-
Рассмотрим функцию <tex>f(t)=sin(t).</tex> Вычислим с помощью сплайн-квадратур приближенное значение интеграла
+
-
 
+
-
::<tex>\int_0^ \pi f(t)dt=\int_0^\pi sin(t) dt=2</tex>
+
-
 
+
-
и исследуем поведение погрешности. Результаты работы [[Media:Splineint.zip|программы]] проиллюстрируем в виде графика зависимости логарифма ошибки от числа разбиений:
+
-
 
+
-
{| class="wikitable"
+
-
|[[Изображение:Splineint100.jpg|thumb|400px| рис.2. Число разбиений: 100]]
+
-
|[[Изображение:Splineint1kk.jpg|thumb|400px| рис.3. Число разбиений: 1 000 000]]
+
-
|}
+
-
 
+
-
При этом красным цветом показан результат рассчетов по формуле {{eqref|9}}, а зеленым - по формуле {{eqref|6}}.
+
-
 
+
-
Как видно из рис.2, в данном примере бессмысленно брать число разбиений большее, чем 30. Более того, как видно из рис.3, при большом числе разбиений погрешность начинает расти.
+
== Рекомендации программисту ==
== Рекомендации программисту ==
-
 
-
=== Пример программы ===
 
-
 
-
Ниже приведен пример программы на языке C++, считающей приближенное значение интеграла с помощью сплайн-квадратур:
 
-
[[Media:Splineint.zip|Splineint.zip [0.7Кб]]]
 
-
 
-
'''Некоторые комментарии по работе с программой:'''
 
-
 
-
В 6-й строке <code>const int N=1000;</code> N - число отрезков <tex>\tau_i.</tex>
 
-
 
-
В 8-й строке <code>const double a=0,b=pi;</code> <tex>a</tex> и <tex>b</tex> - пределы интегрирования.
 
-
 
-
В 18-й строке <code>f[i]=sin(x);</code> f[i] - интегрируемая функция.
 
-
 
-
===Ошибки округления===
 
-
 
-
Следует помнить, что при реализации метода помимо теоретической априорной погрешности из-за суммы большого числа чисел возникает также погрешность округления.
 
-
 
== Заключение ==
== Заключение ==
-
 
-
Итак, мы получили, что погрешность сплайн-квадратуры меньше, чем погрешность метода трапеций. Однако [[алгоритм]] интегрирования с помощью сплайнов сложнее алгоритмов методов трапеций и Симпсона за счет необходимости решения СЛАУ для опрееления коэффициентов сплайнов <tex>c_i.</tex> Также при решении СЛАУ теряется точность. Поэтому рационально использовать сплайн-квадратуры в комплексе, когда сплайны применяются для сглаживания зависимостей, обработки эксперимениальных данных и т.п.
 
-
 
-
Также одним из основных преимуществ метода является возможность работать с произвольным разбиением орезка.
 
-
 
== Ссылки ==
== Ссылки ==
* [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
* [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
-
* [[Интерполяция кубическими сплайнами]]
+
* [[Метод секущих]]
-
 
+
== Список литературы ==
== Список литературы ==
-
* http://www.intuit.ru/department/calculate/calcmathbase/7/1.html
+
*[http://mmphome.1gb.ru/index.php?pid=show&id=79 Н.В.Соснин. Численные методы. Конспект лекций (сост. Д.В.Ховратович, Е.А.Попов)]
-
* http://mathalgo.blogspot.com/2007/11/blog-post.html
+
*Самаский А.А., Гулин А.В. Численные Методы. Учеб. пособие для вузов. - М.:Наука, 1989.
-
* http://myhomepage.h17.ru/Lect06/lect06.htm#02
+
 
-
* А.Е. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск:МП "РАСКО", 1991.
+
<!--{{Stub|}}-->

Версия 11:57, 4 декабря 2008

Содержание

Введение

Пусть задана функция f(x) действительного переменного. Требуется найти корни уравнения

(1)
f(x)=0.

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

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

Метод касательных (Ньютона-Рафсона)

Пусть на отрезке [a,b] существует единственный корень уравнения (1): x*

(2)
f(x*)=0,

а f'(x) существует, непрерывна и отлична от нуля на [a,b]. Перепишем (2) следующим образом:

f(x^k+(x*-x^k))=0

и применим к этому выражению формулу Лагранжа:

f(x^k)+f'(\bar{x})(x*-x^k)=0, \;\bar{x} \in [a,b].

Заменим  \bar x на x^k, а x* - на x^{k+1} и получим формулу итерационного процесса:

f(x^k)+f'(x^k)(x^{k+1}-x^k)=0.

Выразим отсюда x^{k+1}:

(3)
x^{k+1}=x^k-\frac{f(x^k)}{f'(x^k)}.

Метод касательных имеет (когда сходится) квадратичную скорость сходимости: |x^{k+1}-x*|=O((x^k-x*)^2).

Модифицированный метод касательных

Если мы хотим избежать вычисления производной на каждом шаге, то можно взять  f'(x^0) вместо f'(x), где x^0 - начальное приближение:

x^{k+1}=x^k-\frac{f(x^k)}{f'(x^0)}.

В отличие от обычного метода касательных, в модифицированном методе предоставляется меньше требований к выбору начального приближения, а так же гарантировано отсутствие деления на ноль, если f'(x^0) \ne 0.

Однако модифицированный метод имеет один существенный недостаток: линейную скорость сходимости: |x^{k+1}-x*|=O(x^k-x*).

Геометрическая интерпретация

Метод касательных
Метод касательных

Метод Ньютона-Рафсона называют также методом касательных, т.к. новое приближение x^{k+1} является абсциссой точки пересечения касательной, проведенной в точке (x^k,f(x^k)) к графику функции f(x), с осью ОХ.

Сходимость метода

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

x^{k+1}=g(x^k),\; k=0,1,\ldots,

для которого

(4)
g(x)=x-\frac{f(x)}{f'(x)}.

Метод простых итераций сходится тогда и только тогда, когда

|g'(x)|\le q < 1,

Подставим в последнее условие выражение для g(x) (4) и получим условие сходимости метода касательных:

\left|\frac{f(x)f''(x)}{(f'(x))^2}\right| \le q < 1.

Сформулируем теорему о квадратичной скорости сходимости метода касательных:

Теорема. Пусть x* - простой вещественный корень уравнения f(x) = 0, а функция f(x) - дважды дифференцируема в некоторой окрестности U_r(x*), причем первая произодная нигде не обращается в нуль.

Тогда, следуя обозначениям

0 < m_1 = \;\inf_{x\in U_r(x*)}|f'(x)|,\;\; M_2 = \;\sup_{x\in U_r(x*)}|f''(x)|,

при выборе начального приближения x^0 из той же окрестности U_r(x^*) такого, что

\frac{M_2|x^0 - x*|}{2m_1} = q < 1,

итерационная последовательность

x^{k+1} = x^k - \frac{f(x^k)}{f'(x^k)},\; k = 0,1,\ldots

будет сходиться к x*, причем для погрешности на k-м шаге буддет справедлива оценка:

|x^k - x*| \le q^{2^k - 1}|x^0 - x*|.

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

Рассмотрим функцию f(x)=
</p><p>== Рекомендации программисту ==
== Заключение ==
== Ссылки ==
</p>
<ul><li> [[Практикум ММП ВМК, 4й курс, осень 2008|Практикум ММП ВМК, 4й курс, осень 2008]]
</li><li> [[Метод секущих]]
</li></ul>
<p>== Список литературы ==
</p>
<ul><li>[http://mmphome.1gb.ru/index.php?pid=show&id=79 Н.В.Соснин.  Численные методы. Конспект лекций (сост. Д.В.Ховратович, Е.А.Попов)]
</li><li>Самаский А.А., Гулин А.В. Численные Методы. Учеб. пособие для вузов. - М.:Наука, 1989.
</li></ul>
<p><!--{{Stub|}}-->

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