Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 16
Текст из файла (страница 16)
Один из возможных способов оценки ошибки при интерполяции заключается в следующем. Предполагают, что заданные числа у,, являются в действительности значениями некоторой функции ! для данных абсцисс хо Положим, что !" имеет для всех х пц-! непрерывных производных. Пусть р„— единственный полипом степени . и, интерполирующий заданные точки ((хп 4ЗЕ ВЫЧИСЛЕНИЕ ПОЛИНОМОВ 83 у,)), Тогда можно показать, что для любого действительного х а где Ч вЂ” некоторая неизвестная точка на интервале, определяемом точками х„,..., х„и х. Когда известны границы для (ь" и (х), этот результат дает оценку ошибки. Существует много обобщений лагранжевой интерполяции; наиболее употребительна среди ннх зрмипюва интерполяция..
Здесь фиксируются п абсцисс (х;), и заданных значений (у;) и п заданных угловых коэффициентов (у)). Задача состоит в том, чтобы найти полином Р(х) максимальной степени 2п — 1, такой, что для 4=1, 2, ..., п Р(х;)=у4 и Р' (х;)=у). Опять-таки, если все х; различны, то существует единственное решение, и оно может быть построено способом, вполне аналогичным методу Лагранжа. л)ы отсылаем читателя к книгам Вендрофф (!966) или Рэлстон (1966). Помимо вопросов глобальной сходимости, полиномнальная интерполяция имеет и другие недостатки.
Время построения и вычисления интерполяционных полииомов высокой степени может для некоторых приложений оказаться чрезмерным. Полиномы высокой степени могут приводить также к трудным проблемам, связанным с ошибками округлений. 4.2. Вычисление лолиномое Некоторые программы многократно вычисляют определенные полиномы для различных значений их аргументов.
Поэтому важно уметь быстро вычислять полиномы. Вычисление полннома Р (х) =- а,х" + п4х" ' +... + а„+1 посредством фортран-текста Р = А(И+1) 1ЭО Ш1=1, Ы Р = Р + А(1)4Х44(1Ч-1+1) 10 С01ЧТ11Ч11Е требует п(и+!)/2 умножений и п сложений. Простой метод, называемый схемой Горнера, состоит в перезаписи Р(х) в виде Р(Х)=а„э,+Х(а„+Х(а„,+Х(... (а4+а,Х)...))). 4. иитееполяция Это легко программируется. Например, Р А(1) ВО Ш1=1, )Ч Р РеХ+ А(1+1) 10 СО!ЧТ1(ЧОЕ Схема Горнера требует только и умножений и п плавающих сложений. Тот же алгоритм в книгах гю алгебре называется также синтстичесним делением или синтетической подстановкой, Имя Горнера присвоено методу потому, что он изложил это правило в знаменитой статье.(на другую тему), опубликованной в 1819 г. В действительности более чем за сто лет до этого идея переупорядочения была опубликована Исааком Ньютоном. Известно, что схема Горнера является оптимальной среди методов, переупорядочнвающнх полипом для быстрого вычисления и прн этом не делающих значительных вычислений в процессе переупорядочения.
Таким образом, если заданы коэффициенты полипома и аргумент х, то, вообще говоря, нельзя вычислить полинам за меньшее число сложений и умножений, чем для схемы Горнера. Интересное обсуждение вопроса о вычислении полиномов и некоторые обобщения схемы Горнера можно найти в книге Кнут (1969). 4.3. Пример: фуккцмя Рукге Опасности, связанные с полиномиальной интерполяцией, впервые обнаружил Рунге в 1901 г. Он пытался интерполировать полиномами на йнтервале ( — 1, 11 простую функцию 1 у (х) =, при равномерном распределении абсцисс. Выяснилось, что при бесконечном увеличении порядка н интерполяционного полпнома р„последовательность р„(х) расходится в интервалах 0.726...
=1х~(!. Это явление графически показано на рнс. 4.1, Отметим, что глобальная полиномиальная интерполяция дает в примере Рунге оцень хорошее согласие в центральной части интервала. Это привело к идее об интерполировании посредством движущегося иолинома. Например, можно провести полинам 10-й степени через 11 последовательных точек, но использовать его значения только в центральной части этого интервала. Поскольку сплайи-интерполяция дает много лучшие результаты, мы не станем обсуждать подробней идею движущихся полиномов. 4.а. пРимеР: ФУНкпия РУНГЕ: Если заданные абсциссы распределены не равномерно, а так, что вблизи концов интервала они помещены в корни чебышев- ского полинома степеяи и+1, то проблема расходимости для функции Рунге исчезает. Получающиеся интерполяционные многочлены рл(х) сходятся прн и- со к д(х) для всех х из [ — 1, 11.
Разумеется, этот трюк с помещением задаваемых точек в корни чебышевского полпнома работает не для,всех непрерывных 20 — е смепень 5-н спгепень 0.50.. 0.25 0.25 0.50 Ол -0.50 — 0.25 -0.50 Рвс. 4Л. Функция Рунге, интерполированная чногояленаын 5-и н 20-Я степенн прн ранноудаленных абсцнссах. функций. Теорема Фабера отрицает существование какой-либо схемы, подходящей для всех случаев. Однако можно показать, что если функция ~ имеет непрерывную производную в [ — 1, 11, то иитерполяционные многочлены р„, совпадающие с 7 в корнях полиномов Чебышева степени п-,-1, сходятся при и к [ для любой точки х из [ — 1, 1). Вопрос о чебышевских полиномах хорошо изложен в книгах Ланцош (1966), стр.
178, и Вендрофф (1966), стр. 63. вб л инТеРпОЛЯциЯ 4.4. Сплайн-интерполяция Кубические снлайнлфункции — это недавнее математическое изобретение, но они моделируют очень старое механическое устройство. Чертежника издавна пользовались механическими сплайнамп, представляющими собой гибкие рейки из 'какого- нибудь упрутого материала, обычно дерева или (в последнее время) пластика. Механический сплайн закрепляют, подвешивая грузила в точках интерполяции, называемых исторически узлами. Сплайн принимает форму, минимизирующую его потенциальную энергию, и в теории балок устанавливается, что эта энергия пропорциональна интегралу по длине дуги от квадрата кривизны сплайна, Если сплайи представить функцией з(х), то прн малых наклонах вторая производная зл(х) приблизительно равна кривизне, а дифференциал длины дуги можно приближенно заменить на с!х.
Таким образом, энергия подобного линеаризоеанного снлайна пропорциональна интегралу ) (зл(х))'с(х. Если заданы узлы (х„у,), (х„ул), ..., (х„, у„), то линеаризованиый сплайн е(х) есть функция, для которой з(х;)=у; ((=-1, 2, ..., и), и при лл этом интеграл ~ (з'(х))лс(х имеет минимальное значение, Поскольку механический сплайн не разрушается, то следует считать, что з и з' непрерывны на (х„х„). Далее, элементарная теория балок показывает, что з(х) является кубическим полиномом между каждой соседней парой узлов и что соседние поли- номы соединяются непрерывно, так же как и их первые и вторые производные.
Кубическая сплайи-функция, удовлетворя>ощая условиям зл(х,)=-.зл(х„)=0, называется естестаенным кубическим сплайном. С математической точки зрения (Алберг и др. (!967)) доказано, что она является единственной функцией, обладающей свойством минимальной кривизны, среди всех функций, интерполирующих данные точки и имеющих квадратично интегрируемую втору|о производную.
В этом смысле естественный кубический сплайн есть самая гладкая из функций, интерполирующих заданные точки. Стоит подсчитать число параметров в кубической сплайнфункции. Соответственно и†! интервалам между узлами имеется и†! фрагментов кубических кривых, у каждой из них четыре параметра„ всего, таким образом, нужно определить 4п — 4 параметров. Тот факт, что функция з и ее первая и вторая производные непрерывны во всех и — 2 внутренних узлах х;, равносилен 3(н — 2) условиям на е. Далее„ требование, чтобы е(х;)=-у; для 4.4.
сплхйн-интегполяция зг каждого из и узлов, накладывает еще и условий на з; в общей сложности пока получается 4и — 6 условий. Следовательно, нужно еще два условия, чтобы однозначно определить сплайн, Полагая з" (х,) =з" (х„) =О, мы и придем к упомянутому выше естественному сплайну. Иногда при интерполировании посредством кубических сплайнов вместо естественных граничных условий накладывают. какие-либо другие два требования на концах или вблизи концов сплайна: например, предписывают значения наклонов з'(х,) и з'(х„). Такие кубические сплайны минимизируют линеаризованный интеграл энергии при соблюдении наложенных двух условий.
Построение кубического сплайна — простой и численно устойчивый процесс. Рассмотрим подынтервал (х4, х4441, и пусть 64 Х44-1 Хз х х' П4 = —— 84 п4 = 1 — и4. Когда х пробегает этот подынтервал, п4 изменяется от О до 1, а п4 — от 1 до О. Немножко интуиции, и мы придем к представлению сплайна на этом подьштервале формулой з (х) = п4у,+ 4 + й4у4+ й,' [(и' — н4) а;„+ (в* — щ) а,1, где п4 и а4„; — некоторые константы, которые еще предстоит определить.
Первые два члена этого выражения соответствуют стандартной линейной интерполяции, а взятый в скобки член есть кубическая поправка, которая обеспечит нам дополнительную гладкость. Заметим, что поправочный член обращается в нуль на концах подынтервала, так что з(х;) =у; и З(Х4-~4)=У4А4 Таким образом, з(х) иитерполирует заданные значения независимо от выбора чисел о4. Проднфференцируем теперь з(х) трижды как сложную функцию от х, учитывая, что п4'=1~64 и ю'= — 1/й;: з (х)=- ' „'+й;[(Зщ4 — !)П44,— (34эА — 1)п41, з" (х) = бп4а; „., + 6юпц з '(х)= б (а; 4 — а4) э.интвгполяция Заметим, что э"(х) — линейная функция, интерполируюшая значения ба; и богем Следовательно, 5 (Х1) а = — ' 6 Это проясняет смысл коэффициента а,, но не позволяет определить его величину.
Отметим еше, что э"' (х) является константой на каждом подынтервале и что четвертая производная функции э(х) тождественно ранна нулю. Так, конечно, н должно быть, поскольку локально э(х) совпадает с кубической кривой. Вычисление э'(х) в концевых точках подынтервала дает э~ (х;) =- Л, — й; (а; „, + 2о;), где Л;=(у;~,— у;),'йь Приходится временно пользоваться обозначениями э', и э' „так как наша формула для э(х) имеет силу только на [х;, х,„,[, и потому производные в концевых точках являются односторонними.
Чтобы получить желаемую непрерывность э'(х), наложим во внутренних узлах условия Б' (хг) =- э,' (х;), 1 =- 2, ..., и — 1. [Непрерывность э" (х) следует непосредственно из принятого представления для э(х).! Хотя значение э' (х;) выводится из рассмотрения подынтервала [х, „х;[, формула для него получается простой заменой 1 на 1 — 1 в выражении для У (х, „,). Это приводит к равенству Л;,+й;, (2о,+а,,) =Л,— й; (ог ы+2о,) и, следовательно, й;,а;,+2(й;;+й;) о;+й,а,.~,=Ь; — Л, „1=2, „., л — 1. Это система из и — 2 линейных уравнений относительно неизвестных коэффициентов а,, 1=-1,..., л. Нужно указать еще два условия, чтобы однозначно определить интерполяционный сплайн. Среди нескольких различных способов выбора этих двух условий мы предпочли следующий.
Пусть с,(х) и с„(х) — единственные кубические кривые„которые проходят соответственно через четыре первые и четыре последние из заданных точек. Два граничных условия связывают третьи производные функции з(х) с третьими производными этих кубических кривых, именно РФ э (х,) .= с,'" и э"' (х„) -= с„'". 4 4. СПЛАйн-ИНТЕРПОЛЯЦИЯ аэ Константы с,"' и г„"' можно определить прямо из данных задачи, минуя действительное нахождение с,(х) и с„(х). Мы уже ввели величины У4 1 — Ю хз+,— х4 ' являющиеся приближениями к значениям первой производной. Пусть азз, аз 1 х,41 — х; гггг ~гг Л41, х;+ з — хз Эти величины называются разделенными разностями; при этом 2Л';" и 6Л';з' аппраксимируют соответственно вторую и третью производные. В частности, С'," = бЛ1 ' гзз ггз-1 Л~зг з-з.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
