ЭВМ_Семестр3_МетодПособие
.pdfЗадача 2. Записать зависимости 1 (x), 2 (x), 3 (x), 4 (x), интерполирующие точки данных, показанные соответственно на рисунках. 4.4 а), б). в), г). Решение вспомогательных задач:
1 |
(x) |
|
|
(x x2 ) (x x3 ) (x x4 ) |
|
||||||||||||||||||
|
(x1 x2 ) (x1 x3 ) (x1 x4 ) |
||||||||||||||||||||||
|
|
|
|
|
|||||||||||||||||||
|
|
|
|
|
|
|
|
|
(x x ) |
|
|
|
|
|
|
|
|
||||||
|
|
|
(пропущено |
1 |
|
), |
|
|
|
|
|
|
|||||||||||
|
|
|
(x x ) |
|
|
|
|
|
|
||||||||||||||
|
|
|
|
|
|
|
|
|
1 |
|
1 |
|
|
|
|
|
|
|
|
|
|||
2 |
|
(x) |
|
|
(x x1 ) (x x3 ) (x x4 ) |
||||||||||||||||||
|
(x2 x1 ) (x2 x3 ) (x2 x4 ) |
||||||||||||||||||||||
|
|
|
|
||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
(x x2 ) |
|||||||||||||
|
|
|
( пропущено |
|
(x |
2 |
x |
2 |
) |
), |
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
3 |
(x) |
|
|
(x x1 ) (x x2 ) (x x4 ) |
|
||||||||||||||||||
|
(x3 x1 ) (x3 x2 ) (x2 x4 ) |
||||||||||||||||||||||
|
|
|
|
|
|||||||||||||||||||
|
|
|
(пропущено |
|
(x x3 ) |
), |
|
|
|
|
|||||||||||||
|
|
|
|
(x |
|
|
x ) |
|
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
3 |
|
3 |
|
|
|
|
|
|
|
|
|
|
|
4 |
(x) |
(x x1 ) (x x2 ) (x x3 ) |
|
|||||||||||||||||||
|
(x4 x1 ) (x4 x2 ) (x4 x3 ) |
||||||||||||||||||||||
|
|
|
|
|
|
||||||||||||||||||
|
|
|
(пропущено |
|
|
(x x4 ) |
|
). |
|
|
|
||||||||||||
|
|
|
|
(x |
4 |
x |
4 |
) |
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Задача 3. Записать полином Лагранжа L(x) , интерполирующий точки данных, показанные на рис. 4.4 д).
Рис. 4.4.
Решение:
L(x) y1 1(x) y2 2 (x) y3 3 (x) y4 4 (x) . (4.12)
71
Пояснение. Ненулевые ординаты в узлах во вспомогательных задачах с учѐтом масштаба равны соответствующим ординатам полинома на
рис. 4.4 д, а сумма всех ординат в узах xi , i 1, ,4 даѐт все ординаты в узлах xi полинома с рис. 4.4 д. Поэтому сумма вспомогательных
полиномов yi i (x) |
даст полином, интерполирующий узловые точки |
|||||||
данных с рис. 4.4 д. |
|
|
|
|
|
|
|
|
В общем виде интерполяционный полином Лагранжа равен: |
|
|||||||
n |
n |
n |
(x x j ) |
|
||||
L(x) yi |
i (x) yi |
|
|
|
|
|
|
|
(x |
x |
|
) . |
(4.13) |
||||
i 1 |
i 1 |
j 1 ( j i) |
|
|||||
i |
|
j |
|
|
|
Функция, которая вычисляет значение полинома Лагранжа от аргумента x, представлена следующим программным кодом:
Function PLagr( X() As Double, Y() As Double, ByVal N As Long, _ ByVal u As Double) As Double
' N - max индекс в массивах X( ), Y( ) с координатами узлов
Dim P As Double
PLagr = 0#
For i = LBound(X) To N ' ф-ия LBound определяет нижнюю границу индекса
P = Y(i)
For j = LBound(X) To N
If i <> j Then P = P* (u - X(j)) / (X(i) - X(j))
Next j
PLagr = PLagr + P
Next i
End Function
Замечание. Массивы, передаваемые в функции PLagr, могут начинаться как с индекса 1, так и с 0; (в случае нулевого начального индекса число узлов интерполяции равно N+1).
Пример 4.1. Записать в компактном виде |
xi |
0 |
1 |
2 |
3 |
|
полином Лагранжа, интерполирующий за- |
||||||
|
|
|
|
|
||
данные точки: |
yi |
0 |
2 |
1 |
0 |
|
|
|
|
|
|
|
72
Решение примера 4.1 облегчается тем, что ординаты крайних точек равны нулю. С учѐтом этого можно записать:
L(x) 0 1 (x) y2 2 (x) y3 3 (x) 0 4 (x)
0 y2 |
(x x2 ) (x x3 ) (x x4 ) |
|
|
y3 |
(x x1 ) (x x2 ) (x x4 ) |
0 4 (x) |
||||||||||||||||||||
(x x ) (x x ) (x x ) |
|
(x x ) (x x ) (x x ) |
||||||||||||||||||||||||
|
|
1 |
2 |
1 |
3 |
|
1 |
4 |
|
|
3 |
1 |
3 |
|
2 |
2 |
4 |
|
||||||||
0 2 |
(x 0) (x 2) (x 3) |
|
1 |
(x 0) (x 1) (x 3) |
|
0 |
|
|
|
|||||||||||||||||
|
(1 0) (1 2) (1 3) |
(2 0) (2 1) (2 3) |
|
|
|
|||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|||||||||||||||||
0 2 |
(x 0) (x 2) (x 3) |
1 |
(x 0) (x 1) (x 3) |
|
0 |
|
|
|
||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
0 2 |
|
(x 0) (x 2) (x 3) |
1 |
|
(x 0) (x 1) (x 3) |
0 |
|
|
|
|||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
(x 0) (x 3) [2 (x 2) (x 1)] |
|
|
|
x (x 3) (x 3) |
|
|
|
x (x 3)2 |
|
|
|||||||||||||||
|
|
|
|
|
2 |
|
|
|
|
|
|
2 |
|
2 |
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Подстановкой абсцисс заданных точек убеждаемся в выполнении интерполяционных условий.
Представленная выше функция PLagr будет использована в решении примера 4.3, а процедуру решения этого примера можно использовать для контроля только что полученной компактной форму-
лы L(x) x (x 3)2 2 .
4.1.4 Приближение сплайн-функциями
Для гладкой интерполяции часто применяют кубические сплайны. Сплайны второй степени при тех же вычислительных затратах не обеспечивают непрерывность кривизны (вторых производных).
Кубическая сплайн-функция S(x) (или кубический сплайн
S(x) ) на сетке узлов x0 xi 1 xi xn - непрерывная кусочная функция, задаваемая на интервалах [xi 1, xi ] между узлами кубическими полиномами si (x) a0i a1i x a2i x2 a3i x3 , которые гладко (с совпадением первой и вторых производных) стыкуются друг с другом в узловых точках. Такие кубические сплайны имеет разрывы третьей производной (сплайны дефекта 1 =3-2).
73
Сплайн-функции получили своѐ название от гибких деревянных или металлических реек, которые использовались чертѐжниками при проведении плавных линий через заданные точки при разработке (на плазе) чертежей обводов судов или автомобильных кузовов. В настоящее время используют несколько видов математического представления сплайнов. Впервые запись сплайна в «современном» виде
EJ y(x) С2 С1 x M a |
x2 |
Va |
x3 |
P | |
(x x )3 |
Vb | |
(x x )3 |
(4.14) |
|
|
1 |
2 |
|||||
|
|
|
|
|
|
|
||
|
2 |
|
6 |
x x2 |
6 |
x x3 |
6 |
|
как уравнение упругой линии нагруженной балки (рис. 4.5) предло-
жена профессором И.Г. Бубновым в |
|
||
1912 г. в книге «Строительная механика |
|
||
корабля». Прерыватель |
И.Г.Бубнова |
|
|
| |
|
ним |
|
x x указывает, что следующее за |
|
||
i |
|
|
|
выражение добавляется лишь при |
x xi , |
|
|
а иначе оно равно нулю. |
При EJ 1 и |
|
|
большом числе n точек xi |
приложения |
Рис. 4.5. |
сил структура аналитической зависимости y(x) упругой линии такова:
y(x) a |
|
a |
x a |
|
x2 a |
|
x3 in 1 a |
|
| |
(x x )3 |
i3 0 ai xi in 2 ai 2 |
| (x xi )3 (4.15) |
|
0 |
1 |
|
2 |
|
3 |
|
i 2 |
x xi |
i |
|
x xi |
Это и есть уравнение сплайна с прерывателями И.Г. Бубнова. В «современной» записи сплайна И.Г. Бубнова
s(x) 3k 0 ak x3 kn 1 ak 3 (x xk )3 |
|
|
|
(4.16) |
|
просто обозначают прерыватель | (x xi ) |
как (x x |
k |
) |
|
с подстроч- |
x x1 |
|
|
|
||
|
|
|
|
|
ным знаком "+", который имеет точно такой же смысл, что и прерыватель: выражение (x xi ) действует, когда оно положительно, и не действует (равно нулю) при x xi . Это выражается и по другому
| (x xi ) (x xi ) max( 0, x xi ) . |
(4.17) |
x xi |
|
Кубический сплайн однозначно определяется заданием ординат yi в узловых точках xi и условиями на концах: заданием первых производных (как в задании к РГР) или вторых производных (нулевые вторые производные соответствуют распрямлению линии на концах,
74
как в случае реек чертѐжников). Часто используется условие «без узла», когда на двух крайних интервалах сплайн выражается одним полиномом. Это условие даѐт плавное изменение кривизны у краѐв сплайна. Условие периодического повторения обеспечивает в концевых узлах равенство до вторых производных включительно.
При вычислениях значений S(x) вместе с ординатами yi во всех n1=n+1 узлах удобно использовать первые производные pi в узловых точках, которые предварительно вычисляются решением СЛАУ с трѐхдиагональной матрицей:
hi 1 pi 1 2 (hi 1 |
hi ) pi |
hi |
pi 1 3 (hi 1 Ti hi Ti 1 ) , i 1, , n 1, |
(4.18) |
|
где |
hi |
xi |
xi 1 , |
Ti ( yi yi 1 ) hi . |
(4.19) |
Эти уравнения выражают условия равенства вторых производных во внутренних узлах сетки, то есть в местах стыковки сегментов сплайна.
Вычисление недостающих производных pi в узловых точках и расчѐт (по ним и yi ) ординат v S(u) сплайна в точках u выполняется обращением к функции HS (Эримитов сплайн) при Kod≤0, которая описана в пособии [5, библ. №2810]. В задании к данной РГР используют функцию HconS; при Kod=-1 она по заданным p0 , pn (а при иных значениях Kod≤0 по условиям «без узла») вычисляет остальные производные pi в узлах сплайна, заданного на равномерной сетке узлов. Затем она вычисляет ординату S(x) и устанавливается Kod=1 равным номеру очередного обращения для вычисления ординат.
Пример 4.2. Определить для заданных узлов (см. таблицу 4.1 вар. 30)
и вычисленных производных для полинома в |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
xi |
-2 |
-1 |
0 |
1 |
2 |
|
||
концевых узлах ( p0 14 / 3 , p4 |
14 / 3 ) про- |
|
||||||
yi |
2 |
0 |
0 |
0 |
2 |
|
||
изводные кубического интерполяционного |
|
|||||||
|
|
|
|
|
|
|
||
сплайна во внутренних узлах. |
|
|
|
|
|
|
|
|
Решение. Сетка узлов равномерная, и все шаги hi xi xi 1 1. Поэтому уравнения для трѐх внутренних узлов имеют следующий вид:
0 4 p1 1 p2 |
3 ( y1 y0 y2 y1 ) 1 p0 , |
для узла |
i 1, |
|
1 p1 4 p2 |
1 p3 |
3 ( y2 y1 y3 y2 ) , |
для узла i 2 , |
|
1 p2 |
4 p3 |
0 3 ( y3 y2 y4 y3 ) 1 p4 , |
для узла i 3 , |
75
Отметим, что известные концевые производные (даны жирным шрифтом) перенесены в правую часть (в первом и последнем уравнениях). Подстановка ординат и значений производных в эти уравнения
|
4 1 |
0 p1 |
|
|
6 |
|
||||
|
|
|
|
|
|
|
|
|
|
|
даѐт: |
|
1 |
4 |
1 |
p2 |
|
|
0 |
|
систему трѐх уравнений, решая еѐ, |
|
|
0 |
1 |
4 |
|
|
|
6 |
|
|
|
|
p3 |
|
|
|
|
|
p1 |
|
|
1.5 |
|
|
|
|
|
|
|
получаем |
p2 |
|
|
0 |
|
|
|
|
|
1.5 |
|
|
p3 |
|
|
|
По значениям ординат |
y |
i 1 |
, y |
i |
и производных |
p |
, p на концах отрез- |
|||
|
|
|
|
|
i 1 |
i |
|
|||
|
ка [x |
|
, x ] |
|
можно вычислить ординаты |
сегмента |
||||
|
i 1 |
|
i |
|
|
|
|
|
||
|
si (x) |
кубического сплайна S(x) |
либо по формулам: |
|||||||
|
S(x) si (u) yi 1 hi u ( pi 1 |
u (C u D)), (4.20) |
||||||||
|
где |
u (x xi 1 ) / hi , |
|
|
|
|||||
Рис. 4.6 |
D (Ti |
pi 1 ) (Ti pi ) , C D (Ti pi 1 ) , |
(4.21) |
|||||||
либо по всем yi , pi обращением к HconS при зада- |
||||||||||
|
нии Kod>0, указывающего номер обращения к HconS, например, при
Kod=1.
4.1.5 Примеры интерполяции полиномами и сплайнами
Пример 4.3. Интерполировать полиномом и кубическим сплайном n равномерно расположенные на отрезке [0; 3] узловые точки с ординатами, равными округлѐнным до десятых значениями функции
x (x 3)2 2 из предыдущего примера. Проследить на графиках поведение приближающих функций при различных n (3<n<16).
Текст программы решения примера 4.3 включает главную процедуру, функцию PnX для расчета ординат y x (x 3)2 2 и процедуру Graph для вывода на точечную диаграмму узлов интерполяции и приближающих полинома и кубического сплайна. Функция PLagr, представленная выше, вызывается из подключаемого модуля
Math.bas .
Sub INTERPOL( ) |
' Главная процедура решения примера 4.3 |
Const M As Integer = 50 |
' число звеньев на выводимых графиках |
76
Dim Xgr(M) As Double, Ygr1(M) As Double, Ygr2(M) As Double Dim N As Long, Kod As Long, KodErr As Long
Dim N As Long
Dim X(0 To 15) As Double, Y(0 To 15) As Double Dim P(0 To 15) As Double
N = Val(InputBox("Введите N"))
For i = 0 To N
X(i) =0+ 3 * (i / N): Y(i) = Round(PnX(X(i)), 1) ' округление до десятых Cells(i + 2, 1) = X(i): Cells(i + 2, 2) = Y(i)
Next i |
|
|
|
P(0)=0: Kod = -01 |
' Гран.Условия: Sʹʹ(начало)=0, Sʹ(конец)=P(1)=0 |
||
For i = 0 To M |
|
|
|
Xgr(i) = X(0) + (X(N) - X(0)) * (i / M) |
|
|
|
Ygr1(i) = HS(Xgr(i), N, X(0),Y(0), P(0), Kod, KodErr) |
' ф-ию HS (см. в [5]) |
||
If KodErr > 0 Then Stop |
|
' прервать работу из-за ошибки |
|
Cells(i + 2, 3) = Xgr(i): Cells(i + 2, 4) = Ygr1(i) |
|
||
Ygr2(i) = PLagr(X, Y, N, Xgr(i)) |
|
|
|
Cells(i + 2, 5) = Xgr(i): Cells(i + 2, 6) = Ygr2(i) |
|
||
Next i |
|
|
|
' Вызов процедуры построения диаграммы: |
|
||
Call Graph("СПЛАЙН", "Лист1!$C$" & 2 & ":$C$" & (M + 2), _ |
|||
|
|
"Лист1!$D$" & 2 & ":$D$" & (M + 2), _ |
|
"ПОЛИНОМ", "Лист1!$E$" & 2 & ":$E$" & (M + 2), _ |
|
||
|
|
"Лист1!$F$" & 2 & ":$F$" & (M + 2), _ |
|
"УЗЛЫ", "Лист1!$A$" & 2 & ":$A$" & (N + 1),"Лист1!$B$" & 2 & ":$B$" & _ |
|||
(N+1)) |
|
|
|
End Sub |
|
|
|
Function PnX(X As Double) As Double |
' Вычисление ординат полинома |
||
PnX = X * (X - 3) ^ 2 / 2 |
|
|
|
End Function |
|
|
|
Sub Graph(Name1$, Adr1X$, Adr1Y$, Name2$, Adr2X$, Adr2Y$, _ Name3$, Adr3X$, Adr3Y$)
' Процедура вывода на диаграмму трѐх массивов точек (графиков)
77
Name01$ = "=" & """" & Name1$ & """" |
' |
оформление подписи графика 1 |
Name02$ = "=" & """" & Name2$ & """" |
' |
оформление подписи графика 2 |
Name03$ = "=" & """" & Name3$ & """" |
' |
оформление подписи графика 3 |
ActiveSheet.Shapes.AddChart.Select |
|
ActiveChart.SeriesCollection.NewSeries |
|
ActiveChart.SeriesCollection(1).Name = Name01$ |
' "=""СПЛАЙН""" |
ActiveChart.SeriesCollection(1).XValues = Adr1X$ ' "=Лист1!$C$2:$C$32" |
|
ActiveChart.SeriesCollection(1).Values = Adr1Y$ |
' "=Лист1!$D$2:$D$32" |
ActiveChart.SeriesCollection.NewSeries |
|
ActiveChart.SeriesCollection(2).Name = Name02$ |
' "=""ПОЛИНОМ""" |
ActiveChart.SeriesCollection(2).XValues = Adr2X$ |
' "=Лист1!$E$2:$E$32" |
ActiveChart.SeriesCollection(2).Values = Adr2Y$ |
' "=Лист1!$F$2:$F$32" ' ' |
ActiveChart.ChartType = xlXYScatterLinesNoMarkers |
|
ActiveChart.SeriesCollection.NewSeries |
|
ActiveChart.SeriesCollection(3).Name = Name03$ |
' "=""УЗЛЫ""" |
ActiveChart.SeriesCollection(3).XValues = Adr3X$ |
' "=Лист1!$A$2:$A$18" |
ActiveChart.SeriesCollection(3).Values = Adr3Y$ |
' "=Лист1!$B$2:$B$18" |
ActiveChart.SeriesCollection(3).ChartType = xlXYScatter ActiveChart.SeriesCollection(3).Select
Selection.MarkerStyle = 1: Selection.MarkerSize = 6
End Sub
Замечания. В приведѐнной программе сплайн формируется и рассчитывается с помощью универсальной функции HS – она подробно описана в пособии [5, библ. №2810]. В примере задания (см. далее) использована аналогичная функция HconS, вычисляющая сплайн на равномерной сетке узлов.
Процедура Graph - это результат редактирования текста макроса, записанного при «ручном» построении соответствующей точечной диаграммы. Процедура входит в модуль Math.bas (или POL_SPL_mod.bas), который можно импортировать. Вообще говоря, процедуру Graph можно не вызывать, а строить диаграмму вручную.
Для запуска главной процедуры удобно назначить еѐ макрос INTERPOL какой-либо фигуре. После щелчка по этой фигуре в появив-
78
шемся окне указывается число N узлов на отрезке [0; 3], после чего выполняется расчѐт и построение диаграммы. На рисунках 4.7 и 4.8 представлены результаты расчѐта соответственно при N=4 и N=15, из которых видно, что при больших степенях интерполяционный полином имеет значительные осцилляции между узлами на концах отрезка, а сплайн не имеет таких всплесков и при большом числе узлов.
Рис. 4.7. Рис. 4.8.
На поведение участка сплайна существенно влияют лишь близлежащие узлы, а на поведение полинома (степени n) существенно влияет каждый из n+1 определяющих его узлов. (Уместно вспомнить математический анализ, где полином, являясь ограниченным рядом Тейлора, определяется своим значением и значениями своих производных в одной произвольной точке).
По высказанным причинам обычно не используют интерполяционные полиномы высоких степеней, либо у приближающего полинома не используют области вблизи граничных узлов, либо специальным образом выбирают расположение узлов интерполяции, что является сложной задачей.
4.2Задание к расчетно-графической работе №4
4.2.1Пример программы для варианта №30
Вприведѐнной ниже процедуре InterPOL_SPLn для вычисления недостающих производных и ординаты S(u) кубического сплайна на
равномерной сетке с числом N+1 узлов используется импортируе-
мая с модулем Math.bas (или с POL_SPL_mod.bas) функция HconS:
Function HconS(ByVal u As Double, ByVal N1 As Long, _
x() As Double, y() As Double, p() As Double, ByRef Kod As Long) As Double
79
Перед первым обращением к функции следует задать Kod≤0, и также задать: ординату u, число узлов N1=N+1, их абсциссы x и ординаты y; а при Kod=–1 значения производных p(0), p(N) для крайних точек x(0), x(N). Для иных неположительных значений Kod задавать p(0), p(N) не требуется – на краях будут реализованы условие «без узла».
'Предварительно импортируют модуль POL_SPL_mod.bas
Sub InterPOL_SPLn()
Const M As Integer = 31 ' число сегментов, рисуемых на кривых ' Массивы для графиков:
Dim Xgr(0 To M) As Double, Ygr1(0 To M) As Double, Ygr2(0 To M) As Double
Dim N As Long |
' переменная под мах. индекс узла интерполяции |
||
Dim x(0 To 22) |
As Double, y(0 To 22) As Double ' массивы под узлы |
||
Dim p(0 To 22) |
As Double |
' массив под производные |
|
Dim Kod As Long |
|
' код, управляющий формированием сплайна |
|
N = Val(InputBox("Введите N")) |
' ввод наибольшего значения индекса |
||
For i = 0 To N |
|
' расчѐт узловых точек: |
x(i) = -2 + 4 * (i / N): y(i) = Round(POL(x(i)), 1) Cells(i + 2, 1) = x(i): Cells(i + 2, 2) = y(i)
Next i
' Назначение параметров сплайна с производными на краях:
Kod = -1: p(0) = -14 / 3: p(N) = 14 / 3
' Назначение сплайна с условиями «без узла» по краям:
' Kod = -2 ' задать для N>4
' Вычисление параметров сплайна и полинома
For i = 0 To M
Xgr(i) = x(0) + (x(N) - x(0)) * (i / M) Cells(i + 2, 3) = Xgr(i)
Cells(i + 2, 4) = HconS(Xgr(i), N+1, x, y, p, Kod) Ygr2(i) = PLagr(x, y, N, Xgr(i))
Cells(i + 2, 5) = Xgr(i): Cells(i + 2, 6) = Ygr2(i)
Next i
Call Graph("СПЛАЙН", "Лист1!$C$" & 2 & ":$C$" & (M + 2), _
"Лист1!$D$" & 2 & ":$D$" & (M + 2), _
80