close

Вход

Забыли?

вход по аккаунту

?

FreeBASIC7

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. В главе 7 "Аппроксимация. Сплайны" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также преподавателей.
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 7 "Аппроксимация. Сплайны" –
продолжение пособия, к
оторое будет интересно для учащихся школ,
студентов институтов, а также преподавателей.
Фрагмент 7. Аппроксимация. Сплайны.
Начнем, пожалуй, с широко цитируемого, но от этого не становящегося устаревшим знаменитого высказывания Эдсгера Дейкстры:
"И
нформатика не более наука о компьютерах,
чем астрономия —
наука о телескопах".
К сожалению, под рукой не оказалось первоисточника, но вот перевод этого высказывания, приведенный факультетом Информатики Технического Университета Вены:
"В информатике речь ид
ёт настолько мало о компьютерах как в Астрономии о Телескопах".
Конечно, смысл совпадает, но приятно иметь под рукой и "авторский текст" (об этом еще придется вспомнить ниже)...
2
0.
Вместо введения
Новогодние праздники, перемежаемые Рождеством Христовым, как это исстари принято у инженеров, в ночь с 30 декабря на 14 января
–
благополучно завершились! Почему
-
то вспомнилась математическая задачка
писателя
-
сатирика Аркадия Михайловича Арканова: На моём дне рождения папа принял четыреста. Мама -
пятьсот. Елена
Георгиевна семьсот пятьдесят. Дедушка -
чекушку. Ему больше нельзя, так как он старенький и у него язва. Вопрос: Сколько принял товарищ Ерохин, если папа вырубился в час ночи, мама -
в два часа ночи, Елена Георгиевна вообще не вырубалась, а товарищ Ерохин
подрался с дедушкой? А за этими воспоминаниями потянулись цветные картинки:
"Пир у Симона фарисея" (около 1618 года)
Картина написана Питером Паулем Рубенсом (Pieter Paul Rubens, 1577 —
1640) —
фламандским (южнонидерландским) живописцем. Сцена застол
ья... смысл таков -
во время пира у фарисея Симона к Христу подошла грешница и полила его ноги душистым маслом. Симон осудил невозмутимость Христа, сказав: "Грешница посмела осквернить Христа прикосновением". Христос ответил: "Прощаются грехи ее многие за то, что она возлюбила много"...
Очень прошу –
не стесняйтесь, пишите, вы же видите, я "пытаюсь соответствовать" всем вашим запросам! Конечно, пока только в рамках основной цели "Фрагментов" –
записать небольшие, но полезные алгоритмы
(из минимального джен
тльменского набора) в виде доступных для всех текстов на языке FreeBASIC. Надеюсь, настанет время, когда сможем пофантазировать над собственными алгоритмами, которых еще не знает инженерная братия... Вот таков мой краткий ответ на заданный читателем вопрос
. Но он не вполне полон, нет философской подоплеки. А поэтому продолжу.
3
Аркадий Натанович и
Борис Натанович Стругацкие
Есть одна из важнейших каждодневных проблем, которые приходится решать каждому взрослому человеку (возможно, не только инженеру) это пр
облема выбора. Вы читали повесть братьев Стругацких (Аркадий Натанович; 1925
—
1991 и Борис Натанович; 1933
—
2012) "За миллиард лет до конца света"
(1976)? Тема выбора основная в этой повести, герои которой поставлены перед жестокой необходимостью выбирать ме
жду возможностью творить под угрозой смерти, либо отказаться от своих убеждений ради спокойной жизни. Действие происходит в 70
-
x годах XX
-
го века в СССР, в Ленинграде. Целая группа людей неожиданно осознала, что какая
-
то загадочная, но могучая и очень изби
рательная сила останавливает их разработки в самых разных областях науки от биологии до математической лингвистики
. Стоит им заняться определенной научной работой, тут же начинаются странности.
Один из них, математик Вечеровский предполагает:
"Если бы сущ
ествовал только закон неубывания энтропии
, структурность мироздания исчезла бы, воцарился бы хаос. Но, с другой стороны, если бы существовал или хотя бы возобладал только непрерывно совершенствующийся и всемогущий разум, заданная гомеостазисом структура ми
роздания тоже нарушилась бы. Это, конечно, не означало бы, что мироздание стало бы хуже или лучше, оно бы просто стало другим, ибо у непрерывно развивающегося разума может быть только одна цель: изменение природы Природы.
Поэтому сама суть Гомеостазиса Мир
оздания состоит в поддержании равновесия между возрастанием энтропии и развитием разума. Поэтому нет и не может быть сверхцивилизаций, ибо под сверхцивилизацией мы подразумеваем именно разум, развившийся до такой степени, что он уже преодолевает закон неуб
ывания энтропии в космических масштабах. И то, что происходит сейчас с нами, есть не что иное, как первые реакции Гомеостатического Мироздания на угрозу превращения человечества в сверхцивилизацию. Мироздание защищается
"...
В прошлом выпуске пожаловался, что хотел, но не смог написать очередной "Фрагмент" про алгоритмы "аппроксимации". Не смог, но постоянно помнил! Это, действительно очень важная тема.
1.
Из записок сумасшедшего
Конечно, эта глава не является обязательной. Даже не стану агитировать бой
ких читателей и милых читательниц за ознакомление с ней. Просто считаю своим долгом напомнить, что в каждом продукте, который мы потребляем сегодня, есть капля пота и крови ученого и инженера жившего до нас!
Об этом стараются не вспоминать люди, стремящиес
я заполучить "венки первооткрывателей" и просто обыватели, 4
которые таким образом уходят от осознания малости своего вклада во Всеобщую историю Человечества... Поэтому потерпите –
будут ссылки на достойных, на мой взгляд (может не совпадать с "мнением админ
истрации"), людей...
И так, уже упомянут Эдсгер Вибе Дейкстра (Edsger Wybe Dijkstra; 1930 —
2002) —
нидерландский ученый, идеи которого оказали влияние на развитие компьютерной индустрии.
Эдсгер Вибе Дейкстра
Родился в семье ученых: отец —
химик, мать —
математик
. Такое или подобное сочетание (опять же, на мой взгляд) обладает потенцией саморазвития (ниже мы коснемся системы "предиктор
-
корректор"). После школы поступил на факультет теоретической физики Лейденского университета. В 1951 увлекся программиров
анием, поступил на трехнедельные компьютерные курсы в Кембридже, с 1952 работал программистом в Математическом центре Амстердама под руководством профессора Адриана ван Вейнгаардена (1916 –
1987), впоследствии —
автора одного из способов формального описан
ия грамматики формальных языков —
так называемых двухуровневых грамматик Ван Вейнгаардена. Уже в 1952 принял решение окончательно специализироваться на программировании, но курс теоретической физики закончил. В 1956 принял участие в разработке ЭВМ "X1"
. В "X1" впервые были использованы прерывания, позволяющие прервать ввод и вывод ради синхронизации его с вычислениями. "X1" был одним из лучших компьютеров, как по элегантности внутреннего дизайна, так и по скорости выполнения вычислений, но "правила рынка" н
е позволили ему конкурировать с "признанными авторитетами" (теперь
-
то мы догадываемся, кто и как навязывает "авторитеты"
). Концептуально "X1" был предметом диссертации Эдсгера Дейкстры и первой платформой для реализации его полноценного компилятора языка А
лгол 60.
2.
Коротко о методе аппроксимации
Аппроксимация, или приближение —
научный метод, состоящий в замене одних объектов другими, в том или ином смысле близкими к исходным, но более простыми
. Метод позволяет исследовать числовые характеристики и ка
чественные свойства объекта, сводя задачу к изучению более простых или более удобных объектов (например, таких, характеристики которых легко вычисляются, или свойства которых уже известны). В теории чисел изучаются диофантовы (Диофант Александрийский уже б
ыл упомянут в 5 и 6
-
ом выпусках "Фрагментов") приближения, в частности, приближения иррациональных чисел рациональными. В геометрии рассматриваются аппроксимация кривых отрезками прямых. Некоторые разделы математики целиком посвящены аппроксимации, наприме
р, теория приближения функций, численные методы анализа.
Термин "аппроксимация" употребляется в переносном смысле в философии как метод приближения, указание на приблизительный, неокончательный характер. В 5
таком смысле термин "аппроксимация" активно употр
еблялся датским философом и теологом Сёреном Кьеркегором (Soren Aabye Kierkegaard; 1813
—
1855) в трактате "Заключительное ненаучное послесловие к философским крохам" (1846). О Кьеркегоре в статье "О фанатизме, ортодоксии и истине" (1935) русский религиозный
и политический философ, представитель экзистенциализма Николай Александрович Бердяев (1874 —
1948) пишет:
Николай Александрович Бердяев
Тема о фанатизме, связанная с приверженностью к ортодоксальным учениям, очень актуальна... Пафос ортодоксальной докт
рины, которая оказывается полезной для борьбы и для организации, ведет к полной потери интереса к мысли и к идеям, к познанию
, к интеллектуальной культуре, и сравнение со средневековьем очень неблагоприятно для нашего времени. Никакого идейного творчества при этом не обнаруживается. В этом отношении наша нетерпимая эпоха поразительно бездарна и убога, в ней творческая мысль замирает, она паразитарно питается предшествующими эпохами. Мыслители наиболее влиятельные в современной Европе, -
как Маркс, Ницше, Ки
ркегардт, -
принадлежат тому ХIХ веку, против которого сейчас происходит реакция. Единственная область, в которой обнаруживается головокружительное творчество, есть область технических открытий. Мы живем под знаком социальности, и в этой области происходит
много положительного, но никаких социальных идей, социальных теорий сейчас не создается
, все они принадлежат ХIХ веку. Марксизм, прудонизм, синдикализм, даже расизм, -
всё порождение мысли ХIХ века. Главное преимущество нынешнего века в том, что он более обращен к реальностям, разоблачает реальности. Но, разоблачая старые идолы, новый век создает новые идолы.
Процессы аппроксимации приобрели особо актуальное значение в связи с ростом числа исследований сложных систем. Когда системы становятся сложными, го
ворит У. Эшби, то их "
теория практически заключается в том, чтобы найти пути упрощения
" (1964).
Уильям Росс Эшби
(William Ross Ashby; 1903 —
1972) —
английский психиатр, специалист по кибернетике, пионер в исследовании сложных систем. Эшби принадлежит изо
бретение гомеостата
(1948), введение понятия самоорганизации
. Он сформулировал закон о требуемом разнообразии
, названный его именем (закон Эшби): "управление может быть обеспечено только в том случае, если разнообразие средств управляющего (в данном случае
всей системы управления) по крайней мере не меньше, чем разнообразие управляемой им ситуации" (1959). Иногда этот закон "упрощают" и формулируют так: "что управляющий должен иметь большее разнообразие". Это неверно. Регулятор может 6
иметь разнообразие мень
шее, чем управляемая система и это почти всегда имеет место быть на практике. Просто в этом случае регулятор не имеет возможности полностью управлять ситуацией.
Близким понятием к "аппроксимации" является "
абстракция
" (abstractio —
отвлечение) —
отвлечени
е в процессе познания от несущественных сторон, свойств, связей объекта (предмета или явления) с целью выделения их существенных, закономерных признаков; абстрагирование; теоретическое обобщение как результат такого отвлечения. В философии и логике абстраг
ирование трактуется как способ поэтапного продуцирования понятий, которые образуют всё более общие модели —
иерархию абстракций. Наиболее развитой системой абстракций обладает математика. В зависимости от целей и задач, можно рассуждать об одном и том же о
бъекте на разных уровнях абстракции.
3.
Феномен (явление) Рунге
В качестве оселка
(Touchstone, Schleifstein) т.е. средства для апробирования программ ("Деньги и богатство являются оселком, на котором испытывается человек" -
см. Большой толковый словарь
русского языка. Сост. С.А.Кузнецов. 1998, 1536 с.), связанных с аппроксимацией, будем использовать функцию Рунге
. Феноменом Рунге называют эффект нежелательных осцилляций (колебаний), возникающий при интерполяции функции полиномами высоких степеней. Эффек
т был открыт Карлом Рунге при изучении ошибок полиномиальной интерполяции для приближения некоторых функций. Рунге (1901) исследовал интерполяцию полиномами функции (плотность вероятности распределения Коши) следующего вида:
f(x) = 1/(1+25*x^2)
Результат а
ппроксимации ее полиномом представлен на рисунке:
Графики функции Рунге и интерполяционного полинома
7
Карл Давид Тольме Рунге
(Carl David Tolme Runge; 1856 -
1927) —
немецкий математик, физик и спектроскопист. Учился в Берлинском университете у Карла Вей
ерштрасса, в 1880 получил степень доктора философии по математике, с 1886 —
профессор математики в Ганноверском университете. В 1904 по инициативе Феликса Клейна приглашён в Университет Георга
-
Августа в Гёттингене и возглавил вновь открытую кафедру приклад
ной математики. Считается исторически первым немецким математиком по этой дисциплине. В Гёттингене, совместно с М.Куттой разработал методы численного интегрирования
систем обыкновенных дифференциальных уравнений —
методы Рунге
—
Кутты. Исследовал поведение п
олиномиальной интерполяции
при повышении степени полиномов —
Феномен Рунге. Из
-
за стремления уменьшить ошибку феномена Рунге и появилась идея использования интерполяции сплайнами, как метода обработки данных.
Для подтверждения плодотворности этой идеи пре
длагается немного поэкспериментировать с программой аппроксимации функции Рунге, текст которой приведен ниже. В ней используются две процедуры DECOMP() и SOLVE()
, которые были приведены в статье "Фрагмент 6. Линейная алгебра" -
DECOMP выполняет ту часть га
уссова исключения, которая зависит лишь от матрицы. Она сохраняет множители и информацию о ведущих элементах. Подпрограмма SOLVE использует эти результаты, чтобы получить решение для произвольной правой части.
Для наглядности воспроизведения результатов на
"хвост" программы "посажен" блок графического вывода. Конечно, программу можно сделать более красивой и более компактной, но хотелось сохранить "девственность" используемых элементов...
' P R O G R A M "Appr_Runge"
' 12.12.2012
'
-------------------------------------------------------------
' Попытки аппроксимации функции Рунге F(x) = 1/(1+25*X^2)
'
-------------------------------------------------------------
'
#Lang "FB" ' режим FreeBASIC
-
сов
местимости
' Описание используемых подпрограмм (процедур)
Declare Sub DECOMP(NDIM As Integer, N As Integer,_
A() As Single, ByRef COND As Single, IPVT() As Integer,_
WORK() As Single)
Declare Sub SOLVE(NDIM As Integer, N As Integer,_
A() As Single, WORK() As Single, IPVT() As Integer)
' Подпрограмма вычисление значений функции Рунге
Function Runge(X As Single) As Single
' X -
аргумент функции (X <= 35.0);
' FR -
значение
функции
.
Dim As Single FR
FR = 1/(1+25*X^2)
Runge = FR
End Fun
ction
8
' Основной модуль тестовой программы
Dim As Integer NDIM
NDIM = 10 ' резервируемая размерность массива
Dim As Single A(NDIM, NDIM), B(NDIM), WORK(NDIM), COND, CONDPI
Dim As Integer IPVT(NDIM), I, J, N, NVar
N = 3 ' текущий порядо
к системы
Dim As Single X(N), Y(N), XX, YY, ZZ Data 0.000, 0.125, 0.250 ' значения аргументов вар. 1
Data 0.375, 0.500, 0.625 ' значения аргументов вар. 2
Data 0.750, 0.875, 1.000 ' значения аргументов вар. 3
For NVar = 1 To 3 ' рассматриваются три вариа
нта
Screen 0 ' вывод результатов на консоль
Print: Print Using " NVar = ##"; NVar
For I = 1 To N: Read X(I):
Y(I) = Runge(X(I)): Next I
For I = 1 To N
Print Using " X = ##.###### Y = ##.######"; X(I); Y(I)
Next I
' Формирование матрицы коэффицие
нтов уравнений
For I = 1 To N
A(I, 1) = X(I)*X(I)
A(I, 2) = X(I)
A(I, 3) = 1
Next I
Print " Matrix A(I,J):"
For I = 1 To N
For J = 1 To N
Print Using " ###.###### "; A(I, J);
Next J
Print
Next I
' Вызов
процедуры
DECOMP
DECOMP (NDIM, N, A(), COND, IPVT(), WORK())
Print Using " COND = ###.###### ";COND
CONDPI = COND + 1
If (CONDPI = COND) Then
Print " Matrix is Singular!"
GoTo St
EndIf
' Формирование вектора правой части уравнений
For I = 1 To N: B(I) = Y(I): Next I
Print " Vector B(I):"
For I = 1 To N
Print Using " ###.###### "; B(I);
Next I
Print
' Вызов процедуры SOLVE
SOLVE (NDIM, N, A(), B(), IPVT())
Print " Vector X(I):"
For I = 1 To N
Print Using " ###.###### "; B(I);
Next I
Print
9
For I = 1 To N
ZZ = B(1)*X(I)*X(I)+B(2)*X(I)+B(3)
Print Using "##.
###### ##.###### ##.######"; X(I); Y(I); ZZ
Next I
St:
' Вывод графиков функции Рунге и вариантов парабол
Screen 12 ' вывод на граф. экран 640x480
View (1, 1)
-
(478, 478), , 15 ' квадратное физическое окно
' самое большое для данного ре
жима с белым кантиком...
Window (0, 0)
-
(1, 1) ' логическое окно вывода
For XX = 0 To 1 Step 0.125
Line (0, XX)
-
(1, XX), 6, , &HAAAA ' горизонтали
сетки
Line (XX, 1)
-
(XX, 0), 6, , &HAAAA ' вертикали сетки
Next XX
Line (0, 0)
-
(0, 0), 7, , &HAAAA '
линия
оси
X
Line (0, 1)
-
(0, 0), 7, , &HAAAA ' линия
оси
Y
'
For XX = 0 To 1 Step 0.001 ' перебор
аргумента
YY = Runge(XX) ' вычисление функции
PSet(XX, YY), 13 ' вывод розовой точки графика
ZZ = B(1)*XX*XX+B(2)*XX+B(3
)
PSet(XX, ZZ), 11 ' вывод салатовой точки графика
Next XX ' следующая точка
Locate 2, 2: Print "Y" ' метка
оси
Y
Locate 29,61: Print "X" ' метка
оси
X
Sleep
Next NVar
Sleep
Ниже приведены результат
ы выполнения программы для трех групп исходных троек узлов интерполяции, расположенных в I квадранте (учитывая симметричность функции).
Вариант NVar = 1
X = 0.000000 Y = 1.000000
X = 0.125000 Y = 0.719101
X = 0.250000 Y = 0.390244
Matrix A(I
,J):
0.000000 0.000000 1.000000
0.015625 0.125000 1.000000
0.062500 0.250000 1.000000
COND = 181.136642
Vector B(I):
1.000000 0.719101 0.390244
Vector X(I):
-
1.534670 -
2.055357 1.000000
0.000000 1.000000 1
.000000
0.125000 0.719101 0.719101
0.250000 0.390244 0.390244
10
Вариант NVar = 2
X = 0.375000 Y = 0.221453
X = 0.500000 Y = 0.137931
X = 0.625000 Y = 0.092888
Matrix A(I,J):
0.140625 0.375000 1.000000
0.250000 0.500000 1.000000
0.390625 0.625000 1.000000
COND = 322.881439 Vec
tor B(I):
0.221453 0.137931 0.092888
Vector X(I):
1.231343 -
1.745603 0.702897
0.375000 0.221453 0.221453
0.500000 0.137931 0.137931
0.625000 0.092888 0.092888
Вариант NVar = 3
X = 0.750000 Y = 0.066390
X = 0.875000 Y = 0.049651
X = 1.000000 Y = 0.038462
Matrix A(I,J):
0.562500 0.750000 1.000000
0.765625 0.875000 1.000000
1.000000 1.000000 1.000000
COND = 505.129333 Vector B(I):
0.066390 0.049651 0.038462
Vector X(I):
0.177594 -
0.422503 0.283371
0.750000 0.066390 0.066390
0.875000 0.049651 0.049651
1.000000 0.038462 0.038462
Полученные результаты только с очень большой натяжкой можно признать удовлетворительными. Другого и не могло быть -
сказывается выбор элементов для аппроксимации (параболы) и выбор расположения узлов (просим прощения у П.Л.Чебышёва)!
3.1.
Варианты аппро
ксимации
Еще два небольших штриха к главе. Следует подчеркнуть, что, как правило, продуктом экспериментальной работы или теоретических исследований являются большие массивы численной информации:
одномерные массивы (временные сигналы),
двумерные массивы (п
лоские изображения),
многомерные массивы (результаты промышленных экспериментов).
11
Такие данные подлежат обработке, целью которой может быть:
удаление шума, сглаживание данных;
выявление скрытых закономерностей;
сжатие данных для их хранения;
сжатие данны
х для передачи по каналам связи.
Аппроксимация данных может состоять как в интерполяции, так и в экстраполяции. При этом в вычислительной математике под интерполяцией (интерполирование, interpolation
–
внутри, между) понимается способ нахождения промежуто
чных значений величины по имеющемуся дискретному набору известных значений. Как правило, на основании этих наборов требуется построить некоторую функцию, на график которой могли бы с высокой точностью попадать другие, например, получаемые в эксперименте зн
ачения. Интерполирующая функция часто называется "интерполянт". Интерполяционные ряды вошли в математику в основном благодаря И.Ньютону. В XVIII веке бесконечными интерполяционными рядами как инструментом математического анализа широко пользовались Эйлер, Лагранж и Лаплас, в XIX веке —
Гаусс, Абель и Коши. В конце XIX века обобщение задач интерполирования послужило одним из источников проблемы моментов в работах Чебышева, Стилтьеса и Маркова. Экстраполяция (экстраполирование, extrapolation
—
вне, снаружи) —
особый способ аппроксимации, при котором находятся значения функции вне заданного интервала.
Комбинация указанных выше способов аппроксимации составляет основу так называемой схемы "предиктор
-
корректор"
(метод прогноза и коррекции, предсказывающе
-
исправл
яющий метод). В вычислительной математике это семейство алгоритмов численного решения различных задач, которые состоят из двух шагов. На первом шаге (предиктор) вычисляется грубое приближение требуемой величины. На втором шаге при помощи иного метода прибл
ижение уточняется (корректируется). Эти методы (алгоритмы) являются одними из наиболее популярных многошаговых методов. Их неоспоримым достоинством является высокая точность. Недостаток заключается в том, что они не являются самостартующими, т.е. требуют а
ккуратного выбора исходных точек...
3.2.
Теория автоматического управления
Теория автоматического управления (ТАУ) —
научная дисциплина, изучающая процессы автоматического управления объектами разной физической природы
. При помощи математических средс
тв и методов выявляются общие свойства систем автоматического управления и разрабатываются рекомендации по их проектированию. Многие положения Теории Автоматического Управления содержатся в Общей теории регуляторов, которая была разработана русским ученым Иваном Алексеевичем Вышнеградским (1831 —
1895), почетным членом Петербургской АН (1888), в 1887
—
1892 —
министр финансов России. Вначале теория создавалась для изучения статики и динамики процессов автоматического управления техническими объектами -
произв
одственными, энергетическими, транспортными, а в последнее время делаются попытки (не всегда успешные) использования ее выводов и результатов для изучения динамических свойств систем управления экономических, организационных, биологических объектов.
12
Основ
ополагающими трудами Вышнеградского являются: "Об общей теории регуляторов" (1876), "О регуляторах непрямого действия" (1877). В этих работах можно найти истоки современных инженерных методов исследования устойчивости и качества регулирования, "системный п
одход" к проблеме. Так система автоматического управления рассматривается как единая динамическая система, состоящая из двух частей —
"объекта" (управляемая часть системы, как правило, неизменяемая) и "регулятора" (управляющая часть системы, создаваемая дл
я объекта).
Регуляторы следят за изменением некоторых параметров объекта управления и реагируют на их изменение с помощью некоторых алгоритмов управления в соответствии с заданным критерием управления. По существу подход был тесно связан с теорией устойчив
ости "в малом" А.М.Ляпунова, но имел ярко выраженную инженерную направленность.
Регуляторы в подавляющем большинстве работают по принципу отрицательной обратной связи с целью компенсировать внешние возмущения, действующие на объект управления и отработать
заданный извне или заложенный в системе закон управления. Хорошим тоном сегодня признается схема управления, в которой управляющий сигнал вырабатывается не только на основе информации о текущем состоянии системы (объекта), но и на основе прогноза ее дальн
ейшего поведения. Такая схема называется "предиктор
-
корректор" (или "предсказатель
-
поправщик") и позволяет обеспечить наиболее высокое качество управления, поскольку часть контуров циркуляции информации замкнута не через "свершившееся прошлое", а через "пр
огнозируемое будущее". Это обстоятельство позволяет минимизировать запаздывание управления относительно возмущающего воздействия, а при необходимости перейти к упреждающему управлению
, при котором управляющее воздействие упреждает причину, вынуждающую к уп
равлению.
По отношению к социальным системам управление по схеме предиктор
-
корректор, как явствует из истории, осуществлялось уже в древности. Так высшее жречество древнего Египта звалось "иерофантами" ("знающий будущее" —
у древних греков старший пожизне
нный жрец при Элевсинских таинствах), что означало их умение читать судьбу (т.е. матрицу возможных состояний), предвидеть будущее. Последнее есть основа управления, поскольку управлять -
это значит на основе знания возможных состояний приводить систему (в данном случае, —
общество) к избранному определённому варианту из множества возможных
. Естественно, что избрание варианта обусловлено истинной нравственностью и произволом тех, кто поднялся до предвидения и управления на его основе. Интересное з
амечание сделал русский историк Василий Осипович Ключевский (1841 -
1911): "Великоросс часто думает надвое и это кажется двоедушием. Он всегда идёт к прямой цели, но идёт, оглядываясь по сторонам, и потому походка его кажется уклончивой колеблющейся. Ведь лбом стену не прошибёшь и только вороны летают прямо"...
4.
Что такое "сплайн"
Существует большое количество конструкций, которые называют сплайнами. Поэтому внесем некоторую определенность в это многообразие, имея целью выделить те признаки, которые позволят выбрать сплайны годные для конкретной прикладной задачи. В первую очередь необходимо подчеркнуть, что для сплайнов характерны следующие признаки:
13
сплайн состоит из фрагментов, т.е. функций одного класса, которые отличаются только своими параметра
ми,
на соседние фрагменты в точках стыковки накладываются определенные условия, которые сводятся к непрерывности значений и некоторых первых производных.
Основная идея применения сплайнов состоит в следующем. Интервал, на котором восстанавливают функцию, разбивают на фрагменты (подинтервалы), на каждом из которых функцию задают, например, полиномом достаточно низкой степени и обеспечивают непрерывность кривой в точках "склейки" путем приравнивания значений полиномов на границах подинтервалов. При этом важн
ым условием является также непрерывность нескольких производных. Иными словами, "сплайн" (spline –
рейка, гибкая линейка) -
функция, область определения которой разбита на конечное число отрезков, на каждом из которых сплайн совпадает с некоторым алгебраич
еским полиномом. Максимальная степень из использованных полиномов называется степенью сплайна. Разность между степенью сплайна и получившейся гладкостью называется дефектом сплайна... Хотя сплайны и склеены из кусков алгебраических многочленов, они имеют к
онечную гладкость и внутренняя природа их принципиально иная, чем у полиномов, являющихся аналитическими функциями.
4.1.
Вычисление значений интегралов Френеля
Для "разминки" решил посмотреть популярную в прошлом книгу
: Дьяконов В.П. Справочник по алго
ритмам и программам на языке бейсик для персональных ЭВМ. Справочник. М., 1989, 240 с., ил.
В этой книге в Приложении на странице 231 есть раздел:
§ П5.11. Сплайн
-
аппроксимация, интерполяция и экстраполяция.
Полиномиальная интерполяция (аппроксимация) не о
беспечивает непрерывность производных функции у(х) и может давать значительные погрешности в промежутках между узлами функции. Кроме того, она плохо приспособлена для экстраполяции и, как правило, не обеспечивает правильное асимптотическое поведение функци
и у(х) при изменении аргумента х за пределами интервала интерполяции. Нередко с увеличением числа узлов погрешность такой интерполяции не только не уменьшается, но и начинает расти. От этих недостатков свободна аппроксимация и интерполяция с помощью сплайн
-
функций...
Далее в книге приведен текст программы на языке BASIC и контрольный пример, демонстрирующий сплайн
-
аппроксимацию для интеграла Френеля С(х), заданного таблицей из N = 7 пар чисел xi и yi. К слову, в наш "информационный век" оказалось не так пр
осто найти таблицы значений этой функции. Поэтому, для желающих поэкспериментировать с "синус"
-
и "косинус"
-
функциями Френеля, ниже приведен фрагмент таблицы их значений
:
x S(x) C(x)
0.2 0.0041876092 0.1999210576
0.3 0.0141169980 0.29940097
61
0.4 0.0333594327 0.3974807592
0.5 0.0647324329 0.4923442259
0.6 0.1105402074 0.5810954470
0.7 0.1721364579 0.6596523519
14
0.8 0.2493413930 0.7228441700
0.9 0.3397763440 0.7648230210
1 0.438259147000 0.7798934000
1.1 0.5364979111 0.7638066661
1.2 0.6234009
190 0.7154377230
1.3 0.6863332855 0.6385504547
1.4 0.7135250774 0.5430957840
1.5 0.6975049601 0.4452611760
1.6 0.6388876835 0.3654616834
1.7 0.5491959403 0.3238268760
1.8 0.4509387693 0.3336329272
1.9 0.3733473178 0.3944705350
2 0.343415678400 0.4882534061
Тут сработало чисто "инженерное" любопытство (что за зверь этот интеграл Френеля?), которое обошлось автору в 3 дня специальных исследований!
Огюстен Жан Френель
(Augustin
-
Jean Fresnel; 1788 —
1827), французский физик, один из создателей волновой теори
и света.
Основные работы Френеля посвящены физической оптике. Физику изучал самостоятельно после ознакомления с работами Э.Малюса и начал абсолютно самостоятельно проводить небезуспешные эксперименты по оптике. В 1815 году переоткрыл принцип интерференции
, проделав по сравнению с Томасом Юнгом несколько новых опытов (в частности опыт с "бизеркалами Френеля"). В 1816 году дополнил принцип Гюйгенса, введя представление о когерентной интерференции элементарных волн, излучаемых вторичными источниками (принцип Гюйгенса —
Френеля). Исходя из этого принципа, в 1818 году разработал теорию дифракции света, на основе которой предложил метод расчёта дифракционной картины, основанный на разбиении фронта волны на зоны (так называемые зоны Френеля). С помощью этого метод
а рассмотрел задачу о дифракции света на краю полуэкрана и круглого отверстия. В 1821 независимо от Юнга доказал поперечность световых волн. В 1823 установил законы изменения поляризации света при его отражении и преломлении (формулы Френеля). Изобрел неск
олько новых интерференционных приборов (зеркала Френеля, бипризма Френеля, линза Френеля). В 1823 Френель был избран членом Парижской АН, в 1825 стал членом Лондонского королевского общества. Его имя внесено в список величайших учёных Франции
, помещённый н
а первом этаже Эйфелевой башни.
Для предварительного знакомства с особенностями интегралов Френеля, можно воспользоваться фундаментальным справочником: "Справочник по специальным функциям, с формулами, графиками и математическими таблицами". Абрамовиц М.
, Стиган И., 1979, 832 с., ил. (Abramowitz, Stegun. Handbook)
На странице 123 находится раздел: 7.3. Интегралы Френеля.
15
Фрагмент справочника
Обратите внимание, что численные значения интегралов Френеля определяются тремя различными формулами
. Конечно
, наибольшее хождение имеют формулы (7.3.1) и (7.3.2), но среди физиков популярны и формулы (7.3.3) и (7.3.4) в варианте C2(x) и S2(x).
Ниже приведен текст ординарно переработанной подпрограммы SF45R из:
Библиотека численного анализа НИВЦ МГУ
Вычисление значений интегралов Френеля S(x) и C(x)
http://num
-
anal.srcc.msu.ru/lib_na/cat/sf/sf45r.htm
Назначение подпрограммы: Вычисление значений интегралов Френеля S(x) и C(x).
Параметры вызова:
X (X1) -
заданное значение аргумента x (тип: вещественный);
C, S -
в
ещественные переменные, которым в результате работы подпрограммы присваиваются значения C(x) и S(x) соответственно. При вводе в подпрограмму SF45R для отрицательных значений аргумента x < 0 переменным C и S присваиваются значения C(|x|) и S(|x|) соответств
енно.
В примере на сайте приводятся результаты вычисления функций для x = 13:
C = 0.5425104114
S = 0.3982677211
16
' P R O G R A M "IntegrFren_SF45R"
' 10.10.2012
'
--------------------------------
-----------------------------
' Вычисление значений интегралов Френеля S(x) и C(x)
'
-------------------------------------------------------------
'
#Lang "FB" ' режим FreeBASIC
-
совместимости
Sub SF45R(X As Single, ByRef C As Single, ByRef S As Single)
Dim As Single A(52), RK(13), RL(13), Pi, S2Pi
Dim As Integer K, L, I
Dim As Single Z, H, Y, F, D, E, B, R, T
'
A(1) = 0.1E
-
10: A(2) = -
0.366E
-
9: A(3) = 0.10898E
-
7
A(4) = -
0.267681E
-
6: A(5) = 0.527608E
-
5: A(6) = -
0.8105684E
-
4
A(7) = 0.9339901E
-
3: A(
8) = -
0.7651297E
-
2: A(9) = 0.0411409E0
A(10) =
-
0.1271339E0: A(11) = 0.1743607E0: A(12) = -
0.0808111E0
A(13) = 0.5479103E0: A(14) = 0.4E
-
11: A(15) = -
0.128E
-
9
A(16) = 0.4206E
-
8: A(17) = -
0.11507E
-
6: A(18) = 0.2562196E
-
5:
A(19) = -
0.4532192E
-
4: A(20) = 0.617
4202E
-
3: A(21) = -
0.6220184E
-
2
A(22) = 0.0438681E0: A(23) = -
0.2007174E0: A(24) = 0.5386666E0
A(25) = -
0.7996168E0: A(26) = 1.053859E0: A(27) = 0.1E
-
11
A(28) = -
0.4E
-
11: A(29) = 0.14E
-
10: A(30) = -
0.54E
-
10
A(31) = 0.239E
-
9: A(32) = -
0.1176E
-
8: A(33) = 0.65
45E
-
8
A(34) = -
0.42829E
-
7: A(35) = 0.347441E
-
6: A(36) =
-
0.3810219E
-
5
A(37) = 0.6627508E
-
4: A(38) = -
0.2617529E
-
2: A(39)=0.9945488E0
A(40) =0.2E
-
11: A(41)=
-
0.6E
-
11: A(42)=0.18E
-
10:A(43)=
-
0.72E
-
10
A(44) = 0.298E
-
9: A(45) = -
0.1346E
-
8: A(46) = 0.6798E
-
8
A(47)
= -
0.39518E
-
7: A(48) = 0.275996E
-
6: A(49) =
-
0.2475448E
-
5
A(50) = 0.3202967E
-
4: A(51) =
-
0.7552029E
-
3: A(52)= 0.0608819E0
Pi = 3.14159265359 ' число
"
Пи
"
'
S2Pi = 1/(Sqr(2*Pi)): Z = Abs(X)
If Z < 8 Then H = Z/8: K = 0
If Z >= 8 Then H = 8/Z: K = 26
L = K+13: Y = 4*H*H
-
2
RK(1) = A(1+K): RK(2) = Y*RK(1)+A(2+K)
RL(1) = A(1+L): RL(2) = Y*RL(1)+A(2+L)
For I = 3 To 13
RK(I) = Y*RK(I
-
1)
-
RK(I
-
2)+A(I+K)
RL(I) = Y*RL(I
-
1)
-
RL(I
-
2)+A(I+L)
Next I
F = S2Pi: D = F*RK(13): E = F*RL(13)*H: B = Sqr(Z)
If Z < 8 T
hen
C = D*B: S = E*B
Else
R = Sin(Z): T = Cos(Z)
C = 0.5+(D*R
-
E*T)/B
S = 0.5
-
(E*R+D*T)/B
End If
End Sub
'
17
'
-------------------------------------------------------------
' Вычисление значений интегралов Френеля S(x) и C(x)
'
Dim As Single X, X
1, C, S
Const Pi = 3.14159265359 ' число
"
Пи
"
Open "Con" For Output As #2
'Open "IntegrFren_SF45R.res" For Output As #2
Print #2, " X X1 C S"
For X = 1 To 15 Step 2 ' значение
аргумента
X1 = Sqr(2*X/Pi) ' пересчет
аргумента под С2 и S2
SF45R (X, C, S) ' вызов подпрограммы
Print #2, Using " ###.# ###.###### ###.######## ###.########";_
X; X1; C; S
Next X
Close #2
Sleep
' Результаты вычислений:
' X X1 C S
' 1.0 0.797885 0.72170556 0.24755824
' 3.0 1.381977 0.56102014 0.71168500
' 5.0 1.784124 0.32845652 0.46594140
' 7.0 2.111004 0.59011602 0.38119435
' 9.0 2.393654 0.56080395 0.61721361
' 11.0 2.
646284 0.38039190 0.50478631
' 13.0 2.876814 0.54251039 0.39826772
' 15.0 3.090194 0.56933606 0.57580328
Кривые интегралов Френеля
18
Параметрический график S(x) и C(x) даёт кривую на плоскости, называемую спиралью Корню или клото
идой. Спираль Корню была придумана Мари Альфредом Корню для облегчения расчёта дифракции в прикладных задачах.
Корню Мари Альфред
(Cornu Marie Alfred; 1841 -
1902) -
французский физик, иностранный чл.
-
корр. Петербургской АН. Для решения некоторых задач д
ифракции света Корню использовал (1874) кривую, состоящую из двух ветвей, симметричных относительно осей координат (спираль Корню), усовершенствовал метод определения скорости света французского физика И. Физо, уточнил плотность Земли (1873), обнаружил ано
мальный эффект Зеемана (1898).
Любители созерцать результаты в виде кривых на плоскости могут заменить в приведенной выше программе "IntegrFren_SF45R" (10.10.2012) блок вывода на следующие строки:
'
'
---------------------------------------------------
----------
' Вычисление значений интегралов Френеля S(x) и C(x)
'
Dim As Single X, X1, C, S
Const Pi = 3.14159265359 ' число "Пи"
Screen 12 ' вывод на граф. экран 640x480
' подготовка поля вывода
View (1, 1)
-
(478, 478), , 15 ' квад
ратное физическое окно
' самое большое для данного режима с белым кантиком...
Window (0, 0)
-
(1, 1) ' логическое окно вывода
For X = 0 To 1 Step 0.125
Line (0, X)
-
(1, X), 6, , &HAAAA ' горизонтали
сетки
Line (X, 1)
-
(X, 0), 6, , &HAAAA ' вертикали се
тки
Next X
Line (0, 0)
-
(0, 0), 7, , &HAAAA ' линия
оси
X
Line (0, 1)
-
(0, 0), 7, , &HAAAA ' линия
оси
Y
Locate 2, 2: Print "Y" ' метка
оси
Y
Locate 29,61: Print "X" ' метка
оси
X
' вывод спирали Корню (параметрический график S(x) и C(x))
Fo
r X = 0 To 1000 Step 0.005 ' значение аргумента
X1 = Sqr(2*X/Pi) ' аргумент как в примере
SF45R (X1, C, S) ' вызов подпрограммы
PSet(C, S), 13 ' вывод очередной точки графика
Next X ' следующая точка
Sleep
'
19
4.2.
Интерполяция кубическим сплайном табличной функции
Библиотека численного анализа
представляет собой большой пакет вычислительных программ, подготовленный в отделе численного анализа Научно
-
исследовательского Вычислительн
ого Центра Московского Государственного Университета им. М.В.Ломоносова (НИВЦ МГУ) под научным руководством В.В.Воеводина. В данной Библиотеке объединены вычислительные подпрограммы на стандарте языка Фортран, удовлетворяющие единым требованиям, обеспечива
ющим согласованность Библиотеки.
Библиотека доступна на сайте:
http://num
-
anal.srcc.msu.ru/lib_na/libnal.htm
и включает большой раздел:
http://num
-
anal.srcc.msu.ru/lib_na/cat/cat9.htm
Интерполяция, аппроксимация, сглаживание, численное дифференцирование
h
ttp://num
-
anal.srcc.msu.ru/lib_na/cat/cat91.htm Интерполяция, включая:
Вычисление разностей;
Полиномиальная интерполяция;
Интерполяция рациональными функциями;
Сплайн
-
интерполяция.
и далее по теме Фрагментов:
http://num
-
anal.srcc.msu.ru/lib_na/cat/cat914.
htm
Сплайн
-
интерполяция.
http://num
-
anal.srcc.msu.ru/lib_na/cat/i/iis3r.htm
Сплайн
-
интерполяция. Подпрограмма: IIS3R.
Назначение:
Интерполяция кубическим сплайном табличной функции от одной переменной на неравномерной сетке при известных краевых условиях на вторые производные.
На сайте и в сети нетрудно найти "теоретическое подтверждение" приведенному алгоритму (подпрограмме), поэтому ниже приведем лишь текст подпрограммы.
Литература:
Дж.Алберг, Э.Нильсон, Дж.Уолш. Теория сплайнов и ее приложения. M., 1972
, 320 с.
Монография посвящена изложению основ теории кусочно
-
полиномиальных приближений и некоторых ее применений. Это новое направление в теории приближений, которое в настоящее время усиленно развивается главным образом американскими математиками.
Исполь
зование.
Вызов процедуры:
IIS3R (X(), Y(), NX, BPAR(), C(), IERR)
Параметры:
NX
-
заданное число узлов сетки;
X
-
вектор длины NX узлов сетки по возрастанию;
Y
-
вектор длины NX значений интерполируемой функции;
C
-
массив размера 3*(NX
-
1), содержащий вычи
сленные коэффициенты
BPAR -
вектор длины 4, содержащий параметры краевых условий;
IERR –
код ошибки:
IERR=65 -
заданное число узлов сетки меньше 2;
IERR=66 -
элементы вектора X не упорядочены по возрастанию.
20
Вызываемые подпрограммы:
UTII10 -
подпрограмма вывода диагностических сообщений (
заменена заглушкой
).
Ниже приведен текст подпрограммы и вызывающей программы, соответствующий примеру использования на сайте, а также результат их выполнения.
'
' P R O G R A M "Spline_IIS3R"
' 10.10.2012
'
-------------------------------------------------------------
' Интерполяция кубическим сплайном табличной функции
'
-------------------------------------------------------------
'
#Lang "FB" ' режим FreeBA
SIC
-
совместимости
' Подпрограмма (заглушка) диагностических сообщений
Sub UTII10(IERR As Integer, N As Integer)
Print " *** UTII10: IERR, N"; IERR; N
End Sub
'
Sub IIS3R(X() As Single, Y() As Single, NX As Integer,_
BPAR() As Single, C() As Single, IER
R As Integer)
Dim AS Integer NXM1, J, I
Dim As Single DXJ, YPPB, PJ, SIXI, DXJP1, YPPA
Dim As Single DXP, DYJP1, DYJ, DX
Dim As Single ZERO, HALF, ONE, TWO, SIX
'EQUIVALENCE(DXJ,YPPB),(PJ,SIXI),(DXJP1,YPPA)
ZERO = 0.0: HALF = 0.5: ONE = 1.0: TWO = 2.0: SIX
= 6.0
IERR = 0: NXM1 = NX
-
1
If NX < 2 Then GoTo L6
If NX = 2 Then GoTo L2
DXJ = X(2)
-
X(1)
If DXJ <= ZERO Then GoTo L7
DYJ = Y(2)
-
Y(1)
For J = 2 To NXM1
DXJP1 = X(J+1)
-
X(J)
If DXJP1 <= ZERO Then GoTo L7
DYJP1=Y(J+1)
-
Y(J)
DXP=DXJ+DXJP1
C(1,J)
=DXJP1/DXP
C(2,J)=ONE
-
C(1,J)
C(3,J)=SIX*(DYJP1/DXJP1
-
DYJ/DXJ)/DXP
DXJ=DXJP1
DYJ=DYJP1
Next J
L2: C(1,1)=
-
BPAR(1)*HALF
C(2,1)=BPAR(2)*HALF
If NX = 2 Then GoTo L4
For J = 2 To NXM1
PJ=C(2,J)*C(1,J
-
1)+TWO
C(1,J)=
-
C(
1,J)/PJ
C(2,J)=(C(3,J)
-
C(2,J)*C(2,J
-
1))/PJ
Next J
21
L4: YPPB=(BPAR(4)
-
BPAR(3)*C(2,NXM1))/(BPAR(3)*C(1,NXM1)+TWO)
SIXI=ONE/SIX
For I = 1 To NXM1
J=NX
-
I
YPPA=C(1,J)*YPPB+C(2,J)
DX=X(J+1)
-
X(J)
C(3,J)=SIXI*(YPPB
-
YPPA)/DX
C(2,J)=HALF*YPPA
C(1,J)=(Y(J+1)
-
Y(J))/DX
-
(C(2,J)+C(3,J)*DX)*DX
YPPB=YPPA
Next I
GoTo L9
L6: IERR=65
GoTo L8
L7: IERR=66
L8: UTII10(IERR,3)
L9: End Sub
'
'
------------------------------------------------
-------------
' Основной модуль тестовой программы
' X -
вектор длины NX узлов сетки;
' Y -
вектор длины NX значений интерполируемой функции;
' C -
массив размера 3*(NX
-
1), содержащий вычисленные коэффициенты
Dim As Single X(4), Y(4), BPAR(4), C(3, 3)
Dim As Integer NX, I, J, IERR
Data 0.00, 0.33, 1.00, 2.00 ' аргументы узлов
NX = 4
For J = 1 To NX: Read X(J): Next J
For I = 1 To NX: Y(I) = ((X(I)
-
2)*X(I)+3)*X(I)+4: Next I
BPAR(1) = 1
BPAR(2) = 6/(X(2)
-
X(1))*((Y(2)
-
Y(1))/(X(2)
-
X(1))
-
3)
BPAR(3) = 1
BPA
R(4) = 6/(X(4)
-
X(3))*(7
-
(Y(4)
-
Y(3))/(X(4)
-
X(3)))
'
IIS3R(X(), Y(), NX, BPAR(), C(), IERR)
'
Print Using " IERR = ###.#####"; IERR
For J = 1 To 3
For I = 1 To NX
-
1: Print Using " ###.#####"; C(I,J);: Next I
Print
Next J
Sleep
'
' Результаты вычислений:
'
IERR = 0
' 3.00000 -
2.00000 1.00000
' 2.00670 -
1.01000 1.00000
' 2.00000 1.00000 1.00000
22
4.3.
Подпрограммы сплайн
-
интерполяция из любимой книжки
Речь идет
, конечно,
о великолепной книге:
Форсайт Дж., Малькольм М., Моулер К., Машинные методы математических вычислений. М., Мир, 1980, 280 с., илл.
которую рекомендую всем начинающим инженерам
–
она легко доступна в сети.
Открываем страницу 78: Глава 4. Интерполяция.
Предположим, что задано множество вещественных абсцисс:
x1, ... , xn и со
ответствующие ординаты y1, ... , yn.
Здесь x1 < x2 < ... < xn, а каждое yi -
некоторое вещественное число, отвечающее x
i, определенное математически или являющееся результатом какого
-
либо наблюдения. Задача одномерной интерполяции состоит в построении функ
ции f такой, что f(хi) = уi для всех i, и при этом f(x) должна принимать "разумные" значения для x, лежащих между заданными точками. Критерий разумности меняется от задачи к задаче и ему, возможно, никогда не будет дано точного определения...
Листаем до с
траницы 86: Пункт 4.4. Сплайн
-
интерполяция.
Кубические сплайн
-
функции это недавнее математическое изобретение, но они моделируют очень старое механическое устройство. Чертежники издавна пользовались механическими сплайнами, представляющими собой гибкие рей
ки из какого
-
нибудь упругого материала, обычно дерева или (в последнее время) пластика...
И, наконец, на странице 91: Пункт 4.5. Подпрограммы SPLINE и SEVAL
.
Метод вычисления параметров кубической сплайн
-
функции, обсуждавшийся сейчас, реализован в подпрог
рамме SPLINE...
С помощью подпрограммы
-
функции SEVAL можно вычислять значения сплайна после того, как его коэффициенты были определены подпрограммой SPLINE...
Ниже приведен текст подпрограмм, адаптированных под язык FreeBASIC
, максимально приближенный к к
нижному варианту.
' P R O G R A M "Spline_Fors"
' 10.10.2012
'
-------------------------------------------------------------
' Вычисление параметров SPLINE и значений SEVAL сплайна
'
---------------
----------------------------------------------
' Ординарная переработка программ из книги:
' Форсайт Дж., Малькольм М., Моулер К., Машинные методы
' математических вычислений. М., Мир, 1980, 280 с., илл.
#Lang "FB" ' режим FreeBASIC
-
совместимости
Su
b SPLINE(N As Integer, X() As Single, Y() As Single,_
B() As Single, C() As Single, D() As Single)
' Вычисляются коэффициенты B(I),C(I) и D(I) кубического
' интерполяционного сплайна для I = 1, 2,..., N
' S(X) = Y(I)+B(I)*(X
-
X(I))+C(I)*(X
-
X(I))^2
+D(I)*(X
-
X(I))^3
' Входные данные:
' N -
количество заданных точек (узлов) N >= 2
' X -
абсциссы узлов в возрастающем порядке
' Y -
соответствующие ординаты узлов
23
' Выходные данные:
' B(), C(), D() -
массивы коэффициентов сплайна.
' Если обоз
начить через P символ дифференцирования, то
' Y(I) = S(X(I))
' B(I) = SP(X(I))
' C(I) = SPP(X(I))/2
' D(I) = SPPP(X(I))/6 (правосторонняя производная)
' С помощью функции SEVAL можно вычислить значения сплайна.
Dim As Integer NM1, IB, I
Dim As Si
ngle T
NM1 = N
-
1
If (N < 2) Then Return
If (N < 3) Then GoTo L50
'
' Построить трехдиагональную систему, где
' B -
диагональ, D -
наддиагональ, C -
правые части.
D(1) = X(2)
-
X(1)
C(2) = (Y(2)
-
Y(1))/D(1)
For I = 2 To NM1
D(I) = X(I+1)
-
X(I)
B(I) = 2.*(D(I
-
1)+D(I))
C(I+1) = (Y(I+1)
-
Y(I))/D(I)
C(I) = C(I+1)
-
C(I)
Next I
' Граничные условия. Третьи производные в точках X(1) и X(N)
' вычисляются с помощью разделенных разностей
'
B(1) = -
D(1): B(N) = -
D(N
-
1)
C(1) = 0: C(N) = 0
If(N = 3) Then GoT
o L15
C(1) = C(3)/(X(4)
-
X(2))
-
C(2)/(X(3)
-
X(1))
C(N) = C(N
-
1)/(X(N)
-
X(N
-
2))
-
C(N
-
2)/(X(N
-
1)
-
X(N
-
3))
C(1) = C(1)*D(1)^2/(X(4)
-
X(1))
C(N) = -
C(N)*D(N
-
1)^2/(X(N)
-
X(N
-
3))
'
' Прямой ход вычислений
L15: For I = 2 To N
T = D(I
-
1)/B(I
-
1)
B(I) = B(I)
-
T*D
(I
-
1)
C(I) = C(I)
-
T*C(I
-
1)
Next I
'
' Обратная подстановка
C(N)=C(N)/B(N)
For IB = 1 To NM1
I = N
-
IB
C(I) = (C(I)
-
D(I)*C(I+1))/B(I)
Next IB
' C(I) теперь хранят величины SIGMA(I), определенные в п.4.4
'
' Вычисление коэ
ффициентов полиномов
B(N)=(Y(N)
-
Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.*C(N))
24
For I = 1 To NM1
B(I) = (Y(I+1)
-
Y(I))/D(I)
-
D(I)*(C(I+1)+2.*C(I))
D(I) = (C(I+1)
-
C(I))/D(I)
C(I) = 3.*C(I)
Next I
C(N) = 3.*C(N)
D(N) = D(N
-
1)
Return
'
L50: B(1) = (Y(2)
-
Y(1))/(
X(2)
-
X(1)): C(1) = 0: D(1) = 0
B(2) = B(1): C(2) = 0: D(2) = 0
Return
End Sub
'
Function SEVAL(N As Integer, U As Single, X() As Single, Y() As Single, B() As Single, C() As Single, D() As Single) As Single
' Подпрограмма
-
функция вычисляет значе
ние кубического сплайна
' SEVAL = Y(I)+B(I)*(U
-
X(I))+C(I)*(U
-
X(I))^2+D(I)*(U
-
X(I))^3
' где X(I) < U < X(I+1). Используется схема Горнера
' если U < X(1), то берется значение I = 1,
' если U >= X(N), то берется значение I = N.
' Входные данные:
' N -
количество заданных точек
' U -
абсцисса, для вычисле
ния значени
я
сплайна
' X, Y -
массивы заданных абсцисс и ординат
' B, C, D -
массивы вычисленных подпрограммой SPLINE
' коэффициентов сплайна
'
' Если по сра
внению с предыдущим вызовом U не находится
' в том же интервале, то для отыскания нужного интервала
' применяется двоичный поиск.
'
Dim As Integer I, J, K
Dim As Single DX
I = 1 ' было
DATA I/1/
If(I >= N) Then I=1
If(U < X(I)) Then GoTo L10
If(
U <= X(I+1)) Then GoTo L30
'
' Выполнение двоичного поиска
L10: I=1
J=N+1
L20: K=(I+J)/2
If(U < X(K)) Then J=K
If(U >= X(K))Then I=K
If(J > I+1) GoTo L20
'
' Вычислить сплайн
L30: DX = U
-
X(I)
SEVAL = Y(I)+DX*(B(I)+DX*(C(I)+DX
*D(I)))
End Function
25
'
'
-------------------------------------------------------------
' Основной модуль тестовой программы
'
Dim As Single X(10), Y(10), B(10), C(10), D(10)
Dim As Single S, U
Dim As Integer I, N
N = 10
For I = 1 To N
X(I) = I ' зна
чения
аргументов
Y(I) = X(I)^3 ' значения тестовой функции
Next I
SPLINE(N, X(), Y(), B(), C(), D())
Open "CON" For Output As #2
'Open "Spline_Fors.res" For Output As #2
For U = 1.0 To 3.0 Step 0.5
S = SEVAL(N, U, X(), Y(), B(), C(), D())
Print #2, Using
" ##.### ###.######"; U; S
Next U
Close #2
Sleep
' Результаты вычислений:
' 1.000 1.000000
' 1.500 3.375000
' 2.000 8.000000
' 2.500 15.625000
' 3.000 27.000000
Прочитал весь опус и ужаснулся –
нет, праздничные дни вовсе не предназна
чены для подготовки "научных статей"!
Нужно сделать пару итераций, чтобы было "вкусно", но, пожалуй, предоставлю такую возможность читателям...
На этом позвольте закруглиться.
До новых встреч!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
254
Размер файла
864 Кб
Теги
информатика freebasic программирование
1/--страниц
Пожаловаться на содержимое документа