|
|
| (24 промежуточные версии не показаны) |
| Строка 1: |
Строка 1: |
| - | == Введение ==
| |
| | | | |
| - | === Постановка математической задачи ===
| |
| - |
| |
| - | Одной из основных задач численного анализа является задача об интерполяции функций.
| |
| - | Пусть на отрезке <tex>a\le \x\le \b</tex> задана сетка <tex>\omega=\{x_i:x_0=a<x_1<\cdots<x_i<\cdots<x_n=b\}</tex> и в её узлах заданы значения функции <tex>y(x)</tex>, равные <tex>y(x_0)=y_0,\cdots,y(x_i)=y_i,\cdots,y(x_n)=y_n</tex>. Требуется построить ''интерполянту'' — функцию <tex>f(x)</tex>, совпадающую с функцией <tex>y(x)</tex> в узлах сетки:
| |
| - | {{ eqno | 1 }}
| |
| - | <p align="center"><tex>f(x_i)=y_i, i=0,1,\cdots,n.</tex></p>
| |
| - |
| |
| - | Основная цель интерполяции — получить быстрый (экономичный) алгоритм вычисления значений <tex>f(x)</tex> для значений <tex>x</tex>, не содержащихся в таблице данных.
| |
| - |
| |
| - | Интерполируюшие функции <tex>f(x)</tex>, как правило строятся в виде линейных комбинаций некоторых элементарных функций:
| |
| - | <p align="center"><tex>f(x)=\sum_{k=0}^N {c_k\Phi_k(x)},</tex></p>
| |
| - |
| |
| - | где <tex>\{\Phi_k(x)\}</tex> — фиксированный линейно независимые функции, <tex>c_0, c_1, \cdots, c_n</tex> — не определенные пока коэффициенты.
| |
| - |
| |
| - | Из условия {{eqref|1}} получаем систему из <tex>n+1</tex> уравнений относительно коэффициентов <tex>\{c_k\}</tex>:
| |
| - | <p align="center"><tex>\sum_{k=0}^N {c_k\Phi_k(x_i)}=y_i, i=0,1,\cdots,n.</tex></p>
| |
| - |
| |
| - | Предположим, что система функций <tex>\Phi_k(x)</tex> такова, что при любом выборе узлов <tex>a<x_1<\cdots<x_i<\cdots<x_n=b</tex> отличен от нуля определитель системы:
| |
| - | <p align="center"><tex>\Delta(\Phi) = \begin{vmatrix} \Phi_0(x_0) & \Phi_1(x_0) & \cdots & \Phi_n(x_0) \\\Phi_0(x_1) & \Phi_1(x_1) & \cdots & \Phi_n(x_1)\\ \hdotsfor{4} \\\Phi_0(x_n) & \Phi_1(x_n) & \cdots & \Phi_n(x_n)\end{vmatrix}</tex>.</p>
| |
| - |
| |
| - | Тогда по заданным <tex>y_i (i=1,\cdots, n)</tex> однозначно определяются коэффициенты <tex>c_k (k=1,\cdots, n)</tex>.
| |
| - |
| |
| - | == Изложение метода ==
| |
| - | Интерполяция кубическими сплайнами является частным случаем кусочно-полиномиальной интерполцией. В этом специальном случае между любыми двумя соседними узлами функция интерполируется кубическим полиномом. его коэффициенты на каждом интервале определяются из условий сопряжения в узлах:
| |
| - | <p align="center"><tex>\begin{multline}f_i=y_i,
| |
| - |
| |
| - | f'(x_i-0)=f'(x_i+0),
| |
| - |
| |
| - | f''(x_i-0)=f''(x_i+0), i=1, 2, \cdots, n-1.\end{multiline}</tex></p>
| |
| - |
| |
| - | Кроме того, на границе при <tex>x=x_0</tex> и <tex>x=x_n</tex> ставятся условия
| |
| - | {{ eqno | 2 }}
| |
| - | <p align="center"><tex>f''(x_0)=0, f''(x_n)=0.</tex></p>
| |
| - |
| |
| - | Будем искать кубический полином в виде
| |
| - |
| |
| - | {{ eqno | 2 }}
| |
| - | <p align="center"><tex>f''(x_0)=0, f''(x_n)=0.</tex></p>
| |
| - |
| |
| - | === Формула прямоугольников ===
| |
| - |
| |
| - | [[Изображение:Integral_rect.png|thumb|200px|Рис.2]]
| |
| - |
| |
| - | Заменим интеграл {{eqref|3}} выражением <tex>f(x_{i-1/2})h</tex>, где <tex>x_{i-1/2}=x_{i}-0.5h.</tex>
| |
| - |
| |
| - | Тогда получим формулу
| |
| - |
| |
| - | {{ eqno | 5 }}
| |
| - | <p align="center"><tex>\int_{x_{i-1}}^{x_i}{f(x)dx}\approx f(x_{i-1/2})h,</tex></p>
| |
| - |
| |
| - | которая называется ''формулой прямоугольников на частичном отрезке'' <tex>[x_{i-1},x_i].</tex>
| |
| - |
| |
| - | Погрешность метода {{eqref|5}} определяется величиной
| |
| - |
| |
| - | <p align="center"><tex>\psi_{i}=\int_{x_{i-1}}^{x_i}{f(x)dx}-f(x_{i-1/2})h</tex></p>
| |
| - |
| |
| - | которую легко оценить с помощью формулы Тейлора. Действительно, запишем <tex>\psi_{i}</tex> в виде
| |
| - |
| |
| - | {{ eqno | 6 }}
| |
| - | <p align="center"><tex>\psi_{i}=\int_{x_{i-1}}^{x_i}{(f(x)-f(x_{i-1/2}))dx}</tex></p>
| |
| - |
| |
| - | и воспользуемся разложением
| |
| - |
| |
| - | <p align="center"><tex>f(x)=f(x_{i-1/2})+(x-x_{i-1/2})f'(x_{i-1/2})+\frac{(x-x_{i-1/2})^2}{2}f''(\xi),</tex></p>
| |
| - |
| |
| - | где <tex>\xi_i=\xi_i(x)\in [x_{i-1},x_i]</tex>. Тогда из {{eqref|6}} получим
| |
| - |
| |
| - | <p align="center"><tex>\psi_{i}=\int_{x_{i-1}}^{x_i}{\frac{(x-x_{i-1/2})^2}{2}f''(\xi_i)dx}</tex></p>
| |
| - |
| |
| - | Обозначая <tex>M_{2,i}=\underset{x\in [x_{i-1},x_i]}{max}|f''(x)|</tex>, оценим <tex>\psi_i</tex> следующим образом:
| |
| - |
| |
| - | <p align="center"><tex>|\psi_i|\le M_{2,i} \int_{x_{i-1}}^{x_i}{\frac{(x-x_{i-1/2})^2}{2}dx}=\frac{h^3}{24}M_{2,i}</tex></p>
| |
| - |
| |
| - | Таким образом, для погрешности формулы прямоугольников на частичном отрезке справедлива оценка
| |
| - |
| |
| - | {{ eqno | 7 }}
| |
| - | <p align="center"><tex>|\psi_i|\le \frac{h^3}{24}M_{2,i}</tex></p>
| |
| - |
| |
| - | т.е. формула имеет погрешность <tex>O(h^3)</tex> при <tex>h\rightarrow0</tex>.
| |
| - |
| |
| - | Заметим,что оценка (7) является неулучшаемой, т.е. существует функция <tex>f(x)</tex>, для которой (7) выполняется со знаком равенства. Действительно, для <tex>f(x)=(x-x_{i-1/2})^2</tex> имеем <tex>M_{2,i}=2, f(x_{i-1/2})=0</tex> и
| |
| - |
| |
| - | <p align="center"><tex>\int_{x_{i-1}}^{x_i}{f(x)dx}-f(x_{i-1/2})h=\frac{h^3}{24}M_{2,i}</tex></p>
| |
| - |
| |
| - | Суммируя равенства (5) по <tex>i</tex> от <tex>1</tex> до <tex>N</tex>, получим ''составную формулу прямоугольников''
| |
| - |
| |
| - | {{ eqno | 8 }}
| |
| - | <p align="center"><tex>\int_{a}^{b}{f(x)dx}\approx \sum_{i=1}^N{f(x_{i-1/2})h}</tex></p>
| |
| - |
| |
| - | Погрешность этой формулы
| |
| - |
| |
| - | <p align="center"><tex>\Psi=\int_{a}^{b}{f(x)dx}-\sum_{i=1}^N{f(x_{i-1/2})h}</tex></p>
| |
| - |
| |
| - | равна сумме погрешностей по всем частичным отрезкам,
| |
| - |
| |
| - | <p align="center"><tex>\Psi=\sum_{i=1}^N{\psi_i}=\sum_{i=1}^N{\int_{x_{i-1}}^{x_i}{\frac{(x-x_{i-1/2})^2}{2}f''(\xi_i)dx}}</tex></p>
| |
| - |
| |
| - | Отсюда, обозначая <tex>M_2=\underset{x\in [a,b]}{max}|f''(x)|</tex>, получим
| |
| - |
| |
| - | {{ eqno | 9 }}
| |
| - | <p align="center"><tex>|\Psi|\le\frac{M_2Nh^3}{24}=\frac{h^2(b-a)}{24}M_2</tex></p>
| |
| - |
| |
| - | т.е. погрешность формулы прямоугольников на всем отрезке есть велицина <tex>O(h^2)</tex>.
| |
| - |
| |
| - | В этом случае говорят, что квадратурная формула имеет ''второй порядок точнотси''.
| |
| - |
| |
| - | === Формула трапеций ===
| |
| - |
| |
| - | [[Изображение:Integral_trap.png|thumb|200px|Рис.3]]
| |
| - |
| |
| - | На частичном отрезке эта формула имеет вид
| |
| - |
| |
| - | {{ eqno | 10 }}
| |
| - | <p align="center"><tex>\int_{x_{i-1}}^{x_i}{f(x)dx}\approx \frac{f(x_{i-1})+f(x_i)}{2}h</tex></p>
| |
| - |
| |
| - | и получается путем замены подынтегральной функции <tex>f(x)</tex> интерполяционным многочленом первой степени,постоенным по узлам <tex>x_{i-1},x_i</tex>, т.е. функцией
| |
| - |
| |
| - | <p align="center"><tex>L_{1,i}(x)=\frac{1}{h}((x-x_{i-1})f(x_i)-(x-x_i)f(x_{i-1})).</tex></p>
| |
| - |
| |
| - | Для оценки погрешности достаточно вспомнить,что
| |
| - |
| |
| - | <p align="center"><tex>f(x)-L_{1,i}(x)=\frac{(x-x_{i-1})(x-x_i)}{2}f''(\xi_i(x)).</tex></p>
| |
| - |
| |
| - | Отсюда получим
| |
| - |
| |
| - | <p align="center"><tex>\psi_i=\int_{x_{i-1}}^{x_i}{f(x)dx}-\frac{f(x_{i-1})+f(x_i)}{2}h=\int_{x_{i-1}}^{x_i}{(f(x)-L_{1,i}(x))dx}=\int_{x_{i-1}}^{x_i}{\frac{(x-x_{i-1})(x-x_i)}{2}f''(\xi_i(x))dx}</tex></p>
| |
| - |
| |
| - | и,следовательно,
| |
| - |
| |
| - | {{ eqno | 11 }}
| |
| - | <p align="center"><tex>|\psi_i|\le \frac{M_{2,i}h^3}{12}</tex></p>
| |
| - |
| |
| - | Оценка (11) неулучшаема, так как в ней достигается равенство, например, для <tex>f(x)=(x-x_i)^2</tex>.
| |
| - |
| |
| - | ''Составная формула трапеций'' имеет вид
| |
| - |
| |
| - | {{ eqno | 12 }}
| |
| - | <p align="center"><tex>\int_{a}^{b}{f(x)dx}\approx \sum_{i=1}^N{\frac{f(x_{i-1})+f(x_i)}{2}h}=h(0.5f_0+f_1+...+f_{N-1}+0.5f_N),</tex></p>
| |
| - |
| |
| - | где <tex>f_i=f(x_i),i=0,1,...,N,hN=b-a</tex>.
| |
| - |
| |
| - | Погрешность этой формулы оценивается следующим образом:
| |
| - |
| |
| - | {{ eqno | 13 }}
| |
| - | <p align="center"><tex>|\Psi|\le \frac{h^2(b-a)}{12}M_2,M_2=\underset{x\in [a,b]}{max}|f''(x)|</tex></p>
| |
| - |
| |
| - | Таким образом, формула трапеций имеет, так же как и формула прямоугольников, второй порядок точности,<tex>\Psi=O(h^2)</tex>, но ее погрешность оценивается величиной в два раза большей.
| |
| - |
| |
| - | == Числовой пример ==
| |
| - | Вычислим по формулам прямоугольников и трапеций при <tex>n=2</tex> интеграл
| |
| - |
| |
| - | {{ eqno | 14 }}
| |
| - | <p align="center"><tex>I=\int_{0}^{\pi/2}{\sin(x)dx} = 1</tex></p>
| |
| - |
| |
| - | В данном случае
| |
| - |
| |
| - |
| |
| - | <p align="center"><tex>P_2=\frac{\pi}{4}(\sin(\frac{\pi}{8})+\sin(\frac{3\pi}{8}))=1.026172</tex></p>
| |
| - |
| |
| - | <p align="center"><tex>T_2=\frac{\pi}{4}(\frac{1}{2}\sin(0)+\sin(\frac{\pi}{4})+\frac{1}{2}\sin(\frac{\pi}{2}))=0.948059</tex></p>
| |
| - |
| |
| - | Зная точный ответ (14), найдем погрешности
| |
| - |
| |
| - | {{ eqno | 15 }}
| |
| - | <p align="center"><tex>\alpha_2=-0.026172,\beta_2=0.051941</tex></p>
| |
| - |
| |
| - | Вторая производная функции <tex>\sin(x)</tex> на отрезке <tex>[0,\pi/2]</tex> отрицательна, ее модуль не превышает единицы: <tex>M_2=1</tex>. Величина погрешностей (15) удовлетворяет неравенствам (9) и (13):
| |
| - |
| |
| - | <p align="center"><tex>|\alpha_2|\le \frac{1}{96}(\frac{\pi}{2})^3<0.041,|\beta_2|\le \frac{1}{48}(\frac{\pi}{2})^3<0.081</tex></p>
| |
| - |
| |
| - | == Рекомендации программисту ==
| |
| - |
| |
| - | === Автоматический выбор шага интегрирования ===
| |
| - |
| |
| - | Величина погрешности численного интегрирования зависит как от шага сетки <tex>h</tex>, так и от гладкости подынтегральной функции <tex>f(x)</tex>. Например, в оценку (11), наряду с <tex>h</tex>, входит величина
| |
| - |
| |
| - | <p align="center"><tex>M_{2,i}=\underset{x\in [x_{i-1},x_i]}{max}|f''(x)|,</tex></p>
| |
| - |
| |
| - | которая может сильно меняться от точки к точке и, вообще говоря, заранее неизвестна. Если величина погрешности велика, то ее можно уменьшить путем измельчения сетки на данном отрезке <tex>[x_{i-1},x_i]</tex>. Для этого прежде всего надо уметь апостериорно, т.е. после проведения расчета, оценивать погрешность.
| |
| - |
| |
| - | Апостериорную оценку погрешности можно осуществить ''методом Рунге''. Пусть какая-то квадратурная формула имеет на частичном отрезке порядок точности <tex>m</tex>, т.е. <tex>I_i-I_{h,i}\approx c_i h_i^m</tex>. Тогда
| |
| - |
| |
| - | <p align="center"><tex>I_i-I_{h/2,i}\approx c_i (\frac{h_i}{2})^m,</tex></p>
| |
| - |
| |
| - | откуда получим
| |
| - |
| |
| - | {{ eqno | 16 }}
| |
| - | <p align="center"><tex>I_i-I_{h,i}\approx 2^m (I_i-I_{h/2,i}),</tex></p>
| |
| - |
| |
| - | {{ eqno | 17 }}
| |
| - | <p align="center"><tex>I_i-I_{h/2,i}\approx \frac{I_{h/2,i}-I_{h,i}}{2^m-1}</tex></p>
| |
| - |
| |
| - | Возможность апостериорно оценивать погрешность позволяет вычислять интеграл (1) с заданной точностью <tex>\epsilon >0</tex> путем автоматического выбора шага интегрирования <tex>h_i</tex>. Пусть используется составная квадратурная формула
| |
| - |
| |
| - | <p align="center"><tex>I\approx I_h=\sum_{i=1}^N{I_{h,i}}</tex></p>
| |
| - |
| |
| - | где <tex>I_{h,i}</tex> - квадратурная сумма на частичном отрезке, причем на каждом частичном отрезке используется одна и та же квадратурная формула (например, формула трапеций). Проведем на каждом частичном отрезке <tex>[x_{i-1},x_i]</tex> все вычисления дважды, один раз - с шагом <tex>h_i</tex> и второй раз - с шагом <tex>0.5h_i</tex> и оценим погрешность по правилу Рунге (17).
| |
| - |
| |
| - | Если для заданного <tex>\epsilon >0</tex> будут выполняться неравенства
| |
| - |
| |
| - | {{ eqno | 18 }}
| |
| - | <p align="center"><tex>|I_i-I_{h/2,i}|\approx \frac{|I_{h/2,i}-I_{h,i}|}{2^m-1} \le \frac{\epsilon h_i}{b-a},i=1,2,...,N,</tex></p>
| |
| - |
| |
| - | то получим
| |
| - |
| |
| - | <p align="center"><tex>|I-I_h/2|\le \frac{\epsilon}{b-a}\sum_{i=1}^N{h,i}=\epsilon,</tex></p>
| |
| - |
| |
| - | т.е. будет достигнута заданная точность <tex>\epsilon</tex>.
| |
| - |
| |
| - | Если же на каком-то из частичных отрезков оценка (18) не будет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида (18). Заметим, что для некоторой функции <tex>f(x)</tex> такое измельчение может продолжаться слишком долго. Поэтому в соответствующей программе следует предусмотреть ограничение сверху на число измельчений,а также вожможность увеличения <tex>\epsilon</tex>.
| |
| - |
| |
| - | Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции <tex>f(x)</tex> и с мелким шагом - на участках быстрого изменения <tex>f(x)</tex>. Это позволяет при заданной точности <tex>\epsilon</tex> уменьшить количество вычислений значений <tex>f(x)</tex> по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем, что для нахождения сумм <tex>I_{h/2,i}</tex> не надо пересчитывать значения <tex>f(x)</tex> во всех узлах, достаточно вычислять <tex>f(x)</tex> только в новых узлах.
| |
| - |
| |
| - | == Заключение ==
| |
| - |
| |
| - | == Список литературы ==
| |
| - |
| |
| - | * ''А.А.Самарский, А.В.Гулин.'' Численные методы М.: Наука, 1989.
| |
| - | * ''А.А.Самарский.'' Введение в численные методы М.: Наука, 1982.
| |
| - |
| |
| - | {{stub}}
| |
| - | [[Категория:Численное интегрирование]]
| |
