Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 37
Текст из файла (страница 37)
выравнивание данных 2!5 Следующие далее сегменты фортранной программы иллюстрируют использование Я(Р для задач выравнивания по методу наименьших квадратов. Мы не привели полную подпрограмму, чтобы эти сегменты не потонули во множестве прочих деталей, таких, как ввод данных, генерирование базисных функций, вычисление модели в других точках, используя найденные коэффициенты, и т. д.
Матрица плана может быть построена следующим образом: РО 20 1 = 1, М Т(1) = адсиисса 1-ой заданной точки У(1) =ордината 1-ои заданнон точки РО101=1,Х А(1,1) = 1-я Базисная Чзонкннв, 10 СОХТ(Х()Е вычисленная а Т(!) 20 СОХТ1ХОЕ. В случае полиномиальной аппроксимации эффективный способ построения А таков: РО 20 1 = 1, М Т(1) =... УО) =-... А(1,1) = 1.0' РО!01 =. 2, Х А(1,1) = Т(1)чА(1,1-1) 10 СОХТ1Х()Е 20 СОХТ! Х() Е. Вслед за этим можно немедленно обратиться к ЯгР. (Относительно деталей использования ЯГР см. комментарии в ее распечатке.) СНЕ ЗУР(МР1М, М, Х, А, 31ОМА,,ТКРЕч (3, .ТК()Е., У, 1ЕКК, )УОКК) 1Р (1ЕКК,ХЕ.
О) (УК1ТЕ(6,13) !3 РОКМАТ('ОШИБКА. ВОЗВРАТ МЗ БУР') В нашей практике не было случая, чтобы при правильном использовании Я(Р произошел выход с диагностикой ошибки. Следующий сегмент программы находит на ибо льп!ее сингулярное число, которое будет участвовать в выявлении пренебрежимо малых сингулярных чисел. Здесь также присваивается нулевое значение вектору коэффициентов, что будет конечным результатом, если все сингулярные числа окажутся пренебрежимо малыми. 9.
НАИМЕНЬШИЕ КВАДРАТЫ 2!б Я!СМА! =.0. ХЗО ХО Х =- 1, ХЧ 1К (Б!СМА(Х),СТ. Я!СМА!) Б!СМА! = Б1СМА(Х) С(Х) = О. ЗО СОХЧТ!Ы(ХЕ Теперь нужно задать границу относительной погрешности КЕЕЕКК. Например, если данные имеют три верных знака", то подходящим будет значение КЕЕЕКК=10-з. Если же данные рассматриваются как точные, то КЕЕЕКК должно отражать ошибку округлений, ожидаемую при вычислениях внутри самой ЗЧР; эта ошибка имеет порядок умноженного на п машинного эпсилон. Граница абсолютной ошибки т получается теперь как ТА(Х КЕЬЕККеБ1СМА1. Следующий сегмент применяет к исходным данным преобразование (Х, находит преобразованные коэффициенты, а затем применяет преобразование (т, чтобы получить сами коэффициенты ТХО бо Х = 1, Ы 1Е (ХИСМА(Х) .ЬЕ.
ТАЫ) СО ТО бо Я=о. ХХО 40 1 = 1, М Я = Б + (Х(Х,Х)еу(1) 40 СОХЧТТЖХБ Я = Я/БХСМА(Х) ХЗО 501 1, ХЧ С(1) = С(!) + Бе'т'(1,1) Ю СОХЧТ1)Ч(ХБ бо СОХЧТИЧУЕ Заметим, что М вЂ” это число элементов в У и число строк в А и У, в то время как Ах — число элементов в С и число строк в (х. Все матрицы имеют Лх столбцов. Коэффициенты теперь готовы для использования. При желании их можно вывести на печать. ЧУКХТЕ(6,...)ХС(Х), 1=1, ХЧ) Вычисляется модель в произвольной точке ТТ таким образом: УУ =О. ХХО 70 Х = 1, ХЧ УУ = чХУ+ С(Х)Е(Х я баЗиСНая Чаиияция, вычисленная в ТТ) 70 СОХЧТХМ(ХЕ ') Имеется в виду десятичных.— Прим. нерее, эл.
вывлвнивлиие длнных 2|7 Для полиномиальной модели следует использовать схему Гор- пера (см. 6 4.2). ь'т = о. 2)О 70 1В = 1, Э( 1=)я+1-1В ТУ = ТТ.УТ+ С(1) 70 СО|ЧТ1ХОЕ Без труда можно получить квадратный корень из суммы квад- ратов невязок, т. е. ту величину, которая минимизируется. КБО = О. РО 90 1 = 1, М К1 = О. 1)О 80 1 = 1, )Ч К1 = К1 + С(1)*А(1,1) ао сомт1 ше КБО = КБО + (К1 — У(1))э*2 90 СО)ЧПХОЕ К = БОКТ(КБО) Наконец, было сказано, что важной особенностью ЯЧ0 является способность обнаруживать зависимость и отсутствие единственности. Следующий фрагмент находит пренебрежимо малые сингулярные числа.
Для каждого найденного печатается список коэффициентов, называемых нулееыми коэффициентами. Любое кратное а этих коэффициентов можно добавить к напечатанным прежде, не увеличивая г более чем на ат. Коэффициенты нормируютси так, что сумма их квадратов равна 1.0. ЮХ|О()Е = .ТКОЕ. 1)О 1001 1, Х 1Р (Б1ОМА(1) .ОТ. ТА()) ОО ТО 100 1ЛЧ1О()Е = .РА13Е. ЭУК1ТЕ(6,10! ) ЮК1ТЕ(6,...) (Ч(1,1), 1 = 1, |Ч) 100 СО|чТ1Ы()Е 1Р (Пч|О()Е) 17К1ТЕ(6,102) 101 РОКМАТ( НУЛЕВЫЕ КО ЭФФИ(1ИЕНТЪ| ) 102 ЕОКМАТ(КОЗФФИ|АИЕНТЪ| ОПРЕДЕЛЕНЫ ОДНОЗНАЧНО ') Этот тип дополнительной информации для задач на наименьшие квадраты называется сингулярным анализом. Он может быть очень полезен при анализе сложных математических моделей. я пхименьшие квлдгхты з!з Детальное обсуждение сингулярного анализа имеется в книге Лоусон, Хэнсан (1974).
Использование сингулярного анализа для задач выравнивания наименьшими квадратами можно проиллюстрировать на данных Бюро переписей США из упр. 4.5. Заданных точек здесь восемгс т=б. Значения г; суть 1900, 19!О,..., 1970, а соответствующие значения у,. (единицей является миллион людей) примерно 75.99„91.97, ..., 203.2!. Если аппроксимировать эти данные квадратичным многачленом у(!) с,+с,!+с,Р, а вычисления проводить на 1ВМ 360 с удвоенной точностью, то будут получены приблизительно такие результаты с,=0.373х10"', с,=- — 0,402х10', с,=0.108х10 '. Проделав те же вычисления с обычной точностью, получим при- близительно следующее: с,= — 0.372х!О', с.,= — 0.368х!0', с,= — 0.905х10 '.
У этих двух наборов коэффициентов пе совпадают даже знаки. Если такую модель использовать для предсказания численности населения в 1980 г., то для коэффициентов, вычисленных с удвоенной точностью, прогноз будет 227.78 миллиона, в та время как для коэффициентов, полученных с обычной точностью, 145.21 миллиона. Ясно, что второй набор коэффициентов бесполезен. Но что можно сказать о коэффициентах первого набора? Насколько они достоверны? Матрица плана для этого примера имеет размеры 8хЗ: ! 0.190 х !О' 0.36100 х !О' 1 О.!91 Х 10' 0.3648! Х 10' 1 0.197 х !О' 0.39204 х 1О" Ее сингулярные числа можно с достаточной точностью найти в обоих случаях.
Они таковы: а,=0.10бх10", аз=0.648х!О', а,=0.346х10 '. Число обусловленности равно а,(о,=-0.30бх10", и столь большая величина сигнализирует наличие каких-то трудностей. Для 9.Е ВЫРАВНИВАНИЕ ДАННЫХ 2!9 значений ! между 1900 и 1970 гг. три базисные функции 1, Г и 1Р очень близки к линейной зависимости.
Чтобы исправить ситуацию, можно предпринять одну нз следующих двух мер. Во-первых, можно выбрать границу для относительной ошибки, которая бы отражала точность данных и точность арифметики. Если взять границу, например, в интервале от 0.3х!О " до О.бх[0-', то описанная выше программа сингулярного анализа отбросит и,. Тогда будут получены такие наборы коэффициентов: с,= — 0.167к!О ', с.,= — 0.162х!0', с,=-0.87!х10 ' для удвоенной точности и с,= — 0.166м 10-', с,= — 0.162к 10', с,=0.869х 10-а для обычной. Теперь они находятся в гораздо лучшем согласии друг с другом. Кроме того, коэффициенты стали существенно меньше, а это значит, что не будет столь большого, как прежде, взаимного уничтожения слагаемых при вычислении квадратичного многочлена.
Прогноз численности населения на 1980 г. теперь будет соответственно 212.91 и 214.96 миллиона. Эффект обычной точности все еще заметен, однако результаты уже не являются катастрофическими. Программа сингулярного анализа печатает также набор нулевых коэффициентов, соответствующих пренебрежимо малому сингулярному числу.
Вот эти коэффициенты: у,= — 0.9999995, 7,=0.103 х !О--', у, — — — 0.266х 10-'. Для значений ! от 1900 до 1970 величина функции [у,+у,!+7,12! не превосходит 0.0017. Поэтому при любом и коэффициенты можно изменить с с, на с!+ау;, и при этом значения, выдаваемые моделью, изменятся не более чем на 0.0017а. Любой из четырех вычисленных нами наборов коэффициентов можно получить из любого другого подобным изменением. Другая мера, которую можно принять для улучшения ситуации, состоит в изменении базиса. Модели вида у (!)жс„+с, (! — 1900)+с, (! — 1900)' или гораздо более удовлетворительны. Важно прн этом то, что независимое переменное преобразуется из интервала [1900, 19701 в какой-нибудь более приемлемый интервал вроде [О, 70! или, еще лучше, [ — 3.5, 3.6!. Сингулярные числа матриц планов, по- 9, нхнмвньшие квадухты лученных для этих двух моделей, будут сг,=0.884 х 10', а,=0.293 х 10'"', о,=0.119х 10~ и а, = О.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
