Кафедра теплофизики и информатики в металлургии
В. С. Швыдкий, Ф.В. Куваев
Сборник задач
По дисциплине "Дифференциальные уравнения
и численные методы" для студентов специальности
071900 "Информационные системы в технике и технологиях"
всех форм обучения
Екатеринбург 2001
ВВЕДЕНИЕ
Информатика как научная дисциплина определяет методологические принципы информационного моделирования реальной действительности и манипулирования такими моделями средствами вычислительной техники. В контексте этого определения основной задачей инженера-информатика как специалиста является разработка информационных моделей объектов реального мира и получения информации об объекте путём реализации этих моделей программными средствами. Учитывая технологическую направленность специальности 071900, последнее утверждение подразумевает разработку и реализацию математических моделей технологических объектов.
Технические и технологические объекты чёрной и цветной металлургии, теплоэнергетики, химии и т.д. (в принципе, любые реальные объекты) характеризуются взаимосвязанностью и взаимообусловленностью различных параметров, определяющих процесс. Вследствие этого математические модели таких объектов в математическом и физическом планах являются существенно нелинейными, и уравнения этих моделей не могут быть решены средствами традиционной высшей математики. Рабочим инструментом инженера-информатика являются методы и средства вычислительной (или, как ещё говорят, прикладной) математики. Задачи и упражнения, приведённые в данном пособии, предназначены для закрепления теоретических положений основных численных методов, составляющих значительную часть методов прикладной математики.
ИНТЕРПОЛИРОВАНИЕ
При использовании численных методов часто приходится в целях упрощения расчётов (или по другим причинам) заменять некоторые функции y (x) их упрощенными (приближёнными) зависимостями. Отсюда вытекает задача: выполнить эту замену таким образом, чтобы приближённая зависимость была как можно более близкой к исходной функции. Именно для решения сформулированной задачи и предназначены методы интерполяции.
Напомним некоторые положения теории. Пусть функция y (x) известна только в узлах некоторой сетки xi, т.е. задана таблицей. В этих случаях выгодно заменить функцию y (x) (аппроксимировать) приближённой формулой, т.е. подобрать некоторую другую функцию j(х), которая близка в некотором смысле к у (х) и просто вычисляется. Затем при всех значениях аргумента x можно полагать, что y (x)» j(x). Близость функций обеспечивают введением в аппроксимирующую функцию свободных параметров a = { a 1, a 2, …, an } и соответствующим их выбором.
Процедура выбора свободных параметров а определяет тип метода интерполяции, и здесь можно отметить нижеприведённые наиболее распространённые методы.
Линейная интерполяция
Если потребовать, чтобы аппроксимирующая функция j(х, а) совпадала с табличными значениями исходной функции у (х) в n выбранных узлах сетки, то получим систему
j(xi, a 1, a 2, …, an) = y (xi) º yi, 1 £ i £ n, (1.1)
из которой можно определить параметры ai. Этот способ подбора параметров называется лагранжевой интерполяцией. По числу используемых узлов сетки различают интерполяцию одноточечную, двухточечную и т.д.
В простейшем случае аппроксимирующая функция j(х, а) линейно зависит от параметров, т.е. может быть представлена в виде так называемого обобщённого полинома
(1.2)
Функции j k (x) должны быть линейно-независимыми. Кроме того, на систему функций j k (x) налагается ещё одно ограничение. Подставляя (1.2) в (2.1), получим для определения параметров ak следующую систему линейных уравнений:
(1.3)
Для единственности задачи интерполяции необходимо, чтобы при любом расположении узлов (лишь бы среди них не было совпадающих) определитель системы (1.3) был бы отличен от нуля:
(1.4)
Система функций, удовлетворяющих требованию (1.4), называется чебышевской. Таким образом, при линейной интерполяции нужно строить обобщённый полином по какой-нибудь чебышевской системе функций.
Для линейной интерполяции наиболее удобны обычные полиномы, так как они легко вычисляются и на элементарном калькуляторе, и на ЭВМ.
· Интерполяционный многочлен Ньютона. Если в качестве системы функций выбрать j k (x) = xk, 0 £ k £ n, то, используя предыдущие соотношения, придём к интерполяционной формуле Ньютона
(1.5)
Здесь y (x 0, x 1, ¼, xk) – так называемая k-я разделённая разность, определяемая соотношениями
(1.6)
и т.д. Разделённые разности первого, второго и более высоких порядков имеют размерности производных соответствующих порядков. Отметим также, что во многих учебниках для таблиц с постоянным шагом h = xj – xi вводят конечные разности D ny, связанные с разделёнными разностями соотношением D ny = n! y (x 0, x 1, ¼, xn). В конечных разностях формула Ньютона (1.5) записывается так:
(1.5,a)
Формула Ньютона удобна для вычислений и на ЭВМ, и на калькуляторе. Легко составить следующую табл.1.1 разделённых разностей для табулированной функции у (х) и произвести вычисления по формуле (1.5).
Таблица 1.1. Разделённые разности для использования формулы Ньютона (1.5)
| x 0 | y (x 0) | |||
| y (x 0, x 1) | ||||
| x 1 | y (x 1) | y (x 0, x 1, x 2) | ||
| y (x 1, x 2) | y (x 0, x 1, x 2, x 3) | |||
| x 2 | y (x 2) | y (x 1, x 2, x 3) | ||
| y (x 2, x 3) | ||||
| x 3 | y (x 3) |
За точностью расчёта удобно следить, визуально оценивая скорость убывания слагаемых суммы(1.5). Если они убывают медленно, то на хорошую точность рассчитывать нельзя. Если убывание быстрое, то оставляют только те слагаемые, которые больше допустимой погрешности; тем самым определяют, сколько узлов требуется подключить в расчёт.
Пример 1.1. Покажем, как вычислять синус в первом квадранте, используя четыре стандартных значения. Составим табл. 1.2 с четырьмя узлами, причём для удобства вычисления положим у(х) = sin(p×x/6).
Таблица 1.2. К примеру 1.1
| xi | y (xi) | y (xi, xi +1) | y (xi, xi +1, xi +2) | y (xi, ¼, xi +3) |
| 0,000000 | ||||
| 0,500000 | ||||
| 0,500000 | -0,066987 | |||
| 0,366025 | -0,016346 | |||
| 0,866025 | -0,116025 | |||
| 0,133975 | ||||
| 1,000000 |
Для проверки точности, используя разности верхней косой строки, вычислим
у(1,5) = 0 + (1,5 – 0)×0,500 – (1,5 – 0) × (1,5 – 1) × 0,066987 -
- (1,5 – 0) × (1,5 – 1) × (1,5 – 2) × 0,016346 = 0,705890.
Это приближённое значение лишь на 0,17 % отличается от точного значения у(1,5) = sin(45o) = 0,707107. Таким образом, достаточно помнить только верхнюю косую строку табл. 1.2, чтобы вычислить синус с точностью около 0,002.
Интерполяцию можно использовать для решения нелинейных алгебраических уравнений. Для этого применяют процедуру обратного интерполирования. При этом для уравнения y (x) составляют небольшую таблицу значений y (xi) и записывают её в виде xi = x (yi). Далее обычным способом
находят приближённое значение х (0). Этот способ даёт хорошие результаты, если функция достаточно гладкая, а корень лежит между рассчитанными узлами.
Пример 1.2. Решим уравнение
y(x) º (1 + x) e 0,5 x – 2,5 = 0.
Составим таблицу 4-х значений функции:
х 0 0,5 1,0 1,5
у - 1,5 - 0,573962 0,797443 2,792500
Запишем первым столбец значений у, так как в дальнейших вычислениях эта величина будет аргументом. Найдём разделённые разности
Таблица 1.3. Разделённые разности к примеру 1.2
| yi | xi | x(yi, yi +1 ) | x(yi, yi +1, yi +2) |
| - 1,5 | |||
| 0,539935 | |||
| - 0,573962 | 0,5 | - 0,076322 | |
| 0,364590 | |||
| 0,797443 | 1,0 | - 0,033855 | |
| 0,250619 | |||
| 1,792500 | 1,5 |
Произведём вычисления по верхней косой строке:
x(0)» x 0 + (0 – y 0 ) × x(y 0, y 1 ) + (0 – y 0 ) × (0 – y 1 ) × x(y 0, y 1, y 2 ) =
= 0 + 1,5 × 0,539935 – 1,5 × 0,573962 × 0,076322 = 0,744193.
Точное решение есть 0,732965, так что ошибка получается равной 1,53 %, т.е. небольшой. Для повышения точности в этом способе целесообразно взять узлы, близко расположенные к грубо найденному корню, а не увеличивать число узлов.
· Численное дифференцирование. Численное дифференцирование применяется, если функцию у (х) трудно или невозможно продифференцировать аналитически – например, если она задана таблично.
При численном дифференцировании функцию у (х) аппроксимируют легко вычисляемой функцией j(х, а) и приближённо полагают y '(x) = j'(x, a). Вводя в формуле Ньютона обозначение x i = x – xi, запишем этот многочлен и продифференцируем его почленно:
Обрывая ряд на некотором числе слагаемых, получим приближённое выражение для соответствующей производной. Наиболее простые выражения получим, оставляя в предыдущих формулах только первые слагаемые:
(1.7)
Пример 1.3. Для функции, заданной табл.1.4, вычислить её значения в точках 0,04170; 0,529394; 0,707282; значение её производной в точке 0,142933 и решить уравнение у(х) = 1,509393.Результаты получить с той точностью, которую допускает табличное задание функции.
Таблица 1.4. Функция к примеру 1.3
| x | y (x) | x | y (x) |
| 0,0 | 1,534906 | 0,5 | 1,458955 |
| 0,1 | 1,517093 | 0,6 | 1,448162 |
| 0,2 | 1,500482 | 0,7 | 1,438994 |
| 0,3 | 1,485191 | 0,8 | 1,431485 |
| 0,4 | 1,471320 | 0,9 | 1,425652 |
Составим таблицу разделённых разностей функции.
Таблица 1.5. Разделённые разности заданной функции
| xi | y (xi) | y (xi, xi +1) | y (xi, xi +1, xi +2) | y (xi, …, xi +3) | y (xi, …, xi +4) |
| 0,0 | 1,534906 | ||||
| - 0,17867 | |||||
| 0,1 | 1,517093 | 0,0628 | |||
| - 0,16611 | 0,01067 | ||||
| 0,2 | 1,500482 | 0,0660 | 0,015 | ||
| - 0,15291 | 0,01667 | ||||
| 0,3 | 1,485191 | 0,0710 | - 0,005825 | ||
| - 0,13871 | 0,01433 | ||||
| 0,4 | 1,471320 | 0,0753 | 0,0000 | ||
| - 0,12365 | 0,01433 | ||||
| 0,5 | 1,458955 | 0,0796 | - 0,022075 | ||
| - 0,10793 | 0,0055 | ||||
| 0,6 | 1,448162 | 0,08125 | 0,0004175 | ||
| - 0,09168 | 0,005667 | ||||
| 0,7 | 1,438994 | 0,08295 | - 0,00625 | ||
| - 0,07509 | 0,003167 | ||||
| 0,8 | 1,431485 | 0,0838 | |||
| - 0,05833 | |||||
| 0,9 | 1,425652 |
Можно видеть, что пятые разделённые разности изменяются хаотическим образом, поэтому интерполяцию целесообразно строить лишь на основе первых четырёх разделённых разностей.
Первая точка интерполяции х = 0,041470 лежит между первым и вторым узлами сетки. Используя данные верхней косой строки, находим
y(0,041470)» 1,534906 – (0,041470 – 0) × 0,17867 + (0,041470 – 0) ´
´ (0,041470 – 0,1) × 0,0628 + (0,041470 – 0) × (0,041470 – 0,1) × (0,04147 – 0,2) × 0,01067 =
@ 1,527385,
причём погрешность результата не превышает 0,00005.
Вторая точка интерполяции х = 0,529394 лежит между пятым и шестым узлами. Здесь мы воспользуемся нисходящей косой строкой, начинающейся на горизонте х = 0,5. Тогда
y(0,529394)» 1,458955 – (0,529394 – 0,5) × 0,10793 + (0,529394 – 0,5) × (0,529394 – 0,6) ´
´ 0,08125 + (0,529394 – 0,5) × (0,529394 – 0,6) × (0,52394 – 0,7) × 0,005667 @ 1,455621.
Наконец, последняя точка интерполяции лежит между восьмым и девятым узлами. Для вычисления функции в этой точке нисходящей косой строки уже не хватает, поэтому воспользуемся восходящей косой строкой, начинающейся на горизонте х = 0,8. Имеем
y(0,707282)» 1,431485 – (0,707282 – 0,8) × 0,07509 + (0,707282 – 0,8) × (0,707282 – 0,7) ´
´ 0,08295 + (0,707282 – 0,8) × (0,707282 – 0,7) × (0,707282 – 0,6)× 0,005667 @ 1,438908.
Производную в точке х = 0,142933 вычислим по формуле Ньютона, принимая х 0 = 0,1,
y ' (0,142933)» - 0,16611 + [(0,142933 – 0,1) + (0,142933 – 0,2)] × 0,066 +
+ [(0,142933 – 0,1) × (0,142933 – 0,2) + (0,142933 – 0,1) × (0,142933 – 0,3) +
+ (0,142933 – 0,2) × (0,142933 – 0,3)] × 0,01667 @ - 0,167047.
Решение уравнения у(х) = 1,509393 организуем по методу обратной итерации. Преобразуем таблицу разделённых разностей, использовав только те строки, значения функции в которых близки к значению правой части уравнения.
| yi | x (yi) | x (yi, yi +1) | x (yi, yi +1, yi +2) | x (yi, yi +1, yi +2, yi +3) |
| 1,517093 | 0,1 | |||
| - 6,019745 | ||||
| 1,500482 | 0,2 | 16,301486 | ||
| - 6,539795 | - 145,416818 | |||
| 1,485191 | 0,3 | 22,957650 | ||
| - 7,209286 | ||||
| 1,471320 | 0,4 |
x(1,509393)» 0,1 – (1,509393 – 1,517093) × 6,019745 + (1,509393 – 1,517093) × (1,509393 – - 1,500482) × 16,301486 - (1,509393 – 1,517093) × (1,509393 – 1,500482) × (1,509393 – 1,485191) ´ ´ 145,416818 @ 0,145475.
Задание № 1
Для функции, заданной таблично, найти её значение в точках х 1, х 2, х 3, значение её производной в точке у и решить уравнение у (х) = а. Результаты получить с той точностью, которую допускает табличное задание функции.
Решение задачи оформить программно (на любом языке программирования или в любой среде). Результаты выдать на печать по образцу примера 1.3.
Вариант А
Интерполируемая функция
| x | y (x) | x | y (x) |
| 0,0 | 0,946083 | 0,5 | 1,324684 |
| 0,1 | 1,028685 | 0,6 | 1,389181 |
| 0,2 | 1,108047 | 0,7 | 1,449592 |
| 0,3 | 1,183958 | 0,8 | 1,505817 |
| 0,4 | 1,256227 | 0,9 | 1,557775 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,175418 | 0,715878 | 0,464331 | 0,817438 | 1,344273 | |
| 0,090566 | 0,826611 | 0,395142 | 0,329982 | 1,534355 | |
| 0,157369 | 0,826216 | 0,445135 | 0,045464 | 1,273228 | |
| 0,053526 | 0,866027 | 0,610308 | 0,645041 | 1,435993 | |
| 0,109658 | 0,710092 | 0,468515 | 0,199282 | 1,026087 | |
| 0,094189 | 0,755627 | 0,343888 | 0,406476 | 1,058530 | |
| 0,132747 | 0,755108 | 0,447232 | 0,398774 | 1,390474 | |
| 0,058421 | 0,731564 | 0,429815 | 0,549309 | 1,417970 | |
| 0,053946 | 0,807395 | 0,455020 | 0,406727 | 1,373800 | |
| 0,182461 | 0,717093 | 0,525286 | 0,897085 | 1,197522 | |
| 0,173762 | 0,759231 | 0,482714 | 0,248933 | 1,470054 | |
| 0,129140 | 0,766507 | 0,238812 | 0,518578 | 1,183944 |
Вариант Б
Интерполируемая функция
| x | y (x) | x | y (x) |
| 0,0 | 0,962852 | 0,5 | 1,337904 |
| 0,1 | 1,044824 | 0,6 | 1,401593 |
| 0,2 | 1,123512 | 0,7 | 1,461175 |
| 0,3 | 1,199709 | 0,8 | 1,516552 |
| 0,4 | 1,270228 | 0,9 | 1,567630 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,052914 | 0,882964 | 0,391986 | 0,541089 | 1,189116 | |
| 0,147327 | 0,721404 | 0,355038 | 0,051982 | 1,564167 | |
| 0,148489 | 0,822272 | 0,544808 | 0,707369 | 1,573997 | |
| 0,070926 | 0,801307 | 0,213966 | 0,670379 | 1,547732 | |
| 0,124548 | 0,769617 | 0,303341 | 0,473246 | 1,504725 | |
| 0,009510 | 0,844211 | 0,417157 | 0,181704 | 1,094108 | |
| 0,087996 | 0,778350 | 0,546174 | 0,291991 | 1,110465 | |
| 0,051598 | 0,709740 | 0,542017 | 0,726683 | 1,404176 | |
| 0,083288 | 0,887920 | 0,517459 | 0,367139 | 1,484642 | |
| 0,174607 | 0,802430 | 0,240569 | 0,285460 | 1,275640 | |
| 0,172537 | 0,819656 | 0,698720 | 0,830955 | 1,331830 | |
| 0,160270 | 0,851787 | 0,463563 | 0,486052 | 1,307859 |
Вариант В
Интерполируемая функция
| х | y (x) | x | y (x) |
| 0,0 | 1,848652 | 0,5 | 1,833125 |
| 0,1 | 1,851659 | 0,6 | 1,821948 |
| 0,2 | 1,851401 | 0,7 | 1,808622 |
| 0,3 | 1,848081 | 0,8 | 1,793390 |
| 0,4 | 1,841914 | 0,9 | 1,776501 |
Точки интерполирования
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,102067 | 0,783267 | 0,233571 | 0,701138 | 1,775720 | |
| 0,063358 | 0,845391 | 0,479810 | 0,016395 | 1,836223 | |
| 0,076951 | 0,757948 | 0,534155 | 0,758749 | 1,830957 | |
| 0,148489 | 0,728158 | 0,614515 | 0,492288 | 1,846149 | |
| 0,115456 | 0,778255 | 0,568848 | 0,877666 | 1,815639 | |
| 0,151273 | 0,836281 | 0,621168 | 0,319084 | 1,842761 | |
| 0,144289 | 0,757769 | 0,497664 | 0,018630 | 1,808383 | |
| 0,113056 | 0,901257 | 0,492267 | 0,070493 | 1,775438 | |
| 0,029295 | 0,895077 | 0,544469 | 0,786666 | 1,781802 | |
| 0,120383 | 0,762649 | 0,493649 | 0,297184 | 1,850931 | |
| 0,176574 | 0,770618 | 0,413157 | 0,062844 | 1,849610 | |
| 0,152489 | 0,745512 | 0,409483 | 0,386388 | 1,836708 |
Вариант Г
Интерполируемая функция
| x | y (x) | x | y (x) |
| 0,0 | 1,849524 | 0,5 | 1,831074 |
| 0,1 | 1,851862 | 0,6 | 1,819447 |
| 0,2 | 1,850975 | 0,7 | 1,805720 |
| 0,3 | 1,847068 | 0,8 | 1,790137 |
| 0,4 | 1,840359 | 0,9 | 1,772946 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,044427 | 0,823677 | 0,570277 | 0,756676 | 1,804051 | |
| 0,127215 | 0,787367 | 0,353208 | 0,669386 | 1,774512 | |
| 0,137938 | 0,886153 | 0,514296 | 0,812431 | 1,780675 | |
| 0,034406 | 0,874237 | 0,463792 | 0,697030 | 1,830502 | |
| 0,123011 | 0,739399 | 0,308672 | 0,295801 | 1,847906 | |
| 0,020533 | 0,823018 | 0,372752 | 0,848340 | 1,792606 | |
| 0,119622 | 0,790245 | 0,459663 | 0,273649 | 1,820115 | |
| 0,032266 | 0,883370 | 0,273037 | 0,644472 | 1,826426 | |
| 0,177173 | 0,725458 | 0,414551 | 0,131848 | 1,781212 | |
| 0,042385 | 0,748684 | 0,347492 | 0,511049 | 1,829010 | |
| 0,060953 | 0,888796 | 0,514699 | 0,875722 | 1,808055 | |
| 0,101105 | 0,711049 | 0,315694 | 0,705710 | 1,823896 |
Метод наименьших квадратов
Интерполяция позволяет легко аппроксимировать функцию у (х). Однако точность такой интерполяции гарантирована лишь в небольшом интервале порядка нескольких шагов сетки. Для другого интервала приходится заново вычислять коэффициенты интерполяционной формулы. В практической же деятельности всегда желательно иметь единую приближённую формулу y»» j(x), пригодную для большого отрезка a £ x £ b.
Имеется ещё одна причина стремления получить общую формулу интерполяции для всего интервала изменения х. При интерполяции мы приравниваем значения у (х) и j(х) в узлах сетки. Но если значения y (xi) определены неточно – например, из эксперимента, - то точное приравнивание неразумно. Поэтому нередко целесообразней приближать функцию не по точкам, а в среднем для всего интервала изменения х.
В этом и состоит основная идея метода наименьших квадратов. Пусть аппроксимирующая функция j(х) представляет собой множество в линейном пространстве функций
(1.8)
где функции j k (x) – линейно-независимые. В каждом i -том узле таблицы эта функция будет давать несколько отличную от значения y (xi) º yi величину, образуя узловую невязку Ri = j(xi) – yi. Так как узловые невязки могут иметь разные знаки, то суммарная погрешность формулы (1.8) будет характеризоваться выражением
(1.9)
где i = 1, 2, ¼, I – номер узла таблицы, а n характеризует степень аппроксимирующего полинома. Основное уравнение метода наименьших квадратов формулируется как требование выбрать "настроечные параметры" ak таким образом, чтобы был обеспечен минимум квадратичного отклонения, т.е.
R 2 = min. (1.10)
Очевидным следствием этого требования является система уравнений
(1.11)
В том часто используемом случае, когда j k (x) = xk, систему уравнений (1.11) можно переписать в виде
(1.12)
Здесь введена весовая функция r i – позволяющая учесть относительную значимость узлового значения функции. Очевидно, чем точнее это значение, тем больший вес оно должно иметь по сравнению с менее точными данными.
Поскольку степени на любом отрезке образуют чебышевскую систему функций, то задача (1.12) имеет единственное решение. Но система степеней не ортогональна, и при больших значениях n задача (1.12) плохо обусловлена. Чтобы обойти эту трудность, обычно ограничиваются невысокими степенями n» 2 ¸ 5, когда обусловленность задачи (1ю12) удовлетворительна.
Таким образом, задача интерполяции при описанном подходе сводится к задаче решения системы алгебраических уравнений относительно коэффициентов ak с последующей записью аналитического выражения и выполнением всех вычислений на основе этого выражения. Так как решению систем линейных алгебраических уравнений посвящено отдельное задание, то здесь рассмотрим другую возможность.
В настоящее время вопросы интерполяции в плане подбора аппроксимирующего выражения решаются элементарно. Практически в любом математическом пакете обязательно имеется команда, функция или процедура интерполяции экспериментальных данных в смысле метода наименьших квадратов. В качестве примеров использования таких команд (или функций) рассмотрим решение задачи примера 1.3 в рамках пакетов Mathcad, Maple и MATLAB.
В пакете Mathcad сначала формируются векторы узловых координат и узловых значений функции; обозначим их через Х и Y. Дальнейшая процедура показана на рис. 1.2.
Рис. 1.1. Решение задачи интерполяции в рамках пакета Mathcad
Располагая аналитическим описанием функции, теперь нетрудно ответить на все вопросы, поставленные в условиях задачи (см. рис. 1.2). Можно наглядно убедиться в близости аппроксимации к заданным узловым значениям функции, построив график, а можно вычислить коэффициент корреляции (точнее, корреляционное отношение), характеризующий степень близости исходных данных и результатов расчёта по аппроксимирующему полиному. Как видно из данных рис. 1.2, этот коэффициент оказался равен 1, т.е. полином точно описывает табличные данные.
В рамках пакета Maple задачу интерполяции можно решить двумя способами. Во-первых, для этого имеется специальная команда interp(X, Y, v), где Х и Y имеют тот же самый смысл, что и в Mathcad, а параметр v указывает имя аргумента интерполяционного полинома. По умолчанию степень полинома на единицу меньше количества вводимых узлов, т.е. если ввести все значения исходной таблицы, то пакет определит коэффициенты полинома девятой степени.
Рис. 1.2. Решение задачи интерполяции и оценка точности решения
В рамках пакета Maple имеется возможность проверить утверждение о низкой надёжности (плохой обусловленности) полиномов высокой степени. Приведём фрагмент интерактивной работы с пакетом для случая, когда заданы все данные исходной таблицы (рис. 1.3).
Можно видеть, что в данном случае, когда интерполируемая функция изменяется плавно, использование полинома высокой степени не искажает фактической зависимости. Вопрос о том, приводит ли повышение порядка интерполяционного полинома к увеличению точности интерполяции, может быть решён лишь при сопоставительных расчётах. На рис. 1.4 показаны результаты таких расчётов (по крайней мере, в графическом представлении) для тех же исходных данных, когда исключено каждое чётное узловое значение функции. Вряд ли результаты стали хуже; в то же время выполнять расчёты с полиномом четвёртого порядка намного проще.
Рис. 1.3. Решение задачи интерполяции в пакете Maple
Рис. 1.4. Решение той же задачи при меньшем числе точек
Вычисление значений функции в заданных точках интерполяции при известном уравнении функции не вызывает затруднений, поэтому на этом вопросе мы не останавливаемся. Отметим только, что решение уравнения в рамках пакета Maple также осуществляется специальной командой solve. При этом команда выдаёт вообще все решения для заданного уравнения, а не только для интервала таблицы.
В рамках пакета MATLAB задача интерполяции для всей области, заданной таблицей, решается с использованием команды polyfit(X,Y,n), где n – степень полинома. Команда возвращает коэффициенты полинома в порядке убывания степеней. Результаты работы команды сохраняются в рабочей области пакета и могут быть использованы в других командах без явного указания. Вычисления по аппроксимирующему полиному выполняются командой polyval(p, X), где р – вектор коэффициентов полинома. Здесь можно сразу выполнить оценку погрешности аппроксимации, как это показано на рис. 1.5.
Рис. 1.5. Решение задачи интерполяции в пакете MATLAB
Можно также проиллюстрировать результаты интерполяции соответствующим графиком (рис. 1.6):
>> plot(X, Y, 'ob', X, y)
Рис. 1.6. Оценка точности интерполяции
Можно задачу аппроксимации решать и в рамках программы Excel. Для этого, прежде всего, необходимо уяснить себе вычислительную суть метода наименьших квадратов. Предположим, что мы ожидаем квадратичной зависимости табулированной функции y (xi) от аргумента xi, т.е.
y (x) = a + bx + cx 2.
Чтобы обеспечить выполнение условия
R 2 = 
где yi = y (xi), необходимо подобрать неизвестные коэффициенты a, b и c из условий отыскания экстремума
откуда получаем систему уравнений
Дальнейший расчёт по этим уравнениям приведён ниже.
Задание № 2
Для функции, заданной таблично, используя любой из математических пакетов (допускается применение программы Excel Microsoft Office), найти её значения в точках х 1, х 2, х 3, значение её производной в точке у и решить уравнение у (х) = а. Решение получить с точностью до пяти знаков после запятой.
Вариант А
Интерполируемая функция
| х | у (х) | х | у (х) |
| 0,0 | 0.979498 | 0,5 | 1,350965 |
| 0,1 | 1,060831 | 0,6 | 1,413842 |
| 0,2 | 1,138837 | 0,7 | 1,472590 |
| 0,3 | 1,213312 | 0,8 | 1,527116 |
| 0,4 | 1,248076 | 0,9 | 1,577351 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,047475 | 0,805753 | 0,653213 | 0,577700 | 1,423549 | |
| 0,104834 | 0,865748 | 0,634411 | 0,220322 | 1,361976 | |
| 0,095092 | 0,861356 | 0,490247 | 0,355379 | 1,314712 | |
| 0,133412 | 0,879781 | 0,473981 | 0,821980 | 1,472912 | |
| 0,069174 | 0,793965 | 0,465591 | 0,456153 | 1,498836 | |
| 0,076449 | 0,702410 | 0,552699 | 0,020310 | 1,580131 | |
| 0,086081 | 0,711768 | 0,446462 | 0,830259 | 1,479177 | |
| 0,045216 | 0,721249 | 0,670968 | 0,064005 | 1,429159 | |
| 0,182628 | 0,729050 | 0,433820 | 0,743754 | 1,097459 | |
| 0,155091 | 0,892695 | 0,403838 | 0,000455 | 0,981061 | |
| 0,118535 | 0,884719 | 0,589080 | 0,490510 | 1,399221 | |
| 0,041623 | 0,839124 | 0,294268 | 0,902946 | 1,048028 |
Вариант Б
Интерполируемая функция
| х | у (х) | х | у (х) |
| 0,0 | 0,996021 | 0,5 | 1,363865 |
| 0,1 | 1,076705 | 0,6 | 1,425925 |
| 0,2 | 1,154020 | 0,7 | 1,483835 |
| 0,3 | 1,227768 | 0,8 | 1,537508 |
| 0,4 | 1,297768 | 0,9 | 1,586879 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,143689 | 0,915777 | 0,599657 | 0,095444 | 1,018375 | |
| 0,042666 | 0,825305 | 0,584323 | 0,573204 | 1,262410 | |
| 0,165626 | 0,705625 | 0,456150 | 0,346457 | 1,445709 | |
| 0,126742 | 0,748197 | 0,450881 | 0,777682 | 1,467961 | |
| 0,096649 | 0,801334 | 0,629973 | 0,914914 | 1,342491 | |
| 0,064256 | 0,843814 | 0,373073 | 0,759432 | 1,362658 | |
| 0,071993 | 0,742732 | 0,474374 | 0,648813 | 1,567605 | |
| 0,105362 | 0,711920 | 0,524242 | 0,160211 | 1,585812 | |
| 0,159355 | 0,819366 | 0,337411 | 0,563150 | 1,427805 | |
| 0,161147 | 0,900531 | 0,252948 | 0,609420 | 1,146269 | |
| 0,104874 | 0,852465 | 0,470731 | 0,272350 | 1,033643 | |
| 0,175694 | 0,779153 | 0,306156 | 0,767306 | 1,551797 |
Вариант В
Интерполируемая функция
| х | у (х) | х | у (х) |
| 0,0 | 1,850260 | 0,5 | 1,828928 |
| 0,1 | 1,851937 | 0,6 | 1,816862 |
| 0,2 | 1,850428 | 0,7 | 1,802744 |
| 0,3 | 1,845943 | 0,8 | 1,786820 |
| 0,4 | 1,838700 | 0,9 | 1,769337 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,070782 | 0,933883 | 0,532991 | 0,913222 | 1,847893 | |
| 0,057504 | 0,740810 | 0,300079 | 0,419130 | 1,825998 | |
| 0,156507 | 0,759066 | 0,665331 | 0,882467 | 1,808450 | |
| 0,188887 | 0,733400 | 0,521501 | 0,297946 | 1,771221 | |
| 0,199791 | 0,895346 | 0,529348 | 0,163911 | 1,775305 | |
| 0,118051 | 0,747759 | 0,394174 | 0,871601 | 1,771558 | |
| 0,117019 | 0,819226 | 0,504759 | 0,743122 | 1,790246 | |
| 0,045029 | 0,705874 | 0,513963 | 0,437270 | 1,784766 | |
| 0,164030 | 0,728025 | 0,227815 | 0,868064 | 1,811010 | |
| 0,085978 | 0,777671 | 0,494874 | 0,303219 | 1,828668 | |
| 0,097120 | 0,858195 | 0,315239 | 0,377217 | 1,799677 | |
| 0,149490 | 0,705743 | 0,487051 | 0,882789 | 1,841611 |
Вариант Г
Интерполируемая функция
| х | у (х) | х | у (х) |
| 0,0 | 1,850859 | 0,5 | 1,826690 |
| 0,1 | 1,851883 | 0,5 | 1,814195 |
| 0,2 | 1,849763 | 0,7 | 1,799696 |
| 0,3 | 1,844708 | 0,8 | 1,783440 |
| 0,4 | 1,836941 | 0,9 | 1,765676 |
Точки интерполяции
| Вариант | х 1 | х 2 | х 3 | у | а |
| 0,149737 | 0,747740 | 0,566665 | 0,332744 | 1,851871 | |
| 0,097144 | 0,850730 | 0,445502 | 0,910664 | 1,846796 | |
| 0,135971 | 0,868823 | 0,373904 | 0,511079 | 1,766612 | |
| 0,072560 | 0,847703 | 0,502299 | 0,664096 | 1,849372 | |
| 0,150472 | 0,723186 | 0,339706 | 0,133476 | 1,788283 | |
| 0,119440 | 0,826031 | 0,337938 | 0,869169 | 1,811741 | |
| 0,137092 | 0,879423 | 0,550316 | 0,820591 | 1,789388 | |
| 0,090480 | 0,851292 | 0,509452 | 0,685084 | 1,771309 | |
| 0,094842 | 0,783729 | 0,621566 | 0,160136 | 1,831897 | |
| 0,153709 | 0,752197 | 0,357420 | 0,859274 | 1,833945 | |
| 0,048033 | 0,758573 | 0,361878 | 0,418896 | 1,844793 | |
| 0,042689 | 0,831531 | 0,504781 | 0,388134 | 1,812163 |
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
Пусть требуется вычислить определённый интеграл
(2.1)
где функция f (x) непрерывна на отрезке [ a, b ], а весовая функция r(x) непрерывна на интервале (a, b).
