close

Вход

Забыли?

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

?

FreeBASIC15

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. Работа состоит из нескольких глав. В главе 15 "Аппроксимации функций" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также препода
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 15 "Аппроксимации функций" –
продолжение пособи
я, которое будет интересно для учащихся школ, студентов
институтов, а также преподавателей.
Фрагмент 15. Аппроксимации функций
Хочу предложить три цитаты, которые, на мой взгляд,
могут послужить ориентирами в нашей повседневной жизни:
У России нет друз
ей, нашей огромности боятся... У России только два надежных союзника -
ее армия и ее флот.
Александр III
2
В России любят затевать реформы только потому, что так легче скрыть неумение править!
Петр Столыпин
Сегодня существует гораздо более реалистичная о
ценка России –
как региональной силы, но уже не как глобальной сверхдержавы, каковой она действительно больше не является.
Пол Гобл
Справки:
Александр III Александрович (1845 —
1894)
—
Император всероссийский, Царь Польский и Великий Князь Финляндский
, из императорского дома Романовых -
сын Императора Александра II и внук Николая I, отец русского Императора Николая II. В официальной дореволюционной историографии именовался Миротворцем.
Столыпин Петр Аркадьевич (1862 —
1911)
—
государственный
деятель Российской империи. В разные годы занимал посты уездного предводителя дворянства, губернатора Гродненской и Саратовской губернии, министра внутренних дел, премьер
-
министра. В российской истории начала XX века известен в первую очередь как реформат
ор и государственный деятель, сыгравший значительную роль в подавлении революции 1905
—
1907 годов.
3
Пол Гобл -
американский политолог
-
"кремлинолог"
, аналитик ЦРУ и госдепартамента, пишет статьи в газете "Georgian Daily". В России известен своими а
нтироссийскими проектами на Кавказе, в последнее время очень много пишет о Татарии. Вряд ли это случайность -
такие профи случайно ничего не делают. Впрочем, и без того понятно, что Татария и Урало
-
Поволжский регион в целом -
одно из главных направлений уд
ара. Достаточно напомнить о конфликтных ситуациях в Сургуте Ханты
-
Мансийский АО и в городе Первоуральске Свердловской области... За этими событиями, вполне возможно, стоит отнюдь не бытовой конфликт и даже не межнациональная неприязнь. Все громче звучат го
лоса тех, кто предупреждает: сибирские регионы вполне могут стать базой для радикальных исламистов –
ваххабитов. Мнение экспертов о возможности создания так называемого "Сибирского Халифата"...
В начале мая вашингтонский политолог Пол Гобл поделился с Рус
ской службой "Голоса Америки" анализом происходящих в России процессов и дал оценку американо
-
российским отношениям:
http://www.golos
-
ameriki.ru/content/paul
-
goble/1658470.html
Я думаю, мы в США, наконец, движемся в том направлении, где Россия больше не занимает центральное место во внешней политике. Очень долго мы жили отголосками "холодной войны", когда весь мир вращался вокруг происходящего в Вашингтоне и в Москве и развития о
тношений между США и Россией. Сегодня существует гораздо более реалистичная оценка России –
как региональной силы, но уже не как глобальной сверхдержавы, каковой она действительно больше не является. Мне также кажется, что сейчас пришло более четкое понима
ние того, что все наши надежды начала 1990
-
х годов на то, что Россия станет демократической страной, страной свободного рынка, где государство защищало бы права и свободы своих граждан и обеспечило бы долю всех жителей страны в доходах от громадных природн
ых богатств –
все эти надежды перечеркнуты
...
Сейчас в России происходит видимое ослабление позиций этнических русских
на фоне усиления позиций нерусских. Это объективный процесс, обусловленный многими факторами, в числе которых –
низкая рождаемость у рус
ских и высокая рождаемость у нерусских, миграция из стран Центральной Азии и Южного Кавказа –
в комплексе ведущий к резкому падению процентной доли этнических русских в общей численности населения России...
Стало совершенно очевидным, что сочинская Олимпи
ада, которую Владимир Путин планировал как символ и высочайшее достижение своего президентства, скорее всего обернется его крупнейшим провалом. Я говорю это не потому, что там обязательно случится какое
-
то насилие, и не потому, что ожидаю широкого мирового
бойкота –
а исходя из опыта Игр в других странах. Когда вы планируете международное событие, которое привлекает внимание всего мира, люди видят не только то, что вы им показываете, но и то, что вы очень хотели бы скрыть
...
4
1.
"Политика и Искусство".
Вы
ше, как вы понимаете, было представлено мнение специалиста из Америки. В стране же запущена кремлевская идеологическая машина
под названием экспертный "Изборский клуб". Это периодическое собрание известных экспертов, специализирующихся на изучении внешней и внутренней политики России. Клуб был создан в сентябре 2012 во время празднования 1150
-
летия города Изборск. Председатель клуба —
Александр Андреевич Проханов (род. 1938) —
советский и российский политический деятель, писатель, публицист. Член секретариа
та Союза писателей России, главный редактор газеты "Завтра". Вот что он сказал в своем выступлении:
Изборский клуб собрал в себя два десятка народно
-
патриотически мыслящих интеллектуалов, антилибералов, с целью разработать идеологию и стратегию рывка, о к
отором обмолвился Путин и который необходим для того, чтобы ликвидировать чудовищное отставание России во всех областях. В том числе и в создании обновленного оборонно
-
промышленного комплекса и нового русского оружия, поскольку мы сейчас находимся накануне
огромной войны.
Но КПД запущенной машины оказался не слишком высоким
и Кремль наспех сооружает самодельные "идеологические подпорки"! Например, многие политологи рассматривают предстоящие Олимпийские игры в Сочи в 2014 году как часть идеи "возрождения в
еликой России" (см. выступление Пола Гобла). Конечно, интересна и идея создания "Средиземноморской эскадры", которая будет действовать не только в Средиземном море, но и в Мировом океане. Штаб средиземноморской эскадры будет сформирован уже этим летом, он будет находиться на одном из флагманских кораблей в Средиземном море. Состав соединения пока точно не определен, в него могут войти подводные лодки (в том числе АПЛ), госпитальные суда и вертолетоносцы типа "Мистраль"... Теперь держись "Шестой флот США" (U
nited States Sixth Fleet —
оперативный флот американских военно
-
морских сил, дислоцирующийся в Средиземном море)! Но, как говорят: "дорога ложка к обеду", а "обед" в том регионе уже подходит к концу и финал его –
справедливое недовольство "братьев по оружи
ю" (начиная с Югославии и заканчивая Сирией)... Кроме того, эта "героическая" выходка может дорого стоить впоследствии при противостоянии со сторонниками "Сибирского Халифата"...
А самое смешное, что по времени это совпало с двумя знаменательными события
ми Российского флота, которые по разным причинам не слишком афишировались СМИ. Потому и были приведены цитаты Александра III и Петра Столыпина (наверное, подумали, что упомянул "для красного словца")
,
а
две картины традиционной "художественной галереи"
пре
дставлены ниже
...
5
Взятие шведских кораблей в устье Невы.
Картину написал Блинов Леонид Демьянович
(1867 —
1903) художник
-
маринист, родился 4 января 1867 года в Подмонастырской слободе Ярославской губернии. В 1886 Леонид Блинов был принят вольнослуш
ателем на отделение живописи в Академию художеств и в том же году от Академии художеств и Главного морского штаба он был назначен в кругосветное плавание на пароходе "Москва". На картине показан бой
—
морское сражение, состоявшееся 18 мая 1703 года в устье
реки Невы между тремя десятками лодок с солдатами Преображенского и Семеновского лейб
-
гвардейских полков под командованием Петра Алексеевича (Петра I) и его сподвижника Александра Даниловича Меншикова и двумя небольшими кораблями "Гедан" и "Астрильд" швед
ского флота, пришедшими в составе эскадры на помощь крепости Ниеншанц. В ходе непродолжительного боя оба шведских корабля были взяты на абордаж. Условно бой можно назвать первым морским сражением русского флота. Каждый моряк тогда получил именную медаль с надписью "Небываемое бывает". Дата 18 мая считается датой рождения Балтийского флота -
310 лет назад была одержана первая победа России в морском бою.
Вход кораблей русской эскадры в бухту Севастополя.
6
Картину написал Псарев Виктор Пантелеевич
, родился в 1950 году в городе Лениногорске Восточно
-
Казахстанской области. Окончил Московский государственны
й
художественный институт имени В.И.Сурикова в 1980. В Студии военных художников имени М.Б.Грекова с 1980 года. В 1995 году присвоено звание "Народный художн
ик РФ". В 2000 году награжден орденом Русской Православной церкви Преподобного Сергия Радонежского III степени. Народный художник РФ, Член
-
корреспондент Российской академии художеств, Лауреат премии Ленинского комсомола. Содержание картины связано с присое
динением Крыма к России.
В январе 1783 правительство назначило для командования флотом, заводимым на Чёрном и Азовском морях, вице
-
адмирала Федота Алексеевича Клокачёва (1739 —
1783). Ему было предписано обследовать Ахтиарскую бухту, и, при установлении пр
игодности её для образования флота, перевести туда часть Азовской и Днепровской флотилий. В апреле 1783 года рескриптом Екатерины II Крым был включён в состав России. 13 мая в Ахтиарскую, а ныне Севастопольскую бухту вошла эскадра Клокачёва с тех пор 13 ма
я считается датой рождения Черноморского флота -
230 лет назад в Севастопольскую бухту вошла Российская эскадра.
1.
Почтовый ящик.
Несмотря на череду изнуряющих праздников, "жизнь научного сообщества" продолжается! Об этом свидетельствует письмо в поч
товом ящике от "постоянного читателя и скачивателя" (как он себя назвал)
SS
уже знакомому по сообщениям в брошюрах "Фрагмент 12" и "Фрагмент 13". В этом письме он просит продолжить публикацию на языке FreeBASIC небольших программ из книги:
Мудров А.Е. Чис
ленные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль.
Томск: МП "РАСКО", 1991. 272 с: ил.
и, прежде всего, программ предназначенных для аппроксимации функций.
Желание читателя –
закон!
В настоящей брошюре постараемся кратко рассмотреть эту важную те
му, применительно к аппроксимации функции одной переменной.
2.
Метод аппроксимации. Введение.
Аппроксимация
(от лат. approximo -
приближаюсь) или приближение —
научный метод, состоящий в замене одних объектов другими, в том или ином смысле близкими к исходным, но более простыми. Аппроксимация позволяет исследовать числовые характеристики и качественные свойства объекта, сводя задачу к изучению объектов, характеристики которых легко вычисляются или свойства которых уже известны. Близкое по смыслу поняти
е -
абстрагирование
(от лат. abstractio —
отвлечение) или мысленное выделение, вычленение некоторых элементов конкретного множества и отвлечение их от прочих элементов данного множества. Это один из основных процессов умственной деятельности человека. 7
2.1.
Аппроксимация функций.
Под аппроксимацией понимается задача нахождения для некоторой функции более простой приближающей функции близкой, в определенном смысле, к функции. Теория аппроксимации в современном виде разработана гениальным русским матем
атиком и механиком Пафнутием Львовичем Чебышёвым
Его работа "О функциях, наименее уклоняющихся от нуля" (приложение к "Запискам Академии Наук". СПб., 1873, том XXII), как по постановке вопроса, так и по методу решения, уже является поводом для обязательног
о упоминания Чебышева. В течение сорока лет Чебышёв принимал активное участие в работе военного артиллерийского ведомства. Он работал над усовершенствованием приемов артиллерийской стрельбы, как по дальнобойности, так и по точности. В результате, как "побо
чный продукт", была построена новая модель суммирующей машины, переданная (1878) в Парижский музей искусств и ремесел, а затем создана множительно
-
делительная приставка к суммирующей машине. Эта приставка также была передана в музей в Париже (1881). С меха
низмами, большей частью уникальными, созданными Чебышевым можно ознакомиться в виртуальном музее
фонда "Математические этюды" на сайте:
http://tcheb.ru/
Теоретическую основу создает теорема Вейерштрасса
, утверждающая, чт
о для каждой "хорошей" функции на отрезке найдется многочлен, отличающийся от нее в любой точке отрезка на величину, не превосходящую допустимого отклонения. Это действительно так. Беда лишь в том, что теорема Вейерштрасса ничего не говорит о способе постр
оения приближающего, или аппроксимирующего, многочлена. Теорема Вейерштрасса принадлежит к числу так называемых чистых теорем существования; утверждая, что приближающий многочлен существует, она умалчивает о том, как его построить. Не следует думать, что т
еоремы существования бесполезны: они говорят нам, что множество многочленов, отличающихся от интересующей нас функции на заданном отрезке, непусто (заранее это не очевидно). Следовательно, мы можем надеяться, что нам удастся найти приближающий многочлен. Е
сли бы была доказана теорема о том, что приближающий многочлен не существует, то от попыток построить его пришлось бы заранее отказаться.
Если говорить о принципиальной стороне дела, то знание ряда Тейлора
действительно позволяет вычислить значение функци
и в любой точке из области ее определения с любой точностью. Но практическая сторона таких вычислений оставляет желать лучшего: для достижения требуемой точности иногда приходится брать много членов ряда Тейлора (в таких случаях говорят, что ряд Тейлора сх
одится медленно). Возникает вопрос, который П.Л. Чебышев назвал общим для всей практической деятельности человека: как располагать средствами своими для достижения по возможности большей выгоды? Применительно к задаче о вычислении значений функции это озна
чает: как найти многочлен, приближающий заданную функцию с требуемой точностью, но более короткий, чем отрезок ряда Тейлора, обеспечивающий ту же точность приближения
?
При этом в качестве критерия согласия могут быть использованы три условия:
1.
точное совпа
дение значений искомой функции с "экспериментом" в узлах таблицы (критерий интерполяции)
;
2.
сумма квадратов отклонений значений искомой и табличной функций минимальна (критерий среднеквадратической аппроксимации)
;
8
3.
максимальное по абсолютной величине из откло
нений значений искомой и табличной функций минимально (критерий равномерной аппроксимации)
.
Когда и как инженер должен использовать тот или иной критерий –
предмет книги среднего размера, а потому не будем сосредотачиваться на этом...
До начала XIX века учёные не имели опредёленных правил для решения системы уравнений, в которой число неизвестных меньше числа уравнений; до этого времени употреблялись частные приёмы, зависевшие от вида уравнений и от остроумия вычислителей, и потому разные вычислители, исх
одя из тех же данных наблюдений, приходили к различным выводам. Лежандру (1805) и Гауссу (1794)
принадлежит первое применение к решению указанной системы уравнений теории вероятности, исходя из начал, аналогичных с началом арифметической середины, уже изда
вна и, так сказать, бессознательно применяемых к выводам результатов в простейшем случае многократных измерений. Как и в случае арифметической середины, вновь изобретённый способ не даёт, конечно, истинных значений искомых, но даёт зато вероятнейшие значен
ия.
Этот способ распространён и усовершенствован дальнейшими изысканиями Лапласа, Энке, Бесселя, Ганзена и др. и получил название метода наименьших квадратов (
МНК, OLS, Ordinary Least Squares
), потому что после подстановки в начальные уравнения неизвестны
х величин, выведенных этим способом, в правых частях уравнений получаются если и не нули, то небольшие величины, сумма квадратов которых оказывается меньшей, чем сумма квадратов подобных же остатков, после подстановки каких бы то ни было других значений не
известных. Помимо этого, решение уравнений по способу наименьших квадратов даёт возможность выводить вероятные ошибки неизвестных, то есть даёт величины, по которым судят о степени точности выводов. Сущность МНК (обычного, классического) заключается в том,
чтобы найти такие параметры уравнений, при которых сумма квадратов остатков RSS (англ. Residual Sum of Squares) будет минимальной.
2.2.
Метод наименьших квадратов.
Обратимся к книге Анатолия Евстигнеевича Мудрова. На странице 105 читаем:
Глава 4. Мето
д наименьших квадратов.
Здесь представляют интерес три небольшие программы, поясняющие три варианта вычислений методом наименьших квадратов:
4.1. МНК со степенным базисом
4.3. МНК с базисом в виде ортогональных полиномов
4.4. МНК с базисом в виде полиномов
дискретной переменной
2.2.1.
МНК со степенным базисом.
В этом варианте метода наименьших квадратов базисные функции выбираются в виде последовательности степеней аргумента
, которые линейно независимы. В этом случае так же, как и при интерполяции, эксп
ериментальная зависимость аппроксимируется полиномом. Однако степень полинома m выбирается обычно m <= n (при лагранжевой интерполяции m = n). Аппроксимирующая кривая в МНК не проходит через значения исходной функции в узлах, но проводится из условия 9
наиме
ньшего суммарного квадратичного отклонения. Экспериментальные данные "сглаживаются" с помощью базисных функций. В основном модуле программы для заданных значений аргумента X(I), где I = 0, 1,...N, вычисляются значения аппроксимируемой функции F(I) = Sin(Pi
*X(I)/6).
' P R O G R A M "OrPol_MNK_sb"
' 10.10.1996
'
-------------------------------------------------------------
' Вариант метода наименьших квадратов (МНК)
'
-----------------------------
--------------------------------
'
' Алгоритм МНК со степенным базисом. См. книгу:
' Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик,
' Фортран и Паскаль. Стр. 110: Программа 4.1
' N –
максимальный номер узла таблицы;
' M –
максимальный номер базисной функции.
' Используются подпрограммы:
' GRAM -
формирование матрицы Грама A(M, M+1);
' GAUSS –
решение СЛАУ методом Гаусса.
#Lang "FB" ' Режим совместимости с FreeBASIC
' Исходные данные размещаются в общей области:
Commo
n Shared X() As Single
Common Shared F() As Single
'
Sub AFUN(M As Integer, C() As Single, X1 As Single, ByRef P As single)
' Вычисление аппроксимирующей функции по схеме Горнера
Dim As Integer I
P = C(M)
For I = M
-
1 To 0 Step -
1: P = C(I) + P*X1: Next I
E
nd Sub
'
'
Sub GRAM(N As Integer, M As Integer, A() As single)
' Формирование матрицы Грама A(M, M+1)
Dim As Integer I, J
Dim As Single P, S, R, Q
For J = 0 To M
S = 0: R = S: Q = S
For I = 0 To N
P = X(I)^J
S = S + P: R = R+P*F(I): Q = Q+P*X(I)^M
Next I
A
(0,J) = S: A(J,M+1) = R: A(J,M) = Q
Next J
For I = 1 To M
For J = 0 To M
-
1: A(I,J) = A(I
-
1,J+1): Next J
Next I
End Sub
'
'
10
Sub GAUSS(M As Integer, A() As Single, C() As Single)
' Метод Гаусса для СЛАУ
Dim As Integer M1, K, K1, I, J
Dim As Single R, S
M1 = M+1
For K = 0 To M: K1 = K+1: S = A(K,K)
For J = K1 To M1: A(K,J) = A(K,J)/S: Next J
For I = K1 To M: R = A(I,K)
For J = K1 To M1: A(I,J) = A(I,J)
-
A(K,J)*R: Next J
Next I
Next K
For I = M To 0 Step -
1: S = A(I,M1)
For J = I+1 To M: S = S
-
A(I,J)*C(J): Next J
C(I) = S
Next I
End Sub
'
' Основной модуль программы
Dim As Integer N, M, I
Dim As Single X0, X9, X1, H, P, Pi
' Входные параметры:
Data 3, 2, 0, 3, 0.5
Read N, M, X0, X9, H
Dim As Single X(N), F(N), A(10,11), C(10)
Pi = 3.1415
Open Cons For Output As #2
'Open "OrPol.res" For Ou
tput As #2
Print #2,: Print #2, " Initial parameters of the problem:"
Print #2, Using " ## ## ##.### ##.### ##.###"; N; M; X0; X9; H
' формирование таблицы исходных данных X(I), F(I):
Data 0, 1, 2, 3
For I = 0 To N: Read X(I): F(I) = Sin(Pi*X(I)/6): N
ext I
Print #2,: Print #2, " Source data table:"
' таблица
исходных
данных
For I = 0 To N: Print #2, Using " ###.#### ###.######"; X(I); F(I): Next I
GRAM(N, M, A())
GAUSS(M, A(), C())
Print #2,: Print #2, " Coefficients of the approximating function:"
' коэффициенты
аппроксимирующей
функции
For I = 0 To M: Print #2, Using " ###.#######"; C(I): Next I
Print #2,: Print #2, " Computed values of the functions:"
' вычисленные значения функции
For X1 = X0 To X9 Step H
AFUN(M, C(), X1, P)
Print #2, Using " ###.#
### ###.######"; X1; P: Next X1
Close #2
Sleep
'
'
11
' Результат выполнения программы:
' Initial parameters of the problem:
' 3 2 0.000 3.000 0.500
' Source data table:
' 0.0000 0.000000
' 1.0000 0.499987
' 2.0000 0.866010
' 3.0000
1.000000
' Coefficients of the approximating function:
' -
0.0049035
' 0.6110998
' -
0.0914992
' Computed values of the functions:
' 0.0000 -
0.004903
' 0.5000 0.277772
' 1.0000 0.514697
' 1.5000 0.705873
' 2.0000 0.851299
' 2.5000 0.950976
' 3.0000 1.004903
'
Расчет коэффициентов аппроксимирующей функции C(I) следует признать вполне удовлетворительным.
2.2.2.
МНК с базисом в виде ортогональных полиномов
Выбор базисных функций в виде степеней аргумента не является оптимальным с точки зрения решения системы н
ормальных уравнений с наименьшими погрешностями. Приемлемые результаты можно получить, если набор экспериментальных данных аппроксимируется полиномом невысокой степени (m <= 4 –
5). При этом лучшие результаты, как правило, дает использование классических о
ртогональных полиномов
. В программе приведенной ниже переменная L служит для выбора необходимого типа базиса. В основном программа повторяет приведенную выше программу.
' P R O G R A M "OrPol_MNK_op"
' 12.10.1996
'
-------------------------------------------------------------
' Вариант метода наименьших квадратов (МНК)
'
-------------------------------------------------------------
'
' Алгоритм МНК с базисом в виде классических ортогональных
'
полиномов. См. книгу:
' Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик,
' Фортран и Паскаль. Стр. 118: Программа 4.3
' L = 0 -
степени
' L = 1,2 -
полиномы Чебышева
' L = 3 -
полиномы Лежандра
'
'
12
#Lang "FB" ' Режим совместимости
с FreeBASIC
'
Common Shared X() As Single
Common Shared F() As Single
'
Sub BFUN(N As Integer, M As Integer, L As Integer, X1 As Single, T() As Single)
' Базисные функции T(M)
Dim As Integer K
Dim As Single Z, R
Z = 2*(X1
-
X(0))/(X(N)
-
X(0))
-
1: T(0) = 1: T(
1) = Z
For K = 1 To M
-
1: R = Z*T(K)
If L = 1 Then R = R
-
T(K
-
1)/4: GoTo L10
If L = 2 Then R = 2*R
-
T(K
-
1): GoTo L10
If L = 3 Then R = ((K+K+1)*R
-
K*T(K
-
1))/(K+1)
L10: T(K+1) = R
Next K
End Sub
'
'
Sub AFUN(N As Integer, M As Integer, L As Integer, C() As Sing
le, X1 As Single, ByRef S As single)
' Аппроксимирующая
функция
Dim As Integer I
Dim As Single T(10)
S = C(0)
BFUN(N, M, L, X1, T())
For I = 1 To M: S = S+C(I)*T(I): Next I
End Sub
'
'
Sub GRAM(N As Integer, M As Integer, L As Integer, A() As single)
' Формирование матрицы Грама A(M, M+1)
Dim A
s Integer I, J, K
Dim As Single T(10), P(10,50), X1, S, R, Q
For I = 0 To N: X1 = X(I)
BFUN(N, M, L, X1, T())
For J = 0 To M: P(J,I) = T(J): Next J
Next I ' P(M,N) -
матрица
плана
For K = 0 To M
For J = K To M: S = 0: R = 0
For I = 0 To N: Q = P(K,I): S =
S+Q*P(J,I)
If J = M Then R = R+Q*F(I) Next I
A(K,J) = S: A(J,K) = S
Next J
A(K,M+1) = R
Next K
End Sub
'
'
13
Sub GAUSS(M As Integer, A() As Single, C() As Single)
' Метод Гаусса для СЛАУ
Dim As Integer M1, K, K1, I, J
Dim As Single R, S
M1 = M+1
For K = 0 To M: K1 = K+1: S = A(K,K)
For J = K1 To M1: A(K,J) = A(K,J)/S: Next J
For I = K1 To M: R = A(I,K)
For J = K1 To M1: A(I,J) = A(I,J)
-
A(K,J)*R: Next J
Next I
Next K
For I = M To 0 Step -
1: S = A(I,M1)
For J = I+1 To M: S = S
-
A(I,J)*C(J): Next J
C(I) = S
Nex
t I
End Sub
'
' Основной модуль программы
Dim As Integer N, M, L, I
Dim As Single X0, X9, X1, H, S, Pi
' Входные параметры:
Data 3, 2, 3, 0, 3, 0.5
Read N, M, L, X0, X9, H
Dim As Single X(N), F(N), P(10,50), A(10,11), T(10), C(10)
Pi = 3.1415
Open Cons For Output As #
2
'Open "OrPol.res" FOR OUTPUT AS #2
Print #2,: Print #2, " Initial parameters of the problem:"
Print #2, Using " ## ## ## ##.### ##.### ##.###"; N; M; L; X0; X9; H
' формирование таблицы исходных данных X(I), F(I):
Data 0, 1, 2, 3
For I = 0 To N: Re
ad X(I): F(I) = Sin(Pi*X(I)/6): Next I
Print #2,: Print #2, " Source data table:"
' таблица
исходных
данных
For I=0 To N: Print #2, Using " ###.#### ###.######"; X(I); F(I): Next I
GRAM(N, M, L, A())
GAUSS(M, A(), C())
Print #2,: Print #2, " Coefficients of the approximating function:"
' коэффициенты
аппроксимирующей
функции
For I = 0 To M: Print #2, Using " ###.#######"; C(I): Next I
Print #2,: Print #2, " Computed values of the functions:"
' вычисленные значения функции
For X1 = X0 To X9 Step H
AFUN(N, M
, L, C(), X1, S)
Print #2, Using " ###.#### ###.######"; X1; S: Next X1
Close #2
Sleep
'
'
14
' Результат выполнения программы:
' Initial parameters of the problem:
' 3 2 0 0.000 3.000 0.500
' Source data table:
' 0.0000 0.000000
' 1.0000 0.499987
' 2.0000 0.866010
' 3.0000 1.000000
' Coefficients of the approximating function:
' 0.7058732
' 0.5049034
' -
0.2058732
' Computed values of the functions:
' 0.0000 -
0.004903
' 0.5000 0.277772
' 1.0000 0.514697
' 1.5
000 0.705873
' 2.0000 0.851300
' 2.5000 0.950976
' 3.0000 1.004903
' Initial parameters of the problem:
' 3 2 3 0.000 3.000 0.500
' Source data table:
' 0.0000 0.000000
' 1.0000 0.499987
' 2.0000 0.866010
' 3.00
00 1.000000
' Coefficients of the approximating function:
' 0.6372488
' 0.5049034
' -
0.1372487
' Computed values of the functions:
' 0.0000 -
0.004903
' 0.5000 0.277772
' 1.0000 0.514697
' 1.5000 0.705873
' 2.0000 0.851299
' 2.5000 0.950976
' 3.0000 1.004903
'
Тоже довольно неплохой результат.
2.2.3.
МНК с базисом в виде полиномов дискретной переменной
В зависимости от распределения погрешности обрабатываемых данных можно построить полиномы дискретной переменной, ортогональные с соответствующими дискретными весовыми функциями. Из классических ортогональных полиномов дискретной переменной наибольшую известность получили полиномы Чебышева, Кравчука, Шарлье, Мейкснера, Хана.
Для ознакомления с ними можно воспользоваться книгой:
15
Никифоров А.Ф., Суслов С.К., Увар
ов В.Б. Классические ортогональные полиномы дискретной переменной. М.: Наука, 1985. 215 с.
В книге на основе разработанного авторами простого подхода построена теория классических ортогональных полиномов дискретной переменной на равномерных и неравномерных
сетках. Здесь следует различать системы многочленов двух видов:
системы ортогональные на множестве равноотстоящих точек.
системы ортогональных на произвольном множестве точек.
Последовательность многочленов, ортогональных на конечном множестве точек дейс
твительной прямой R, впервые была введена и исследована в целом ряде работ П.Л.Чебышева в связи с задачами математической статистики. В дальнейшем было показано, что основные аналитические соотношения, определяющие свойства различных ортогональных базисов дискретной переменной из числа классических, являются дискретными аналогами соответствующих формул, задающих свойства классических ортогональных базисов непрерывной переменной. Это указывает на то, что многие характеристики ортогональных систем непрерывной
и дискретной переменной практически совпадают.
Обращение к ортогональным базисам дискретной переменной обусловлено, прежде всего, необходимостью избавления цифровых ЭВМ от численного интегрирования сложных функциональных зависимостей и заменой его суммиро
ванием при вычислении коэффициентов разложения.
Здесь нельзя не упомянуть книги:
Бейтмен Г., Эрдейи Э. Высшие трансцендентные функции. Том 2. М: Наука, 1966.
В отличие от других справочных пособий, книга содержит не только все формулы по теории специальн
ых функций, полученные к середине 40
-
годов, но и сжато изложенную теорию этих функций. По полноте охвата материала издание уникально.
Лоусон Ч., Хенсон Р. Численное решение задач метода наименьших квадратов. М.: 1986. 232 с.
Книга посвящена изложению числ
енных решений линейных задач метода наименьших квадратов. Достоинством книги являются: отбор наиболее устойчивых методов, полный анализ устойчивости, рассмотрение среднеквадратичных задач с линейными ограничениями, обзор методов перестройки ортогональных р
азложений при добавлении или удалении одного или нескольких наблюдений.
Ниже представлена программа
,
реализующая вариант метода наименьших квадратов для полиномов Чебышева ортогональных на системе точек.
' P R O G R A M "OrPol_MNK_dp"
' 12.10.1996
'
-------------------------------------------------------------
' Вариант метода наименьших квадратов (МНК)
'
-------------------------------------------------------------
'
' Полиномы Чебышева, орт
огональные на системе точек
' Алгоритм МНК с ортогональным базисом, см. книгу:
' Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик,
' Фортран и Паскаль. Стр. 125: Программа 4.4
16
' Пример 2: Демидович, Марон, Шувалова, с.37
' Используются подпрограммы:
' PCHE(N,M,X,Y,A,B,C) -
вычисление коэффициентов
' FUNC(M,A,B,C,X1,Y1) -
вычисление функции
'
#Lang "FB" ' Режим совместимости с FreeBASIC
Sub FUNC(M As Integer, A() As Single, B() As Single, C() As Single, X1 As Single, B
yRef Y1 As Single)
' Подпрограмма вычисления аппроксимирующей функции
Dim As Integer K
Dim As Single P, T, U, S
P = 1: T = 0: U = P: Y1 = C(0)
For K = 1 TO M
S = (X1 -
A(K)) * U -
B(K)*T
Y1 = Y1 + C(K)*S
T = U: U = S
Next K
End Sub
'
Sub PCHE(N As Integ
er, M As Integer, X() As Single, Y() As Single, A() As Single, B() As Single, C() As Single)
' Подпрограмма вычисления коэффициентов полиномов
'
Dim As Integer I, K, K1
Dim As Single P, Q, R, D, S, T(N), U(N)
P = 1
For I = 1 To N: T(I) = 0: U(I) = P: Nex
t I
A(0) = 0: B(0) = 0
For K = 0 To M
Q = 0: R = Q: S = Q
For I = 1 To N
D = U(I): R = R + D*Y(I): D = D*D: Q = Q + D: S = S + D*X(I)
Next I
C(K) = R/Q
If K = M Then Exit For K1 = K + 1
D = S/Q: A(K1) = D
R = Q/P: B(K1) = R
P = Q
For I = 1 To N
S = (X(I) -
D)*U(I) -
R*T(I)
T(I) = U(I): U(I) = S
Next I
Next K
End Sub
'
' Основной модуль программы
Dim As Integer N, M, I, K
Dim As Single X(99), Y(99), A(9), B(9), C(9), X1, Y1, D, S
'
Data 10
'
17
Data 0.0, 1.300
Data 0.3, 1.245
Data 0.6, 1.095
Data 0.9, 0.855
Data 1.2, 0.514
Data 1.5, 0.037
Data 1.8, -
0.600
Data
2.1, -
1.295
Data 2.4, -
1.767
Data 2.7, -
1.914
Read N
For I = 1 TO N
Read X(I), Y(I)
Next I
'
Open Cons For Output As #2
'Open "OrPol.res" For Output As #2
'
--
Вычисление значений коэффициентов полиномов
M = N -
1
PCHE (N, M, X(), Y(), A(), B(), C())
Pr
int #2,
Print #2, " M = N -
1 = "; M
Print #2,
Print #2, " K A B C"
Print #2,
For K = 0 To M
Print #2, Using " ## "; K;
Print #2, Using " ###.######"; A(K); B(K); C(K)
Next K
'
--
Вычисление значений аппроксимирующей функци
и
For M = 2 To N -
1
Print #2,
Print #2, " M = "; M
Print #2,
Print #2, " I X Y0 Y1 Y0
-
Y1"
Print #2,
S = 0
For I = 1 To N
X1 = X(I)
FUNC(M, A(), B(), C(), X1, Y1)
D = Y(I)
-
Y1
S = S + D*D
Print #2, Using " ## "; I;
P
rint #2, Using " ###.######"; X1; Y(I); Y1; D
Next I
Print #2,
Print #2, Using " S2 = ###.########"; S
Next M
'
Close #2
Sleep
18
' Результат выполнения программы:
' M = N -
1 = 9
' K A B C
' 0 0.000000 0.000000
-
0.053000
' 1 1.350000 10.000000 -
1.349535
' 2 1.350000 0.742500 -
0.315530
' 3 1.350000 0.576000 0.299714
' 4 1.350000 0.526500 0.261068
' 5 1.350000 0.480000 0.085840
' 6 1.35000
0 0.426136 -
0.093164
' 7 1.350000 0.362518 -
0.135554
' 8 1.350000 0.288346 -
0.072925
' 9 1.350000 0.203294 -
0.003904
' M = 2
' I X Y0 Y1 Y0
-
Y1
' 1 0.000000 1.3
00000 1.428100 -
0.128100
' 2 0.300000 1.245000 1.250421 -
0.005421
' 3 0.600000 1.095000 1.015947 0.079053
' 4 0.900000 0.855000 0.724677 0.130323
' 5 1.200000 0.514000 0.376612 0
.137388
' 6 1.500000 0.037000 -
0.028249 0.065249
' 7 1.800000 -
0.600000 -
0.489905 -
0.110095
' 8 2.100000 -
1.295000 -
1.008356 -
0.286644
' 9 2.400000 -
1.767000 -
1.583603 -
0.183397
' 10 2.7000
00 -
1.914000 -
2.215645 0.301645
' S2 = 0.28171551
' ...
' M = 9
' I X Y0 Y1 Y0
-
Y1
' 1 0.000000 1.300000 1.300000 0.000000
' 2 0.300000 1.245000 1.245000 0.000000
' 3 0.600000 1.095000 1.095001 -
0.000001
' 4 0.900000 0.855000 0.854998 0.000002
' 5 1.200000 0.514000 0.514000 0.000000
' 6 1.500000 0.037000 0.037001 -
0.000001
' 7 1.800000 -
0.6
00000 -
0.600002 0.000002
' 8 2.100000 -
1.295000 -
1.294999 -
0.000001
' 9 2.400000 -
1.767000 -
1.767001 0.000001
' 10 2.700000 -
1.914000 -
1.914000 -
0.000000
' S2 = 0.00000000
2.3.
Метод сингулярного разложения.
Редкая брошюра может обойтись без упоминания мой любимой книжки:
Форсайт Дж., Малькольм М., Моулер К., Машинные методы математических вычислений. М., Мир, 1980, 280 с.
Здесь она пригодится для изложения наиболее надежного метода для вычисления коэффициентов в общей задаче наименьших квадратов основан на матричной факторизации, называемой сингулярным разложением. Есть другие методы, требующие меньше машинного времени и па
мяти, однако они менее эффективны в 19
том, что касается учета ошибок заданной информации, ошибок округления и линейной зависимости. Но предварительно нужно привести еще одну маленькую программу из популярной книги А.Е.Мудрова. Речь идет о программе 2.3. для вычисления определителя системы линейных уравнений. Для вычисления определителя квадратной матрицы ее приводят к треугольному виду и находят произведение диагональных элементов. При этом следует учесть, что выбор главного элемента осуществляется перестанов
кой строк матрицы, каждая из перестановок изменяет знак определителя на противоположный. Ниже приведен ее текст на языке FreeBASIC.
' P R O G R A M "Appr_MDET"
' 12.08.2002
'
--------------------------
-----------------------------------
' Демонстрационная программа вычисления определителя
'
-------------------------------------------------------------
'
#Lang "FB" ' режим совместимости с FreeBASIC
Sub MDET(N As Integer, A() As Single, ByRef D
ET As Single)
' Подпрограмма вычисления определителя квадратной матрицы
' Для вычисления определителя квадратной матрицы она
' приводится к треугольному виду и находится произведение
' ее диагональных элементов. См. Мудров А.Е. Численные методы
' для ПВЭМ, Томск, 1991, 270с. ил.
' Стр. 51: Программа 2.3. Вычисление определителя по Гауссу
' N -
размерность квадратной матрицы;
' A() -
элементы исходной матрицы;
' DET -
величина определителя матрицы.
Dim As Integer I, J, K, K1, P
Dim As Single
R
P = 1
For K = 1 To N -
1
K1 = K + 1: DET = A(K,K): J = K
For I = K1 To N
R = A(I,K)
If Abs(R) > Abs(DET) Then
DET = R: J = I
EndIf
Next I
'
If DET = 0 Then Exit Sub
If J <> K Then For I = K To N
R = A(K,I)
A(K,I) = A
(J,I)
A(J,I) = R
Next I
P = -
P
EndIf
For J = K1 To N: A(K,J) = A(K,J)/DET: Next J
For I = K1 To N: R = A(I,K)
For J = K1 To N: A(I,J) = A(I,J) -
A(K,J)*R: Next J
Next I
20
P = P*DET
Next K
DET = P*A(N,N)
End Sub
'
' Основной модуль прог
раммы
Dim As Integer N, I, J
Data 3
' Матрица из примера программы 2.3. книги Мудрова
' ---
det(A) = -
52.000000
Data 0.0, 3.0, 1.0
Data 2.0, 8.0, 2.0
Data 4.0, -
1.0, 7.0
' Матрица из примера на сингулярное разложение (см. далее)
' --
-
det(A) = 0.000000
'Data 32.0, 14.0, 74.0
'Data -
24.0, -
10.0, -
57.0
'Data -
8.0, -
4.0, -
17.0
Read N
Dim As Single A(N, N), DET
For I = 1 To N
For J = 1 To N: Read A(I,J): Next J
Next I
'
MDET(N, A(), DET) ' вызов
подпрограммы
'
Open Cons For Ou
tput As #2
'Open "MDET.RES" For Output As #2
Print #2,: Print #2, Using " ---
det(A) = ###.######"; DET
Close #2
Sleep
'
Дело не в том, что программа выдает правильные результаты, а в том, что довольно часто встречаются системы, основные матрицы которых им
еют нулевой определитель. Определитель (или детерминант) —
одно из основных понятий линейной алгебры. Вырожденная система —
это система, описываемая матрицей с нулевым определителем (сингулярной матрицей). Поскольку некоторые уравнения, входящие в такую си
стему, представляются линейной комбинацией других уравнений, то фактически сама система является недоопределенной. Несложно сообразить, что, в зависимости от конкретного вида вектора правой части уравнений, существует либо бесконечное множество решений, ли
бо не существует ни одного.
Метод SVD начинается с составления входной матрицы A, получившей в статистическом анализе экспериментов название "матрица плана". Это прямоугольная матрица из m строк и n столбцов. Подпрограмма SVD, исходя из входной матрицы А,
получает на выходе три матрицы SIG, U и V. Матрица SIG диагональная с неотрицательными диагональными элементами, называемыми сингулярными числами матрицы А. Левая матрица U и правая матрица V могут использоваться для преобразования исходных уравнений в эк
вивалентные диагональные системы.
Сегодня уже забыто, что сингулярное разложение имеет довольно давнюю историю. Основополагающими в большой мере были работы Голуба и его коллег 21
Кахана, Бизингера и Райнша. Хорошим дополнением к книге Форсайт, Малькольм, Мо
улер, Машинные методы... является книга:
Уилкинсон Дж.Х., Райнш С. Справочник алгоритмов на языке Алгол. Линейная алгебра. М.: 1976. 390 с.
В книге приведены алгоритмы решения всех основных задач линейной алгебры, реализованные в виде процедур на языке Ал
гол
-
60. Для специалистов по теории управления представляют интерес алгоритмы решения проблемы собственных значений для произвольных матриц.
2.3.1.
SVD для решения линейных уравнений.
Ниже приведена программа, в которой подпрограмма MSVD используется дл
я решения системы алгебраических уравнений размера n = 3 и m =3 основная матрица которой изначально сингулярна: первая строка равна сумме второй и третьей с обратным знаком.
' P R O G R A M "Appr_MSVD1"
' 12.08.1994
'
-------------------------------------------------------------
' Решение линейных уравнений с нулевым определителем
'
-------------------------------------------------------------
'
' Ординарная переработка программы из книги:
' Форсай
т Дж., Малькольм М., Моулер К., Машинные методы
' математических вычислений. М., Мир, 1980, 280 с., илл.
' Алгоритм и пример: Форсайт, Малькольм, с.248
' Сингулярное
разложение
(singular value decomposition, SVD)
' это разложение прямоугольной матрицы.
..
' I = 1,...M -
номер строки (число элементов Y(I))
' J = 1,...N -
номер столбца (число элементов C(J))
' C(J) -
корни уравнений;
' COND -
число обусловленности;
' TAU -
показатель точности;
#Lang "FB" ' режим совместимости с FreeBASIC
' О
бъявление
используемой
подпрограммы
Declare Sub MSVD(M As Integer, N As Integer, A() As Single, SIG() As Single, U() As Single, V() As Single, IER As Integer)
'
' Основной модуль программы
Dim As Integer M, N, I, J, IER
Dim As Single A(10, 10), Y(10), Y0(
10), SIG(10), U(10, 10), V(10, 10)
Dim As Single SIGmax, SIGmin, TAU, COND, SIG1, S, C(10)
' Число строк M равно числу столбцов N матрицы A:
Data 3, 3
' Исходная матрица:
Data 32.0, 14.0, 74.0
Data -
24.0, -
10.0, -
57.0
Data -
8.0, -
4.0, -
17.0
22
' Векто
р правой части:
Data -
14.0, 13.0, 1.0
' Ввод элементов исходной матрицы
Read M, N
For I = 1 TO M
For J = 1 TO N: Read A(I,J): Next J
Next I
' Ввод вектора правой части
For I = 1 TO M: Read Y(I): Next I
'
Open Cons For Output As #2
'Open "Appr_MSVD.res
" For Output As #2
Print #2,: Print #2, " Singular Value Decomposition"
Print #2,: Print #2, " ---
Initial matrix A(I,J):"
Print #2,
For I = 1 To M
For J = 1 To N: Print #2, Using " ###.######"; A(I, J);: NEXT J
Print #2,
Next I
Print #2,: Print #2, " ---
Right parts of the equations Y(I):"
Print #2,
For I = 1 TO M
Print #2, Using " ###.######"; Y(I)
Next I
'
MSVD (M, N, A(), SIG(), U(), V(), IER)
'
If IER <> 0 Then Print #2, " SVD ---
> Error", IER
Print #2,: Print #2, " ---
Singular values SIG(J):"
Print #2,
For J = 1 TO N
Print #2, Using " ###.######"; SIG(J)
Next J
SIGmax = SIG(1): SIGmin = SIG(1)
For J = 2 To N
If SIG(J) > SIGmax Then SIGmax = SIG(J)
If SIG(J) < SIGmin And SIG(J) > 0 Then SIGmin = SIG(J)
Next J
COND = SIGmax/SIGmin
Print #2,: Pr
int #2, Using " ---
Indicator conditioning COND = #.#####^^^^"; COND
SIG1 = 0
For J = 1 To N
If SIG(J) > SIG1 Then SIG1 = SIG(J) C(J) = 0
TAU = SIG1*0.1E
-
6
Next J
For J = 1 To N
If SIG(J) <= TAU Then GoTo L60 S = 0
For I = 1 To M
23
S = S + U(I,J)*Y(I)
Next
I
S = S/SIG(J)
For I = 1 To N
C(I) = C(I)+S*V(I,J)
Next I
L60: Next J
Print #2,: Print #2, " ---
Roots of the system of equations C(J):"
Print #2,
For J = 1 TO N
Print #2, Using " ###.######"; C(J)
Next J
' Проверка вычисленных корней подстановкой в ура
внения
For I = 1 To M
S = 0
For J = 1 To N
S = S + A(I,J)*C(J)
Next J
Y0(I) = S
Next I
Print #2,: Print #2, " ---
Values of the right
-
hand sides Y0(I):"
Print #2,
For I = 1 TO M
Print #2, Using " ###.######"; Y0(I)
Next I
Print #2,: Print #2, " ---
Left m
atrix U(I,J):"
Print #2,
For I = 1 To M
For J = 1 To N: Print #2, Using " ###.######"; U(I,J);: Next J
Print #2,
Next I
Print #2,: Print #2, " ---
Right matrix V(I,J):"
Print #2,
For I = 1 To N
For J = 1 To N: Print #2, Using " ###.######"; V(I,J);: Ne
xt J
Print #2,
Next I
Close #2
Sleep
' Результат выполнения программы:
' Singular Value Decomposition
' ---
Initial matrix A(I,J):
' 32.000000 14.000000 74.000000
' -
24.000000 -
10.000000 -
57.000000
' -
8.000000 -
4.000000 -
17.000000
' ---
Right parts of the equations Y(I):
' -
14.000000
' 13.000000
24
' 1.000000
' ---
Singular values SIG(J):
' 104.825485
' 1.271749
' 0.000001
' ---
Indicator conditioning COND = 0.12085E+09
' ---
Roots of the system of equations C(J):
' 1.21539
4
' 1.821740
' -
1.059419
' ---
Values of the right
-
hand sides Y0(I):
' -
14.000019
' 13.000010
' 1.000005
' ---
Left matrix U(I,J):
' -
0.780617 0.239384 0.577350
' 0.597620 0.556343 0.577349
' 0.182997 -
0.795725 0.577351
' ---
Right matrix V(I,J):
' -
0.389090 0.529858 0.753564
' -
0.168249 0.763391 -
0.623640
' -
0.905705 -
0.369439 -
0.207880
'
Sub MSVD(M As Integer, N As Integer, A() As Single, SIG() As Single, U() As Single, V() As Single, IER As Integ
er)
' Подпрограмма сингулярного разложения матрицы А=U*SIG*V^Т
' A -
действительная матрица размера M*N (M >= N)
' M -
количество строк матрицы А
' N -
количество столбцов матрицы А
' SIG -
вектор сингулярных чисел
' RV1 -
вспомогательный масс
ив
' IER -
код завершения подпрограммы /0 -
нормальное/
' Внешних процедур нет!
'
Dim As Integer I, I1, J, K, K1, L, L1, MN, ITS
Dim As Single RV1(10), G, G1, S, SCALE, ANORM, F, H, C, X, Y, Z
'
' Подготовительные операции
IER = 0
For I = 1 To M
For J =
1 To N: U(I, J) = A(I, J): Next J
Next I
'
' Приведение к двухдиагональной форме
G = 0: SCALE = 0: ANORM = 0
'
For I = 1 To N
L = I + 1
RV1(I) = SCALE*G
G = 0: S = 0: SCALE = 0
If I > M Then GoTo L210
25
For K = I To M: SCALE = SCAL
E + ABS(U(K,I)): Next K
If SCALE = 0 Then GoTo L210
For K = I To M
U(K,I) = U(K,I)/SCALE
S = S + U(K,I)*U(K,I)
Next K
F = U(I,I)
G = ABS(SQR(S))
If F < 0 Then G = -
G
G = -
G
H = F*G -
S
U(I,I) = F -
G
If I = N THEN GOTO L190
For J = L TO N
S = 0
For K = I TO M
S = S + U(K,I)*U(K,J)
Next K
F = S/H
For K = I TO M
U(K,J) = U(K,J) + F*U(K,I)
Next K
Next J
L190:For K = I To M
U(K, I)
= SCALE*U(K,I)
Next K
L210:SIG(I) = SCALE * G
G = 0: S = 0: SCALE = 0
If I > M Or I = N Then GoTo L290
For K = L To N
SCALE = SCALE + Abs(U(I,K))
Next K
If SCALE = 0 Then GoTo L290
For K = L To N
U(I,K) = U(I
,K)/SCALE
S = S + U(I,K)*U(I,K)
Next K
F = U(I,L)
G = Abs(Sqr(S))
If F < 0 Then G = -
G
G = -
G
H = F*G -
S
U(I,L) = F -
G
For K = L To N
RV1(K) = U(I,K)/H
Next K
If I = M Then GoTo L270
For J = L TO M
S = 0
For K = L To N
S = S + U(J,K)*U(I,K)
Next K
26
For K = L To N
U(J,K) = U(J,K) + S*RV1(K)
Next K
Next J
L270:For K = L To N
U(I,K) = SCALE*U(I,K)
Next K
L290:If ANORM < Abs(SIG(I)) + Abs(RV1(I)) Then ANORM = Abs(SIG(I)) + ABS(RV1(I))
Next I
'
' Формирование матрицы правых преобразований
For I = N To 1 Step -
1
If I = N Then GoTo L390
If G = 0 Then GoTo L360
For J = L To N
V(J,I) = (U(I,J)/U(I,L))/G
Next J
For J = L To N
S = 0
For K = L To N
S = S + U(I,K)*V(K,J)
Next K
For K = L To N
V(K,J) = V(K,J) + S*V(K,I)
Next K
Next J
L360:For J = L TO N
V(I,J) = 0: V(J,I) = 0
Next J
L390:V(I,I) = 1
G = RV1(I): L = I
Next I
'
' Формирование матрицы левых преобразований
MN = N
If M < N THEN MN = M
For I = MN TO 1 STEP -
1
L = I + 1: G = SIG(I)
If I = N THEN GOTO L430
For J = L TO N
U(I, J) = 0
Next J
L430
:If G = 0 THEN GOTO L475
If I = MN THEN GOTO L460
For J = L TO N
S = 0
For K = L TO M
S = S + U(K, I) * U(K, J)
Next K
F = (S/U(I,I))/G
For K = I TO M
27
U(K,J) = U(K,J) + F*U(K,I)
Next K
Next J
L460:Fo
r J = I To M
U(J,I) = U(J,I)/G
Next J
GoTo L490
'
L475:For J = I To M
U(J,I) = 0
Next J
L490:U(I,I) = U(I,I) + 1
Next I
'
' Диагонализация двухдиагональной формы
For K = N To 1 Step -
1
K1 = K -
1: ITS = 0
L520:For L = K TO 1 STEP -
1
L1 = L -
1
If Abs(RV1(L)) + ANORM = ANORM Then GoTo L565
If Abs(SIG(L1)) + ANORM = ANORM Then GoTo L540
Next L
L540:C = 0: S = 1
For I = L TO K
F = S*RV1(I)
RV1(I) = C * RV1(I)
If Abs(F) + AN
ORM = ANORM Then GoTo L565
G = SIG(I)
H = Sqr(F*F + G*G)
SIG(I) = H: C = G/H: S = -
F/H
'
For J = 1 To M
Y = U(J,L1): Z = U(J,I)
U(J,L1) = Y*C + Z*S
U(J,I) = Z*C -
Y*S
Next J
Next I
'
' Проверка сходи
мости
L565:Z = SIG(K)
If L = K Then GoTo L650
'
' Формирование сдвига для QR
-
преобразования
If ITS = 30 Then IER = -
1: Exit For
ITS = ITS + 1
X = SIG(L): Y = SIG(K1)
G = RV1(K1): H = RV1(K)
F = ((Y -
Z)*(Y + Z) + (G -
H)*(G +
H))/(2*H*Y)
G = Sqr(F*F + 1)
G1 = G
If F < 0 Then G1 = -
G1
F = ((X -
Z)*(X + Z) + H*(Y/(F + G1)) -
H)/X
28
'
' Очередной шаг QR
-
преобразования
C = 1: S = 1
If L > K1 Then GoTo L610 ' ???
For I1 = L To K1
I = I1 + 1
G = RV1(I): Y = SIG(I)
H = S*G: G = C*G
Z = Sqr(F*F + H*H)
RV1(I1) = Z
C = F/Z: S = H/Z
F = X*C + G*S
G = G*C -
X*S
H = Y*S: Y = Y*C
'
For J = 1 To M
X = V(J,I1): Z = V(J,I)
V(J, I1) = X*
C + Z*S
V(J, I) = Z*C -
X*S
Next J
'
Z = Sqr(F*F + H*H)
SIG(I1) = Z
'
' Если Z равно 0 -
произвольное вращение
If Z <> 0 Then C = F/Z: S = H/Z
F = C*G + S*Y
X = C*Y -
S*G
'
For J = 1 To M
Y = U(J,I1)
Z = U(J,I)
U(J,I1) = Y*C + Z*S
U(J,I) = Z*C -
Y*S
Next J
'
Next I1
'
L610:RV1(L) = 0
RV1(K) = F: SIG(K) = X
GoTo L520
'
' Формирование массива из неотрицательных элементов
L650:If Z < 0 Then
SIG(K) = -
Z
For J = 1 To N
V(J,K) = -
V(J,K)
Next J
EndIf
Next K
End Sub
'
29
Как видно, практически без хлопот получено решение системы уравнений, при проверке дающее значения правых частей уравнений с приемлемой точностью:
Y(1) -
14.000019
Y(2) 13.000010
Y(
3) 1.000005
Конечно, текст подпрограммы MSVD следует сохранить –
подпрограмма будет использоваться в следующем разделе для решения другой задачи!
2.3.2.
SVD для решения задач аппроксимации.
Следующим примером использования подпрограммы MSVD из книг
и Форсайт, Малькольм, Моулер, Машинные методы... будет попытка экстраполяции статистических данных
из упражнения на странице 95, где приведена таблица данных Бюро переписей населения Соединенных Штатов Америки по восьми годам (m = 8):
Год
Население
1900
75 994 575
1910
91 972 266
1920
105 710 620
1930
123 203 000
1940
131 669 275
1950
150 697 361
1960
179 323 175
1970
203 211 926
Ниже приведена программа аппроксимирующая эти данные квадратичным многочленом вида:
Y(T) = С(1) + С(2)*T + С(3)*T^2,
и
использование многочлена для предсказания численности населения в 1980 году
.
Полученное "число обусловленности" равное COND = SIGmax/SIGmin = 0.30418E+11 сигнализирует о наличии каких
-
то трудностей. Но, "предупрежден, значит –
вооружен" и можно принять н
екоторые дополнительные меры. Приведенная программа проверялась при задании действительных чисел в форме с плавающей точкой как формата SINGLE, так и формата DOUBLE
. Большой разницы замечено не было.
' P R O G R A M "Appr_MSVD2"
' 14.08.1994
'
-------------------------------------------------------------
' Сингулярный анализ в задачах выравнивания (МНК)
'
-------------------------------------------------------------
'
' Ординарная переработка програ
ммы из книги:
' Форсайт Дж., Малькольм М., Моулер К., Машинные методы
' математических вычислений. М., Мир, 1980, 280 с., илл.
30
' Иллюстрирующий пример со страницы 218:
' Использование подпрограммы SVD для предсказания численности
' населения на данных
Бюро переписей США... ' Сингулярное
разложение
(singular value decomposition, SVD)
' это разложение прямоугольной матрицы...
' I = 1,...M -
номер строки (число элементов Y(I))
' J = 1,...N -
номер столбца (число элементов C(J))
' C(J) -
корни уравне
ний;
' COND -
число обусловленности;
' TAU -
показатель точности;
' TF -
прогнозируемый год;
' YF -
вычисленное значение численности;
'
#Lang "FB" ' режим совместимости с FreeBASIC
' Объявление
используемой
подпрограммы
Declare Sub MSVD(M
As Integer, N As Integer, A() As Single, SIG() As Single, U() As Single, V() As Single, IER As Integer)
'
' Основной модуль программы
Dim As Integer M, N, I, J, IER
Dim As Single T(10), TF, Y(10), YF, A(10, 10), Y0(10), SIG(10), U(10, 10), V(10, 10)
Dim As Single SIGmax, SIGmin, TAU, COND, SIG1, S, C(10)
' Число строк = M; число столбцов = N
Data 8, 3
' Исходные данные:
' Годы Население (млн.чел.)
Data 1900, 75.995
Data 1910, 91.972
Data 1920, 105.711
Data 1930, 123.203
Data 1940, 131.669
D
ata 1950, 150.697
Data 1960, 179.323
Data 1970, 203.212
Data 1980 ' дата
прогноза
' Ввод исходных данных и вычисление элементов матрицы
Read M, N
For I = 1 To M
Read T(I), Y(I)
A(I,1) = 1.0
For J = 2 To N
A(I,J) = T(I)*A(I,J
-
1)
Next J
Next I
' Ввод
даты
прогноза
Read TF
'
Open Cons For Output As #2
'Open "Appr_MSVD.res" For Output As #2
Print #2,: Print #2, " Singular Value Decomposition"
Print #2,: Print #2, " ---
Initial matrix A(I,J):"
31
Print #2,
For I = 1 To M
For J = 1 To N: Print #2, Using " ##.####^^^^"; A(I, J);: Next J
Print #2,
Next I
Print #2,: Print #2, " ---
Right parts of the equations Y(I):"
Print #2,
For I = 1 TO M
Print #2, Using " ##.####^^^^"; Y(I)
Next I
'
MSVD (M, N, A(), SIG(), U(), V(), IER)
'
If IER <> 0 Then Print #2, " SVD ---
> Error", IER
Print #2,: Print #2, " ---
Singular values SIG(J):"
Print #2,
For J = 1 TO N
Print #2, Using " ##.####^^^^"; SIG(J)
Next J
SIGmax = SIG(1): SIGmin = SIG(1)
For J = 2 To N
If SIG(J) > SIGmax Then SIGmax = SIG(J)
If SIG(J) < SIGmin And
SIG(J) > 0 Then SIGmin = SIG(J)
Next J
COND = SIGmax/SIGmin
Print #2,: Print #2, Using " ---
Indicator conditioning COND = #.#####^^^^"; COND
SIG1 = 0
For J = 1 To N
If SIG(J) > SIG1 Then SIG1 = SIG(J) C(J) = 0
TAU = SIG1*0.1E
-
6
Next J
For J = 1 To N
If SIG(J) <= TAU Then GoTo L60 S = 0
For I = 1 To M
S = S + U(I,J)*Y(I)
Next I
S = S/SIG(J)
For I = 1 To N
C(I) = C(I)+S*V(I,J)
Next I
L60: Next J
Print #2,: Print #2, " ---
Roots of the system of equations C(J):"
Print #2,
For J = 1 TO N
Print #2, Using " ##.####^^^^"; C(J)
Next J
'
32
' Проверка вычисленных корней подстановкой в уравнения
For I = 1 To M
S = 0
For J = 1 To N
S = S + A(I,J)*C(J)
Next J
Y0(I) = S
Next I
Print #2,: Print #2, " ---
Values of the right
-
hand sides Y0(I):"
Print #2,
For I = 1 TO M
Print #2, Using " ##.####^^^^"; Y0(I)
Next I
' Прогноз численности населения в 1980 году:
Print #2,: Print #2, " ---
Forecast of the number of population:"
YF = C(1) + C(2)*TF + C(3)*TF^2
Print #2,
Print #2, Using " TF = ##### YF = ##.####^^^^"; TF;
YF
' Вывод на печать левой и правой матриц:
Print #2,: Print #2, " ---
Left matrix U(I,J):"
Print #2,
For I = 1 To M
For J = 1 To N: Print #2, Using " ###.######"; U(I, J);: Next J
Print #2,
Next I
Print #2,: Print #2, " ---
Right matrix V(I,J):"
Prin
t #2,
For I = 1 To N
For J = 1 To N: Print #2, Using " ###.######"; V(I, J);: Next J
Print #2,
Next I
Close #2
Sleep
' Результат выполнения программы -
вариант SINGLE:
' Singular Value Decomposition
' ---
Initial matrix A(I,J):
' 1.0000E+00 1.90
00E+03 3.6100E+06
' 1.0000E+00 1.9100E+03 3.6481E+06
' 1.0000E+00 1.9200E+03 3.6864E+06
' 1.0000E+00 1.9300E+03 3.7249E+06
' 1.0000E+00 1.9400E+03 3.7636E+06
' 1.0000E+00 1.9500E+03 3.8025E+06
' 1.0000E+00 1.9600E+03 3
.8416E+06
' 1.0000E+00 1.9700E+03 3.8809E+06
' ---
Right parts of the equations Y(I):
' 7.5995E+01
' 9.1972E+01
' 1.0571E+02
33
' 1.2320E+02
' 1.3167E+02
' 1.5070E+02
' 1.7932E+02
' 2.0321E+02
' ---
Singular values SIG(J):
' 1.0595E+0
7
' 6.4777E+01
' 3.4831E
-
04
' ---
Indicator conditioning COND = 0.30418E+11
' ---
Roots of the system of equations C(J):
' -
1.6705E
-
03
' -
1.6161E+00
' 8.7048E
-
04
' ---
Values of the right
-
hand sides Y0(I):
' 7.1748E+01
' 8.8751E+01
' 1.0593E+
02
' 1.2328E+02
' 1.4081E+02
' 1.5851E+02
' 1.7638E+02
' 1.9443E+02
' ---
Forecast of the number of population:
' TF = 1980 YF = 2.1265E+02
' ---
Left matrix U(I,J):
' -
0.340736 0.542688 -
0.538995
' -
0.344332 0.393157 -
0.089
139
' -
0.347947 0.242132 0.226345
' -
0.351581 0.089495 0.384927
' -
0.355233 -
0.064771 0.383711
' -
0.358905 -
0.220530 0.248785
' -
0.362596 -
0.378099 -
0.080733
' -
0.366305 -
0.537115 -
0.535249
' ---
Right matrix V(I,
J):
' -
0.000000 0.001034 -
1.000000
' -
0.000517 0.999999 0.001034
' -
1.000000 -
0.000517 -
0.000000
'
' Результат выполнения программы -
вариант DOUBLE:
' ---
Singular values SIG(J):
' 1.0595E+07
' 6.4775E+01
' 3.4620E
-
04
' ---
In
dicator conditioning COND = 0.30603E+11
' ---
Roots of the system of equations C(J):
' -
1.6706E
-
03
' -
1.6162E+00
' 8.7056E
-
04
' ---
Values of the right
-
hand sides Y0(I):
' 7.1980E+01
' 8.8987E+01
34
' 1.0617E+02
' 1.2352E+02
' 1.4105E+02
' 1.5
875E+02
' 1.7663E+02
' 1.9468E+02
' ---
Forecast of the number of population:
' TF = 1980 YF = 2.1291E+02
'
2.4.
Метод итерационного R
-
алгоритма.
И на десерт любопытная программа из Библиотеки численного анализа НИВЦ МГУ:
http://num
-
anal.srcc.msu.ru/lib_na/cat/i/ia83r.htm
Подпрограмма IA83R -
построение алгебраического полинома наилучшего равномерного приближения для таблично заданной функции.
Поскольку указанный на сайте университетский сборник практически не доступен, возможно заменить его книгой:
Демьянов В.Ф., Малоземов В.Н. Введение в минимакс. М.: Наука, 1972. 368 с.
Минимакс (минимизация максимального уклонения) —
принцип оптимального выбора параметров. В первых двух
главах книги рассматривается простейшая (и исторически первая) линейная минимаксная задача —
построение алгебраического полинома наилучшего приближения. В остальных четырех главах развивается общая теория нелинейных минимаксных задач. Отдельно рассматрива
ются дискретный и непрерывный случаи, отсутствие и наличие ограничений на параметры. Основные вопросы: дифференцируемость функции максимума по направлениям, необходимые условия минимакса, достаточные условия локального минимакса, методы последовательных пр
иближений для нахождения стационарных точек. Большое количество примеров и рисунков иллюстрируют основные результаты теории.
В представленной ниже задаче на отрезке [
-
1, 1] ищется полином наилучшего приближения 7
-
ой степени вида:
P7(T) = A(1) + A(2)*T + ... + A(N)*T^(N
-
1) + A(N+1)*T^N
' P R O G R A M "Appr_IA83R"
' 10.10.2012
'
-------------------------------------------------------------
' Аппроксимация рациональными функциями (равномерная)
'
----
---------------------------------------------------------
'
' Построение алгебраического полинома наилучшего равномерного
' приближения для таблично заданной функции.
' См. Демьянов В.Ф., Малоземов В.Н. Введение в минимакс.1972.
' Поиск полинома наилуч
шего приближения 7
-
ой степени...
35
' M -
число узлов сетки;
' N -
степень полинома.
#Lang "FB" ' режим FreeBASIC
-
совместимости
'
Sub UTIA10(IERR As Integer, N As Integer)
' Подпрограмма выдачи сообщений при работе IA83R
' IF N = 83 Then GoT
o L83
If IERR = 65 Then Print " Er65: Greater the degree polynomial"
If IERR = 66 Then Print " Er66: Set points not been completed"
If IERR = 67 Then Print " Er67: Set points not sorted order"
If IERR = 68 Then Print " Er68: Exhausted limit of iterations"
End Sub
'
Sub IA83R(M As Integer, N As Integer, X() As Single, Y() As Single, IN() As Integer, ByRef ITER As Integer, A() As Single, AX() As Single, JN() As Integer, RAB() As Single, IERR As Integer)
Dim As Integer I,J,J1,J2,K,IZ,R1,R2,K1,K2,N1,N2,INDEX
Di
m As Single RHO,ZERO,H,Z,W,RMAX,RMIN,RMAX1,RMIN1,U
Dim As Integer D
ZERO = 0: IERR = 0: N1 = N+1: ITER = 0: N2 = N+2
L30:RHO = ZERO
If N2 > M Then GoTo L25
L1: J = -
1
ITER = ITER+1
J2 = IN(1)
For I = 1 To N2
J1=J2
J2=IN(I)
AX(I)
=X(J2)
RAB(1,I)=Y(J2)
J=
-
J
RAB(2,I)=J
If I = 1 Then GoTo L2
If J1 >= J2 Then GoTo L26
If AX(I) <= AX(I
-
1) Then GoTo L27
L2: Next I
For K = 2 To N2
For J = K To N2
I=N2
-
J+K
IZ=I
-
K+1
U=AX(I)
-
AX(IZ)
RAB(1,I)=(R
AB(1,I)
-
RAB(1,I
-
1))/U
RAB(2,I)=(RAB(2,I)
-
RAB(2,I
-
1))/U
Next J
Next K
IZ=1
H=RAB(1,N2)/RAB(2,N2)
For I = 1 To N1
A(I)=RAB(1,I)
-
H*RAB(2,I)
Next I
For J = 1 To N
36
I=N1
-
J
U=AX(I)
For K = I To N
A(K)=A(K)
-
A(K+1
)*U
Next K
Next J
U=RHO
RHO=ABS(H)
If RHO > U Then GoTo L6
RAB(1,1) = RHO
GoTo L29
L6: If ITER = 30 Then GoTo L31
IN(1)=1
IN(N2)=M
D = 1 ' .TRUE.
Z=ZERO
RMAX1=H/2.
RMIN1=
-
RMAX1
For I = 1 To N1
RMIN=ZERO
RMAX=ZERO
J1=IN(I)
R1=J1
R2=J1
J2=IN(I+1)
For J = J1 To J2
W=ZERO
U=X(J)
For K1 = 1 To N1
K=N1+1
-
K1
W=W*U+A(K)
L7: Next K1
W=W
-
Y(J)
If W <= RMAX Then GoTo L9
L8: RMAX=W
R1=J
L9: If W >= R
MIN Then GoTo L11
L10:RMIN=W
R2=J
L11:Next J
If R1 <= R2 Then GoTo L13
L12:J=R1
R1=R2
R2=J
U=RMAX
RMAX=RMIN
RMIN=U
L13:
K1 = 1: If RMIN1 < 0 Then K1 = -
1
K2 = 1: If RMAX < 0 Then K2 = -
1
If K1 <> K2 Then GoTo L18
37
L14
:If ABS(RMIN1) >= ABS(RMAX) Then GoTo L15
JN(I)=R1
JN(I+1)=R2
RMAX1=RMAX
RMIN1=RMIN
If ABS(Z) > ABS(RMAX1) Then JN(I
-
1)=IZ
If D < 0 Then GoTo L17
H=RMAX
D = -
1
GoTo L17
L15: RMAX1=RMIN1
If ABS(Z) <
= ABS(RMIN) Then GoTo L16
JN(I+1)=IZ
RMIN1=Z
GoTo L17
L16: JN(I+1)=R2
RMIN1=RMIN
L17: Z=ZERO
GoTo L21
L18: RMAX1=RMIN1
If Abs(Z) <= Abs(RMAX) Then GoTo L19
JN(I+1)=IZ
RMIN1=Z
GoTo L20
L19: JN(I+1)=R1
RMIN1=
RMAX
L20: IZ=R2
Z=RMIN
L21: Next I
If ABS(Z) <= Abs(H) Then GoTo L23
For J = 1 To N1
JN(J)=JN(J+1)
L22: Next J
JN(N2)=IZ
L23: INDEX=0
For I = 1 To N2
IN(I)=JN(I)
If IN(I) >= IN(I
-
1) Or I = 1 Then GoTo L24
INDEX=
IN(I)
IN(I)=IN(I
-
1)
IN(I
-
1)=INDEX
L24: Next I
If INDEX > 0 Then GoTo L30
GoTo L1
L25: IERR=65
GoTo L28
L26: IERR=66
GoTo L28
L27: IERR=67
GoTo L28
L31: IERR=68
L28: 38
IF IERR <> 0 Then UTIA10(IERR,83)
L29:
End Sub
'
' Осн
овной модуль программы
' M -
число узлов сетки;
' N -
степень полинома, N <= M
-
2;
' IN -
вектор длины N+2 номеров начального набора точек;
' JN -
вектор длины N+2 номеров вычисленных значений точек альтернанса;
' ITER -
число выполненных при
решении задачи итераций;
' A -
вектор N+1 вычисленных значений коэффициентов полинома;
' AX -
вектор длины N+2 вычисленных значений точек альтернанса;
' H -
отклонение полинома наилучшего приближения.
Dim As Integer M, N, I, K, IN(9), JN(9), IE
RR, ITER
Dim As Single X(51), Y(51), H, A(8), AX(9), RAB(2,9)
Data 1, 6, 11, 16, 21, 26, 31, 36, 51
M = 51: N = 7: H = 2/(M
-
1)
Open Cons For Output As #2
'Open " Appr_IA83R.res" For Output As #2
For I = 1 To N+2: Read IN(I): Next I
For K = 1 To M
X(K) = -
1
+ (K
-
1)*H
Y(K) = Abs(X(K))
Next K
IA83R (M, N, X(), Y(), IN(), ITER, A(), AX(), JN(), RAB(), IERR)
Print #2, Using " IERR = #### ITER = ####"; IERR; ITER
Print #2, Using " RAB(1,1) = ####.######"; RAB(1,1)
Print #2,: Print #2, " Calculated values of coef
ficients"
For I = 1 To N+1
Print #2, Using " ### ####.######"; I; A(I)
Next I
Print #2,: Print #2, " Value of alternation points"
For I = 1 To N+2
Print #2, Using " ### ####.######"; I; AX(I)
Next I
Print #2,: Print #2, " Rooms of the calculated values o
f alternation points"
For I = 1 To N+2
Print #2, Using " ####"; JN(I)
Next I
Close #2
Sleep
' Результат выполнения программы:
' IERR = 0 ITER = 6
' RAB(1,1) = 0.045876
' Calculated values of coefficients
' 1 0.045876
' 2 0.000001
39
' 3 2.869770
' 4 -
0.000001
' 5 -
4.181892
' 6 -
0.000002
' 7 2.312120
' 8 0.000000
' Value of alternation points
' 1 -
1.000000
' 2 -
0.880000
' 3 -
0.560000
' 4 -
0.200000
' 5 -
0.000000
' 6 0.200000
' 7 0.560000
' 8 0.880000
' 9 1.000000
' Rooms of the calculated values of alternation points
' 1, 4, 12, 21, 26, 31, 40, 48, 51
'
Следует отметить существование двух алгоритмов советского математика Евгения Яковлевича Ремеза
(1896 -
1975
): алгоритм для случая полиномиальной аппроксимации (он известен как второй алгоритм Ремеза) и первый итеративный алгоритм, менее удобный для ручного счета. На этих словах, как и ранее, очень трудно поставить точку. Думаю, на дополнительных занятиях препод
аватель продолжит эту тему... "Недостаточно, если ученик усвоит теорию, необходимо, чтобы ученик этой теорией овладел, а этого можно достигнуть только ее приложениями к практике и решением многочисленных задач и упражнений. Ученик должен иметь возможность видеть приложение элементов высшей математики, что достигается введением преподавания элементов механики и оптики".
Эти слова принадлежат П.Л.Чебышёву...
Как всегда расстаюсь не без сожаления
-
материал этой главы требует более детального, как минимум –
тщательного, рассмотрения. Как любит говорить моя хорошая знакомая: небрежность в подаче материала почти всегда возбуждает интерес к теме!
Будем надеяться, что это так (хотя знакомая имела ввиду совсем другой случай).
При этом конечно не должны быть забыт
ы такие важные для инженеров вопросы:
аппроксимация функции рядом Фурье
-
Лагерра;
прямое и обратное преобразования Лапласа.
Но сегодня лимит времени и объема брошюры исчерпан!
На этом позвольте закруглиться.
До новых встреч!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
145
Размер файла
630 Кб
Теги
программа, freebasic, абстрагирование, аппроксимация
1/--страниц
Пожаловаться на содержимое документа