Сплайн. Аппроксимация данных — Архив WASM.RU
Вводная часть.
В данном документе речь пойдет об интерпретации экспериментальных данных.
Под экспериментальными данными понимается некий набор измерений, являющийся результатом опыта.
Так, например вы стоите на выходе из метро и считаете число людей, проходящих мимо вас за каждую минуту. Допустим, в результате у вас получились следующие измерения: 45,23,55,87,53,48,61... Обобщив измерения в виде гистограммы, вы получаете некоторое распределение:
|Количество измерений данного числа людей в минуту |30 ---- |29 | | |28 -- -- |27 | | |26 -- | |25 | -- |24 | | |23 | -- | -- | Число людей за минуту -----------------------------------> 44 45 46 47 48 49 50 51 52Для экспериментатора может быть интересно понять, не может ли данное распределение быть описано некоторой аналитической формулой. Иначе говоря, определить природу процесса: с какой вероятностью можно ожидать выхода определенного числа людей в произвольном измерении, среднее число людей, выходящих из метро, наиболее вероятное и т.д.
У вас может возникнуть вопрос, зачем такие сложности, если все эти вопросы можно снять прямым измерением. Т.е. пусть перед нами стоит задача определить, сколько людей проходит за 8 часов через выход за рабочий день. Да, можно стоять все эти восемь часов и старательно записывать каждого отдельного человека. Именно так и поступает наше гос-во, когда хочет переписать всех мемберов нашей большой тусовки. При этом на исходе восьмого часа у вас уже слезятся глаза, у вас трижды проверяли документы и дважды подходили какие-то амбалы...
Однако анализируя результаты получасового опыта, в течении которого вы замеряли число человек за минуту, вы имеете другую возможность решить поставленную задачу. Допустим, вы предполагаете, что данное распределение (вероятность выхода определенного числа людей в минуту) описывается определенной формулой, проверяете свою гипотезу, убеждаетесь, что она более-менее достоверна. Формула дает Вам среднее значение распределения. Вам стало известно среднее число людей в минуту - достаточно умножить его на число минут в восьми часах для получения результата. Более того, вы также можете оценить ошибку ваших расчетов, поскольку подразумевается кратковременное (неполное) измерение, очевидно содержащее некоторую неопределенность результата в сравнении с абсолютно точным измерением, т.е. формула также может дать точность оценки, например:
Число человек в минуту = 47 ± 5;где 5 - точность оценки за минуту, или примерно 10%. С такой же точностью можно будет определить число людей за восемь часов:
Число человек за 8 часов = (47 ± 5)×(8×60) = 22560 ± 2400;Пример аппроксимации данных опыта.
Перейдем от абстрактного опыта к более жизненному. Предположим, вы анализируете свойства генератора 8-ми байтных ключей. Насколько случайными получаются ключи? Т.е. насколько близка псевдослучайная последовательность к истинно случайной.
Критерием качества псевдослучайной последовательности может служить число нулевых/единичных бит, пар бит 00-01-10-11, троек бит... Все они должны быть равновероятны, т.е. вероятность выпадения в каждом ключе определенного числа нулевых бит должно описывается так называемым биномиальным распределением:
k ƒ(k) = C / 2n; nгде n - полное число бит в ключе (в нашем случае 8×8 = 64), k - определенное число бит, Ckn - число вариантов размещения k нулевых бит по n свободным местам:
k C = n! / k! × (n-k)!; n2n - число всех битовых комбинаций в числе из n бит, ƒ(k) же вероятность выпадения k бит в числе из n бит.
Таким образом, появляется возможность проверки генератора ключей на уровне числа битов, пар битов, троек... Для этого необходимо получить достаточное число ключей (провести опыт), получить данные опыта - число нулевых бит в каждом ключе и попытаться аппроксимировать, т.е. описать полученные данные некоторой зависимостью (формулой). Если распределение числа нулевых бит недостоверно описывается биномиальным распределением, что будет наблюдается, например, в случае следующего алгоритма генерации ключа:
LOOP INDEX=1 TO 8 KEY[INDEX] = PASSWORD[INDEX] XOR MAGIC_NUMBER END LOOPто можно будет сделать определенные выводы о свойствах генератора.
В некоторых случаях аппроксимацию данных называют "сплайном". Данное наименование можно встретить в таких популярных программах, как Excel, в котором можно попытаться выполнить сглаживание (splain) зависимости, например полиномом степени от 2 до 6, правда весьма примитивно. Также можно упомянуть визуальное отображение звука в некоторых CD/mp3 плеерах в виде огибающей кривой. В качестве более-менее профессиональных средств анализа численных данных рекомендую Origin.
Еще один пример аппроксимации (реальный опыт).
Приведу также пример одного реального опыта, выполненного учеными во льдах Антарктики.
Как известно, полярная шапка нашей планеты представляет собой огромную толщу льда (до нескольких километров), лежащую на твердой (каменистой) поверхности Земли. Она образовывалась тысячи лет: вода замерзала и увеличивала эту толщу.
Ученые провели следующий опыт: они бурили небольшую скважину в толще льда и исследовали особенности намерзания (вспоминается файл "Лёд" из "Секретных материалов"). Оказалось, что намерзание происходило этапами: например, на некоторой глубине обнаружилось намерзание толщиной 10 метров, потом 50, потом 30 метров - подобно годичным кольцам деревьев.
Как же ученые интерпретировали полученные ими результаты толщин слоев намерзшего льда? Они сделали следующее достаточно простое предположение: допустим, толщина намерзания льда в зависимости от времени описывается суммой гармоник (синусоид), т.е.:
Намерзло льда (t) = S(ai×sin(wi×t + bi));где S означает сумму синусоид, ai - амплитуда i-ой синусоиды, wi - ее период, bi - ее начальная фаза.
Проведя соответствующие расчеты (называемые разложением в ряд Фурье, когда некоторую зависимость ƒ(x) представляют как сумму синусоид с различными частотами и амплитудами), ученые определили, что процесс намерзания хорошо описывается суммой всего трех гармоник с тремя различными периодами ! Это означает, что на нашей планете происходили (и происходят) некие периодические процессы, влияющие в том числе и на намерзание льда на полярных шапках. Так, один из периодов составил около 10 000 лет - следовательно, с такой периодичностью на Земле наступает некоторый процесс похолодания.
Приближение экспериментальных данных аналитической зависимостью.
Критерий Хи-квадрат.
Как же осуществить проверку достоверности выбранной аппроксимирующей функции? Допустим, у вас есть следующая зависимость:
X Y 1 2 2 4 3 6 4 7 5 10 6 13 ...Разумно полагая, что это что-то вроде y=ax, причем скорее всего a=1, только точки x=4 и x=6 не совсем вписываются в эту гипотезу, вы хотели бы проверить свою теорию и получить что-то вроде примерного значения "a" и некоторую оценку достоверности своего предположения.
Каким образом можно осуществить такую проверку? Сформулируем задачу в наиболее общем виде: обозначим аргументы (входные данные) за Xi, значения зависимости (выходные данные) за Yi, погрешность значений за Eyi. Аппроксимирующую функцию (допускаемый аналитический вид зависимости) за YiТеор. Тогда назовем критерием достоверности аппроксимации экспериментальной зависимости Yi = function (Xi) теоретической зависимостью YiТеор = function (Xi) функцию:
Cr = function(Xi,(Yi,Eyi),YiТеор);значение которой говорит о том, насколько точно хорошо теоретическое распределение YiТеор = function (Xi) описывает реальное распределение Yi = function (Xi).
В нашем простом случае у входных данных (X=1,2,3,4...) нет никаких погрешностей, равно как нет их и у выходных данных (Y=2,4,6,7,10,13). Теоретическое распределение имеет простой вид YТеор = aX с одним неизвестным параметром (a). Предполагаем a=2 и получаем ряд значений теоретического распределения:
YТеор(X) = 2X = 2,4,6,8,10,12Полученные значения YТеор необходимо подставить в критерий-функцию Cr и по ее значению оценить достоверность выбранной аппроксимации.
Какой же вид должна иметь критерий-функция Cr? Допустим, чем меньше ее значение, тем достовернее аппроксимация. Очевидно, она должна возрастать тем больше, чем больше получается отклонение экспериментальных значений Y от аппроксимирующих их YТеор, учитывая все точки распределения.
Таких критерий-функций существует несколько. Приведу пример функции Xq ("Хи-квадрат"). Вот ее вид:
Xq = S( (Yi-YТеорi)2 / Eyi2 )i=1..Nгде S означает суммирование квадратов разницы между экспериментальными и теоретическими значениями, деленными на квадраты погрешностей измерения.
Оценим поведение этой функции. Чем больше отклонение теоретического приближения от экспериментального, тем больше значение Xq, т.е. чем больше значение Xq, тем менее достоверным является приближение.
А вот как будет выглядеть запись Xq для нашего случая с линейной (y=ax) зависимостью:
X Y YТеор=ax YТеор-Y (Y-YТеор)2 1 2 2 0 0 2 4 4 0 0 3 6 6 0 0 4 7 8 1 1 5 10 10 0 0 6 13 12 -1 1 ...Итого в этом случае Xq=1+1=2.
Тонкий момент. Обратите внимание на то, что каждый квадрат отклонения делится на квадрат погрешности измерения Y - Eyi. Смысл в том, что разные точки распределения могут иметь различные ошибки, другими словами - давать различное количество информации о характере искомого распределения: чем больше погрешность точки, тем меньше ее вклад в оценку достоверности приближения. В нашем случае Eyi отсутствуют у всех точек, что эквивалентно тому, что все точки вносят одинаковый вклад в искомый вид теоретического распределения.
Оценку достоверности аппроксимации при помощи Xq можно определить по так называемым "таблицам Стьюдента", в которых приводится максимально допустимое значение Xq для аппроксимаций распределений с различным числом входных данных и различным числом параметров аппроксимирующего распределения. Однако должно быть ясно, что чем больше параметров имеет приближение и чем больше входных данных, тем большим будет минимально допустимое значение Xq, определяющее достоверность аппроксимации.
В нашем случае говорят, что зависимость Y(x) = 2,4,6,7,10,13.. аппроксимируется функцией y=ax,a=2 с точностью Xq = 2 для 6 точек, число параметров аппроксимирующей функции равно 1 (это a=2).
Применение критерия Хи-квадрат для аппроксимации.
Разберем более детально пример с аппроксимацией распределения линейной функцией y=ax и покажем, как можно найти искомый параметр распределения "a" не предполагая его значение, а вычислив его с помощью функции Xq.
Способ первый. Аналитическое решение.
Замечаем, что в случае приближения y=ax Xq - это функция от одного аргумента:
Xq = function(a);поскольку множества Y и X фиксированы на момент расчета, а YТеор вычисляется как функция от X. Тогда для того, чтобы найти минимум функции Xq достаточно решить уравнение:
dXq --- = S(2×(Yi-aXi)×(-Xi) ) = 0; daнайдя неизвестное "a".
В более общем случае теоретическое распределение может быть записано в виде:
YТеор = function( a0 , a1 , a2... ),где a0..am - параметры искомого теоретического приближения. Очевидно, уравнение нахождения неизвестных коэффициентов aj можно записать в виде:
dXq --- = 0, da1 dXq --- = 0, da2 ... dXq --- = 0, damт.е. получается система из m уравнений, образованных записью частной производной по каждому из a1..am параметру. Решая эти уравнения совместно, можно получить неизвестные параметры a1..am приближающего распределения.
Способ второй. Численные методы.
Несмотря на то, что Способ 1 может дать однозначное и наиболее точное значение параметров приближающего распределения aj, существует и другой способ определить неизвестные aj. Он заключается в поиске на некотором поле значений ajтаких, при которых значение Xq минимально. Поясню мысль рисунком:
a1 minimum -- a1 minimum+step_a1 -- a1 minimum+step_a1×2 -- ... a1 max a2 minimum -- a2 minimum+step_a2 -- a2 minimum+step_a2×2 -- ... a2 max ... am minimum -- am minimum+step_am -- am minimum+step_am×2 -- ... am maxТ.е. вначале задаются некие граничные (максимальные и минимальные) значения каждого из неизвестных параметров aj и шаг разбиения интервала (aj max - aj min) - step_aj. Каждое значение aj пробегает значения от минимума до максимума. Для каждой комбинации параметров aj (очевидно, всего возможно ((число шагов разбиения))число параметров aj вариантов) рассчитывается значение Xq и ищется такая комбинация, которая отвечает минимальному Xq
Позволю себе назвать такой способ поиском "на решетке"(lattice) значений неизвестных параметров aj. Значение step_aj назовем "шагом решетки".
После того, как на данной "решетке" значений параметров будут найдены параметры, отвечающие минимальному Xq, необходимо будет выполнить модификацию решетки с целью локализации значения каждого искомого параметра в области истинного значения параметра:
- Если значение параметра, отвечающего минимальному Xq будет обнаружено внутри решетки (т.е. он будет между минимальным и максимальным значением для данного параметра), то новые границы будут определены следующим образом:
Новое значение минимума/максимума = Значение параметра ± Шаг решетки
- Если значение параметра, отвечающего минимальному Xq будет обнаружено на краю решетки (т.е. он будет равен либо минимальному либо максимальному значению для данного параметра), то новые границы будут определены следующим образом:
Новое значение минимума/максимума = Значение параметра ± ширина решетки.
Рассмотрим пример того, как это бы происходило в случае нашей простой линейной зависимости Y(x) = 2,4,6,7,10,13.. Мы предполагаем, что у аппроксимирующей зависимости один параметр: a. Найдем его. Сначала зададим максимальное и минимальное начальное значение a:
amax = 100, amin = -100;число разбиений решетки пусть будет равно 5, тогда в первом приближении (на первой решетке) будет проверены следующие значения a:
-100 -50 0 +50 +100,
при этом шаг решетки равен:
(amax - amin) / (5-1) = ( 100 - (-100) ) / 4 = 50,
и мы вычислим пять значений функции Xq для всех значений a:
a -100 -50 0 +50 +100 Xq(a) 9.47E+0005 2.46E+0005 3.74E+0002 2.09E+0005 8.73E+0005Минимум Xq в данном случае получается если a=0. Поскольку точка лежит внутри решетки (a -100 и a 100), то новые максимальные и минимальные значения "a" необходимо выбрать следующим образом:
amax = 0 + 50, amin = 0 - 50;и перейти к следующему расчету Xq на новом интервале.
Чем больше будет выполнено таких приближений, тем точнее будут найдены искомые параметры, но до известного предела, заложенного в самом критерии Xq.
Программная реализация именно такого способа будет рассмотрена ниже. В чем его достоинства/недостатки по сравнению с Способом 1? Достоинством является универсальность относительно вида приближающей функции YТеор: нет необходимости решать сложную систему уравнений в частных производных по aj, недостатком является некоторая неточность результатов определения параметров aj, что будет показано ниже.
Реализация аппроксимации данных полиномом численными методами.
Реализуем программно решение следующей задачи: пусть имеется набор значений X,Y, предположительно образованных полиномом вида:
Y(x) = anxn + an-1x(n-1) + ... a0Попытаемся построить алгоритм, способный аппроксимировать заданную зависимость Y(X) полиномом вида Y(x) = anxn + .. +a0 : найти неизвестные коэффициенты an,an-1,..,a0 и дать оценку достоверности приближения с помощью критерия Xq.
Реализация алгоритма будет проведена на win32-ассемблере с использованием сопроцессора, все особенности работы с которым будут по возможности прокомментированы.
Заголовок программы на MASM: стандартное подключение библиотек и т.д.
.386 .model flat,stdcall option casemap:none include \masm32\include\windows.inc include \masm32\include\kernel32.inc includelib \masm32\lib\kernel32.lib include \masm32\include\user32.inc includelib \masm32\lib\user32.libОписываем данные и константы программы.
; Флаги вывода промежуточных результатов (Debug): ; выводить значение степени полинома из файла ShowPolinomPower_NOT_DEFINED equ TRUE ; выводить значения X,Y из файла ShowXYAtFileReading_NOT_DEFINED equ TRUE FloatMask equ 14 ; Маска для вывода чисел с плавающей точкой: ; 14 цифр после десятичной точки MaxPower equ 10 ; Максимальная степень входного полинома: ; y(x)=a0+a1×x+a2×x2+...aMaxPower×xMaxPower ; Размер решетки для вычислений: LatticeSize equ 4 ; (число разбиений интервала поиска) ImpLattices equ 500 ; Максимальное число итераций MaxPoints equ 40 ; Максимальное число входных точек .data ; Максимальное (по модулю) начальное значение коэффициентов MaxStartCoffValue dd 1.0E+10 BreakValueOfXq dd 0.1 ; Величина, контролирующая момент, когда ; приближение считается достоверным coff struc ; Параметры полинома: его коэффициенты Cofficient dd (MaxPower+1) dup (?) coff ends Point struc ; Описание одной точки: пара (X,Y) X dd? Y dd? Point ends .data? Power dd? ; Предполагаемая степень из файла EntryData db MaxPoints×sizeof(Point) dup(?) ; Входные точки (X,Y) NumberDataPoints dd? ; Число точек в файле CurrLattice dd? ; Текущий номер решетки TempReal dd? ; Временная переменная ; Массивы коэффициентов полинома: amin coff <> ; Набор минимальных коэффициентов amax coff <> ; Набор максимальных коэффициентов awork coff <> ; Набор промежуточных коэффициентов afoundmin coff <> ; Набор коэффициентов, отвечающих минимальному Xq Xq dd? ; Значение Xq Xqmin dd? ; значение минимального Xq на очередном шаге вычислений ; Набор индексов для прохождения узлов решетки LatticeIndexes dd (MaxPower+1) dup(?) ; Индексы, отвечающие меньшему Xq NewIndexes dd (MaxPower+1) dup(?)Константы ShowPolinomPower и ShowXYAtFileReading определяются для вывода отладочных сообщений при чтении параметров из файла. Они начнут работать, если удалить "_NOT_DEFINED".
FloatMask управляет выводом чисел с плавающей точкой: в данном случае эта величина определяет вывод 14 знаков после десятичной точки.
MaxPower определяет максимально допустимую степень аппроксимирующего полинома, поддерживаемого программой - 10. MaxPoints определяет максимально возможный размер входных данных (точек x,y)
LatticeSize является количеством разбиений интервала [максимальное значение - минимальное значение] каждого коэффициента полинома. ImpLattices определяет число итераций в численном поиске минимального Xq (число анализируемых решеток). MaxStartCoffValue определяет начальные значения всех коэффициентов полинома: вначале полагаем макисмально/минимально возможными значения коэффициентов равными ± MaxStartCoffValue. BreakValueOfXq определяет точность определения коэффициентов полинома: чем больше этот порог, тем хуже будет определены коэффициенты полинома: на самом деле эта величина должна быть согласована с таблицей Стьюдента.
.data ErrTitle db "Error? ;)",0 ; Заголовок сообщения отладки/ошибки .code start: .data AppName db "Splain test program",0 Welcome db "Example of Approximate Data by Polinom. Coded by Chingachguk, 2002.",10,10 db "Entry Data(must be in file 'splain.dat'):",10 db "<polinom power>",10 db "<X1>,<Y1>",10 db "...",10 db "<Xn>,<Yn>",0 .code invoke MessageBox,NULL,addr Welcome,addr AppName,MB_OK+MB_ICONINFORMATION .data IniFileName db "splain.dat",0 ; Имя файла данных .data? IniDescr dd? ; Дескриптор файла .code ; ×××××××××××××××× Read Entry File ××××××××××××××××××××××××××××× ; Open entry file (not Create !) invoke CreateFile,offset IniFileName,GENERIC_WRITE or GENERIC_READ,\ FILE_SHARE_WRITE or FILE_SHARE_READ,0,OPEN_EXISTING,FILE_ATTRIBUTE_NORMAL,0 mov IniDescr,eax ; Save file descriptor cmp eax,0ffffffffh jz FileNotOpen .data? CurrString dd? FileStringSize equ 80 FileString db FileStringSize dup(?) .code ; Ведем счет строкам входного файла mov dword ptr CurrString,1 ; Читаем предполагаемую степень полинома (max: MaxPower) push FileStringSize ; размер приемника push offset FileString ; Адрес приемника строки push IniDescr ; Дескриптор файла call ReadFileString ; Переводим строку со степенью в число (real) push offset Power push offset FileString call StringToVal jc InvalidNumericFormat ; Переводим степень из real в int fld dword ptr Power ; Задаем режим округления (биты 11 - отбросить дробную часть) call Set_FPUTruncType frndint fistp dword ptr Power ; Проверяем степень на валидность cmp dword ptr Power,0 jl InvalidPowerValue cmp dword ptr Power,MaxPower ja InvalidPowerValue ; ××× Отладочное сообщение: значение предполагаемой степени полинома ××× ifdef ShowPolinomPower .data PolinomPowerName db "Power: " .code mov edi,offset FileString mov esi,offset PolinomPowerName mov ecx,sizeof PolinomPowerName cld rep movsb mov eax,Power mov ebx,10000000 call DecChar mov byte ptr [edi+8],0 invoke MessageBox,0,offset FileString,offset ErrTitle,MB_OK endif ; Читаем пары точек (X,Y) из файла mov dword ptr NumberDataPoints,0 ; Стартовое значение базы массива EntryData (ebp) xor ebp,ebp @@ReadPoints: ; Ведем счет строкам входного файла inc dword ptr CurrString ; Читаем X push FileStringSize ; размер приемника push offset FileString ; Адрес приемника строки push IniDescr ; Дескриптор файла call ReadFileString test eax,eax jz @@ReadPointsDone ; Переводим строку в число (real) lea eax,EntryData.Point[ebp].X push eax push offset FileString ; Адрес строки call StringToVal jc InvalidNumericFormat ; Читаем Y push FileStringSize ; размер приемника push offset FileString push IniDescr call ReadFileString test eax,eax jz @@ReadPointsDone ; Переводим строку в число (real) lea eax,EntryData.Point[ebp].Y push eax push offset FileString call StringToVal jc InvalidNumericFormat ; ××× Отладочное сообщение: значения X,Y ××× ifdef ShowXYAtFileReading cld mov edi,offset FileString mov eax,' :X ' stosd push FloatMask ; Маска вывода lea eax,EntryData.Point[ebp].X push eax ; Число push edi ; Изображение числа call ValToString push FloatMask lea eax,EntryData.Point[ebp].Y push eax ; Число mov esi,offset FileString call Str_Len lea edi,[esi+ecx] mov eax,' :Y ' stosd push edi ; Изображение числа call ValToString invoke MessageBox,0,offset FileString,offset ErrTitle,MB_OK endif ; Максимальное число входных точек не должно быть больше MaxPoints inc dword ptr NumberDataPoints cmp dword ptr NumberDataPoints,MaxPoints ja InvalidPointsValue ; Индексируем следующие точки в массиве add ebp,sizeof Point jmp @@ReadPoints @@ReadPointsDone: ; ×××××××××× approximate Data ××××××××××××× ; Прочитав все точки из файла, задаем начальные значения коэффициентов ; set initial coefficient's values (равные ±MaxStartCoffValue) mov ecx,Power inc ecx xor esi,esi @@SetStartCoffs: fld MaxStartCoffValue fchs ; Сделать отрицательным fstp amin.Cofficient[esi] fld MaxStartCoffValue fstp amax.Cofficient[esi] ; Следующий коэффициент add esi,sizeof coff / (MaxPower+1) loop @@SetStartCoffs ; ××××××× Основной цикл: поиск коэффициентов, отвечающих минимальному Xq ××××××× ; В этом цикле мы вычисляем Xq на всех узлах решетки, ; находим минимальный Xq и запоминаем отвечающие ему ; значения коэффициентов решетки, ; затем сжимаем/сдвигаем решетку минимального Xq в зависимости от того, ; лежат ли коэффициенты на краю (внутри) решетки. ; Если значение Xq показывает, что приближение достоверно, ; мы заканчиваем основной цикл поиска ; Номер решетки приближения mov dword ptr CurrLattice,1 NextLattice: ; find Xq minimum at the current Lattice ; Получить стартовое значение Xq ; (Получить значение полинома с минимальными коэффициентами (массив amin)) fldz fstp Xqmin xor esi,esi ; Индекс - esi - номер точки из файла mov ecx,NumberDataPoints ; Число точек в файле @@GetStartXq: push dword ptr EntryData.Point[esi].X ; значение аргумента (X) push Power ; размер (степень) полинома push offset amin ; адрес массива коэффициентов полинома call polinom ; Результат: значение полинома в ST(0) ; вычислить разницу -(Yi-polinom(Xi)) fsub dword ptr EntryData.Point[esi].Y ; Возвести в квадрат разницу -(Yi-polinom(Xi)), она в ST(0) mov eax,2 ; степень - 2 call get_pow ; вернется в ST(0) квадрат fadd dword ptr Xqmin ; Прибавить к накопителю Xqmin fstp dword ptr Xqmin ; Запомнить накопитель add esi,sizeof Point loop @@GetStartXq ; Обнуляем стартовые коэффициенты для прохода по значениям коэффициентов mov ecx,Power inc ecx mov edi,offset LatticeIndexes cld xor eax,eax rep stosd ; выполняем поиск Xq на текущих значениях коэффициентов FindXqMinimum: ; Определяем коэффициенты текущего приближенного полинома (awork) xor esi,esi ; Индексы - esi, edi xor edi,edi mov ecx,Power ; По всем степеням inc ecx @@GetCurrentCoffs: ; вычисляем разницу между максимальным и минимальным значением коэффициента fld amax.Cofficient[esi] fsub amin.Cofficient[esi] ; Делим разницу на число (узлов решетки-1), получая шаг решетки mov TempReal,LatticeSize-1 fidiv dword ptr TempReal ; Умножаем на номер узла текущего коэффициента fimul dword ptr LatticeIndexes[edi] ; Получаем значение коэффициента в узле, добавляем минимальное значение fadd dword ptr amin.Cofficient[esi] ; вычисляем текущие коэффициенты полинома (awork) fstp dword ptr awork.Cofficient[esi] ; Следующий коэффициент add esi,sizeof coff / (MaxPower+1) add edi,4 ; Следующий номер узла loop @@GetCurrentCoffs ; вычисляем Xq для промежуточных значений коэффициентов (awork) fldz fstp Xq xor esi,esi ; Индекс - esi - номер точки из файла mov ecx,NumberDataPoints ; Число точек в файле @@GetCurrentXq: push dword ptr EntryData.Point[esi].X push Power ; размер (степень) полинома push offset awork ; адрес массива коэффициентов полинома call polinom ; Результат: значение полинома в ST(0) ; вычислить разницу -(Yi-polinom(Xi)) fsub dword ptr EntryData.Point[esi].Y ; Возвести в квадрат разницу -(Yi-polinom(Xi)), она в ST(0) mov eax,2 ; степень - 2 call get_pow ; вернется в ST(0) квадрат fadd dword ptr Xq ; Прибавить к накопителю Xqmin fstp dword ptr Xq add esi,sizeof Point loop @@GetCurrentXq ; Сравниваем Xq и Xqmin fld Xq fcomp dword ptr Xqmin ; Сравнить Xq с Xqmin и вытолкнуть из стека ST(0) fstsw ax ; Перенести флаги в регистр ax sahf ; Дублировать их в регистре флагов jnc @@CurrXqNoLow ; Если Xq < Xqmin, то присвоить новое значение Xqmin ; Т.е. мы нашли коэффициенты с меньшим Xq, запомним их fld dword ptr Xq fstp dword ptr Xqmin ; Также запомнить текущие коэффициенты в afoundmin mov ecx,Power inc ecx mov esi,offset awork mov edi,offset afoundmin cld push ecx rep movsd pop ecx ; Запомнить текущие номера узлов решетки mov esi,offset LatticeIndexes mov edi,offset NewIndexes rep movsd @@CurrXqNoLow: ; Инкрементировать номера узлов решетки, ; перебирая все узлы решетки xor esi,esi @@NextApprox: inc dword ptr LatticeIndexes[esi*4] cmp dword ptr LatticeIndexes[esi*4],LatticeSize jb FindXqMinimum ; Следующая решетка mov dword ptr LatticeIndexes[esi*4],0 inc esi cmp esi,dword ptr Power jbe @@NextApprox ; Мы прошли по всем узлам решетки ; Анализируем, где был найден минимум Xq mov ecx,Power inc ecx xor esi,esi @@AnaliseLattice: mov eax,dword ptr NewIndexes[esi] test eax,eax jz @@XqMinLieOnRange cmp eax,LatticeSize-1 jnz @@XqMinNoLieOnRange @@XqMinLieOnRange: ; point lie on the rage of lattice: no change size of lattice ! ; Строим решетку симметрично относительно найденного узла ; Новый минимум/максимум=найденный узел±ширина решетки ; Получаем в ST(0) новое значение минимального коэффициента fld dword ptr afoundmin.Cofficient[esi] fsub dword ptr amax.Cofficient[esi] fadd dword ptr amin.Cofficient[esi] ; Получаем в ST(0) новое значение максимального коэффициента, мин.->ST(1) fld dword ptr afoundmin.Cofficient[esi] fadd dword ptr amax.Cofficient[esi] fsub dword ptr amin.Cofficient[esi] jmp @@NextLatticePoint @@XqMinNoLieOnRange: ; point lie inside of lattice: minimize size of lattice ; Строим решетку симметрично относительно найденного узла ; Новый минимум/максимум=найденный узел±шаг решетки ; Получаем в ST(0) новое значение минимального коэффициента fld dword ptr amin.Cofficient[esi] fsub dword ptr amax.Cofficient[esi] ; Делим разницу на число (узлов решетки-1), получая шаг решетки mov TempReal,LatticeSize-1 fidiv dword ptr TempReal fadd dword ptr afoundmin.Cofficient[esi] ; Получаем в ST(1) новое значение максимального коэффициента, ST(0)->ST(1) fld dword ptr amax.Cofficient[esi] fsub dword ptr amin.Cofficient[esi] ; Делим разницу на число (узлов решетки-1), получая шаг решетки fidiv dword ptr TempReal fadd dword ptr afoundmin.Cofficient[esi] @@NextLatticePoint: ; ST(0) - новый максимум, ST(1) - новый минимум fstp dword ptr amax[esi] ; <-ST(0) fstp dword ptr amin[esi] ; <-ST(0) add esi,4 loop @@AnaliseLattice ; Переходим к следующей решетке, выполняя не более ImpLattices приближений inc dword ptr CurrLattice ; Проверяем величину текущего минимума Xq ; Возможно, мы уже нашли коэффициенты fld dword ptr BreakValueOfXq fmul dword ptr Power fadd dword ptr BreakValueOfXq fimul dword ptr NumberDataPoints ; ST(0)=NumberDataPoints × BreakValueOfXq × (Power+1) ; Если величина Xqmin =< NumberDataPoints×BreakValueOfXq×(Power+1), ; то приближение считается достоверным ; Сравнить NumberDataPoints×BreakValueOfXq×(Power+1) с Xqmin и вытолкнуть из стека ST(0) fcomp dword ptr Xqmin fstsw ax ; Перенести флаги в регистр ax sahf ; Дублировать их в регистре флагов jnc @@ApproximateValid cmp dword ptr CurrLattice,ImpLattices jbe NextLattice @@ApproximateValid: ; ×××××××××× Цикл поиска коэффициентов закончен ××××××××××××××××× ; выводим полученные коэффициенты, значение Xq для анализа достоверности ; нашего приближения данных из файла полиномом .data ResultTitle db "Calculations resume",0 ResultHeader0 db "Assume that you data was approximated by polinom with power????.",10,10 ResultHeaderD db "a[?]=?.???????????????????? " ResultHeaderE db "Xq =?.???????????????????? for????? points.",10 .data? ResultHeader db 1000 dup(?) .code cld mov edi,offset ResultHeader ; выводим степень полинома push edi mov eax,Power mov ebx,1000 mov edi,offset ResultHeader0+60 call DecChar pop edi mov esi,offset ResultHeader0 mov ecx,sizeof ResultHeader0 rep movsb ; выводим полученные коэффициенты mov ecx,Power ; По всем степеням inc ecx xor esi,esi @@OutCoffs: pushad mov eax,Power inc eax sub eax,ecx mov ebx,1 mov edi,offset ResultHeaderD+2 call DecChar push FloatMask ; Маска вывода lea eax,afoundmin.Cofficient[esi] push eax push offset ResultHeaderD+6 ; Изображение числа call ValToString popad push esi push ecx mov esi,offset ResultHeaderD mov ecx,sizeof ResultHeaderD @@CopyCoff: lodsb test al,al jnz @@UsalSym mov al,10 mov ecx,1 @@UsalSym: stosb loop @@CopyCoff pop ecx pop esi ; Следующий коэффициент add esi,sizeof coff / (MaxPower+1) loop @@OutCoffs mov al,10 stosb ; выводим Xqmin push edi push FloatMask ; Маска вывода push offset Xqmin ; Число push offset ResultHeaderE+5 ; Изображение числа call ValToString mov byte ptr ResultHeaderE+27,' ' mov eax,NumberDataPoints mov ebx,10000 mov edi,offset ResultHeaderE+32 call DecChar pop edi mov esi,offset ResultHeaderE mov ecx,sizeof ResultHeaderE rep movsb xor al,al stosb invoke MessageBox,0,offset ResultHeader,offset ResultTitle,MB_OK+MB_ICONINFORMATION Exit: invoke CloseHandle,IniDescr invoke ExitProcess,NULL ;×××××××××××× Подпрограмма вычисления значения полинома ×××××××××××××××××××××××× ; Get value of LocCoff[i]×xi, i=0..power ; [ebp+8] - адрес массива коэффициентов полинома ; [ebp+0Ch] - размер (степень) полинома ; [ebp+10h] - значение аргумента (X) ; Результат: значение полинома в ST(0) polinom proc push ebp mov ebp,esp ; result: [ebp-4], локальная степень полинома: [ebp-8] sub esp,2*4 push esi ; Определяем начальный результат fldz fstp dword ptr [ebp-4] ; Загружаем адрес коэффициентов полинома mov esi,[ebp+8] ; Степень полинома fldz fstp dword ptr [ebp-8] @@polinom: ; Получить XТекущая степень (0..[ebp+0Ch]) fld dword ptr [ebp+10h] ; Загрузить X mov eax,dword ptr [ebp-8] ; степень call get_pow ; вернется в ST(0) ; Умножить на текущий коэффициент fmul dword ptr [esi] ; Добавить сумматор fadd dword ptr [ebp-4] ; Запомнит новый результат fstp dword ptr [ebp-4] add esi,sizeof coff / (MaxPower+1) ; Следующий коэффициент inc dword ptr [ebp-8] mov eax,dword ptr [ebp-8] cmp eax,[ebp+0Ch] jbe @@polinom ; Возвращаем в ST(0) результат fld dword ptr [ebp-4] pop esi leave ret 3*4 polinom endp ;×××××××××××××××××××××××××××××××× Сообщения об ошибках ××××××××××××××××××××××××× FileNotOpen: .data IniOpOK db "Entry Data File (splain.dat) not open (not exist?) !",0 .code invoke MessageBox,0,offset IniOpOK,offset ErrTitle,MB_OK or MB_ICONERROR jmp Exit InvalidNumericFormat: .data InvalidNumericFormatMsg db "Invalid numeric format in string " ErrStrNum db "00000000",0 .code mov eax,CurrString mov ebx,10000000 mov edi,offset ErrStrNum call DecChar invoke MessageBox,0,offset InvalidNumericFormatMsg,offset ErrTitle,MB_OK or MB_ICONERROR jmp Exit InvalidPowerValue: .data InvalidPowerValueMsg db "Invalid Power Value",0 .code invoke MessageBox,0,offset InvalidPowerValueMsg,offset ErrTitle,MB_OK or MB_ICONERROR jmp Exit InvalidPointsValue: .data InvalidPointsValueMsg db "Too much data points in file !",0 .code invoke MessageBox,0,offset InvalidPointsValueMsg,offset ErrTitle,MB_OK or MB_ICONERROR jmp ExitКратко о вспомогательных программах. Для реализации алгоритма нам потребуется несколько подпрограмм:
- семейство Set_FPU..Type. Подпрограммы для задания режима округления сопроцессора при операциях конвертирования чисел с плавающей точкой в числа типа inteher.
- get_pow. Подпрограмма возведения в целую степень числа с плавающей точкой. Реализована примитивно умножением аргумента самого на себя степень раз.
- StringToVal: Подпрограмма деформатирования строки в число типа Real. Позволяет читать из файла числа с плавающей запятой в виде строки и преобразовывать их в формат real (4 байта). Реализована следующим образом: накопитель инициализируется числом 0, затем читается цифра из строки, добавляется к приемнику, приемник предварительно доумножается на 10. Десятичная точка учитывается подсчетом числа цифр после нее и финальным делением на 10 в степени "число цифр после десятичной точки".
- ValToString: Подпрограмма форматирования числа типа Real в строку. Позволяет выполнять форматный вывод чисел с плавающей точкой с заданным числом знаков после десятичной точки. Алгоритм вывода довольно прост и, возможно, допускает некоторую погрешность: определяется максимальная степень десяти в числе при помощи команды fyl2x сопроцессора: фактически вычисляется max_power = int( log10(x) ). Далее последовательным делением на 10 в степени max_power, max_power-1, ... определяются коэффициенты разложения числа по степеням десяти.
- ReadFileString: парсер входного файла. Позволяет разделять строки с числовыми данными из входного файла, отделенные стандартными разделителями (10 или 13+10 или пробелами) и читать их в заданный буфер.
; ×××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××× ; Вспомогательные подпрограммы ; ×××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××× .data DecNumber dd 10.0 ; ×××××××××××× Set_FPU...Type: Задать режим округления сопроцессора ×××××××××××× ; Set_FPURoundType: округление до ближайшего целого (00) ; Set_FPUTruncType: отбрасывание дробной части (11) .data? control dw? .code Set_FPURoundType proc fstcw word ptr control and word ptr control,0F3FFh fldcw word ptr control ret Set_FPURoundType endp Set_FPUTruncType proc fstcw word ptr control or word ptr control,0C00h fldcw word ptr control ret Set_FPUTruncType endp ; ×××××××××××× get_pow: Подпрограмма возведения в целую степень числа ×××××××××× ; ST(0) - аргумент, eax - степень ; result: ST(0) ; Для возведения в произвольную стпень числа следует воспользоваться ; более сложным алгоритмом с использованием fyl2x: ; см., например, П.И. Рудаков, К.Г. Финогенов ; "Программируем на языке ассемблера IBM PC" .code get_pow proc test eax,eax jnz @@NoZeroPower fstp ST(0) ; вытолкнуть аргумент из стека fld1 ; Загрузить 1.0 ret ; Вернуть 1.0 @@NoZeroPower: push ecx mov ecx,eax ; Сделать степень положительной @@NegLoop: neg ecx js @@NegLoop fld1 ; Загрузить 1.0, ST(0) - накопитель @@get_pow: fmul ST(0),ST(1) ; Накопитель=Накопитель x arg loop @@get_pow fxch ST(1) ; Поменять накопитель и аргумент fstp ST(0) ; вытолкнуть аргумент из стека ; Если степень отрицательна, необходимо Накопитель=1/Накопитель cmp eax,0 jge @@NoNegPow fld1 ; Загрузить 1 fxch ST(1) ; Обменять ST(1) и ST(0) fdivp ST(1),ST(0) ; Разделить 1 на Накопитель @@NoNegPow: pop ecx ret get_pow endp ; ×××× StringToVal: Подпрограмма деформатирования строки в число типа Real ××××× ; [ebp+8] - Адрес строки(Zero-terminated) ; [ebp+0Сh] - Адрес приемника (real или dd) ; Десятичным разделителем считается символ "." ; Пробелы перед цифрами пропускаются ; Ошибка - флаг CF установлен .code StringToVal proc push ebp mov ebp,esp sub esp,1*4 ; Число из строки [ebp-4] mov ebx,[ebp+0Ch] ; Адрес приемника храним в ebx fldz ; Обнуляем приемник fstp dword ptr [ebx] cld ; Запоминаем знак числа xor edi,edi mov esi,dword ptr [ebp+8] lodsb dec esi cmp al,'-' jnz @@NoSignBefore ; Устанавливаем 31-ый бит знака числа inc edi ror edi,1 inc esi @@NoSignBefore: xor edx,edx ; Флаг десятичной точки и число знаков после нее @@ParseString: xor eax,eax lodsb test al,al jz @@EndParse cmp al,' ' ; Пропускаем пробелы слева jz @@ParseString cmp al,'.' jnz @@NoDecPoint test dh,dh jnz @@ErrorSign mov dx,256 jmp @@ParseString @@NoDecPoint: sub al,'0' js @@ErrorSign cmp al,9 ja @@ErrorSign mov dword ptr [ebp-4],eax fld dword ptr [ebx] ; Загружаем текущее значение результата fmul DecNumber ; Умножаем результат на 10 fiadd dword ptr [ebp-4] ; Добавляем к результату текущее число fstp dword ptr [ebx] ; Запоминаем получившееся значение ; Учитываем число знаков после десятичной точки inc dl jmp @@ParseString @@EndParse: ; Учесть десятичную точку test dh,dh jz @@ThereAintPoint mov dh,0 ; Разделить результат на 10edx fld dword ptr DecNumber ; Загрузить аргумент (10) mov eax,edx ; Целая степень call get_pow ; 10edx вернется в ST(0) fld dword ptr [ebx] ; Загрузить результат без точки fxch ST(1) ; Обменять ST(1) и ST(0) fdivp ST(1),ST(0) fstp dword ptr [ebx] ; Res'=Res/10edx @@ThereAintPoint: ; Учесть знак числа or dword ptr [ebx],edi clc leave ret 2*4 @@ErrorSign: stc leave ret 2*4 StringToVal endp ; ××××× ValToString: Подпрограмма форматирования числа типа Real в строку ×××××× ; [ebp+8] - Адрес строки(будет Zero-terminated) в виде 12345.89 ; [ebp+0Сh] - Адрес числа типа real (или dd) ; [ebp+10h] - Маска вывода числа (byte) ValToString proc .code push ebp mov ebp,esp ; Локальные переменные: ; [ebp-4] - модуль числа ; [ebp-08] - максимальная степень 10-ти в числе ; [ebp-0Ch] - временная sub esp,4*3 cld mov edi,dword ptr [ebp+8] ; Адрес выходной строки mov ebx,[ebp+0Ch] ; Загружаем адрес числа в ebx mov edx,[ebx] ; Читаем его в edx ; Определяем знак, удаляя его одновременно (31 бит) shl edx,1 jnc @@NoSign mov al,'-' stosb @@NoSign: shr edx,1 mov [ebp-4],edx ; Запоминаем модуль числа в локальной переменной ; Проверка на ноль test edx,edx jnz @@NonZero mov al,'0' stosb jmp @@DoneValToStr @@NonZero: ; Определяем старшую степень десяти - Log10(x) = Log2(x)/Log2(10) ; выполняется в два приема, поскольку нет операции Log10 в FPU, ; а есть y×Log2(x) ; Получаем 1/Log2(10) fld1 ; Загрузить 1.0 fldl2t ; Загрузить log2(10) fdivp ST(1),ST(0) ; Разделить 1/log2(10), результат в ST(0) fld dword ptr [ebp-4] ; Загрузить наше число, 1/log2(10) протолкнуть в ST(1) fyl2x ; ST(0)=ST(1)×log2(ST(0)) ; Задаем режим округления (биты 11 - отбросить дробную часть) call Set_FPUTruncType frndint fistp dword ptr [ebp-8] ; выгрузить степень как целое со знаком ; Задаем режим округления (биты 00 - округлить до ближайшего целого) call Set_FPURoundType ; Сохранить степень числа push dword ptr [ebp-8] ; вывести число, деля его на 10... ; вывести целую часть числа, точку, дробную часть mov ecx,dword ptr [ebp+10h] ; Сколько цифр после десятичной точки выводить inc ecx ; Цифра перед десятичной точкой xor edx,edx ; Флаг вывода точки @@OutNumbers: fld dword ptr DecNumber ; Загрузить 10 - основание mov eax,dword ptr [ebp-8] ; Степень call get_pow ; ST(0) = 10Степень fld dword ptr [ebp-4] ; Загрузить число ; ST(0)=Число, ST(1) = 10Степень fdiv ST(0),ST(1) ; ST(0)=ST(0)/ST(1) ; Задаем режим округления (биты 11 - отбросить дробную часть) call Set_FPUTruncType frndint fistp dword ptr [ebp-0Ch] ; выгружаем как целое без округления ; Задаем режим округления (биты 00 - округлить до ближайшего целого) call Set_FPURoundType ; выводим цифру mov eax,dword ptr [ebp-0Ch] mov ebx,1 call DecChar ; В формате - @n1 inc edi ; выводим точку, если это необходимо test edx,edx jnz @@PointAlreadyOut mov al,'.' stosb inc edx @@PointAlreadyOut: ; вычесть из числа 10Степень(сейчас в ST(0)) x Текущая цифра fimul dword ptr [ebp-0Ch] fld dword ptr [ebp-4] ; Загрузить число fsub ST(0),ST(1) fstp dword ptr [ebp-4] ; выгрузить число с вычетом fstp ST(0) ; Очистить стек ; Уменьшить степень на 1 dec dword ptr [ebp-8] loop @@OutNumbers ; Восстановить степень числа pop dword ptr [ebp-8] ; вывести степень 10-ти mov ax,'+E' mov edx,[ebp-8] cmp edx,0 jge @@NoNegPower neg edx mov ah,'-' @@NoNegPower: stosw mov eax,edx mov ebx,1000 call DecChar ; В формате @n4 add edi,4 @@DoneValToStr: mov byte ptr [edi],0 leave ret 3*4 ValToString endp ; ××××××××××××× ReadFileString: Подпрограмма чтения строки из файла ×××××××××××× ; [ebp+8] - Дескриптор файла ; [ebp+0Ch] - Адрес приемника строки ; [ebp+10h] - размер приемника ; -> eax реальный размер строки в приемнике и ноль в конце приемника ; Разделителями считаются символы 0A или 0D,0A, пробел. ; Пробелы слева читаются до любого другого символа, ; затем считываются символы числа, пробел служит знаком конца строки ReadFileString proc .code push ebp mov ebp,esp ; Переменная для байта из файла размером в 1 байт [ebp-8] ; Число байт, прочтенных ReadFile [ebp-4] ; Число байт в строке [ebp-0Ch] sub esp,3*4 cld mov edi,[ebp+0Ch] mov dword ptr [ebp-0Ch],0 ; Флаг чтения любого символа, кроме пробела xor esi,esi @@NextByte: push NULL lea eax,[ebp-4] push eax push 1 lea eax,[ebp-8] push eax push dword ptr [ebp+8] call ReadFile cmp dword ptr [ebp-4],0 jz @@EndRead mov al,byte ptr [ebp-8] cmp al,0Ah jz @@EndRead cmp al,0Dh jz @@NextByte ; Пропускаем байт 0Ah cmp al,' ' ; Пробел? jnz @@NoFiller test esi,esi jnz @@EndRead ; Пробел после строки - break jmp @@AfterFiller ; Пробелы перед строкой заносим в буффер @@NoFiller: inc esi ; Установим флаг, что был символ <> пробел @@AfterFiller: mov ebx,dword ptr [ebp-0Ch] inc ebx cmp ebx,dword ptr [ebp+10h] jae @@EndRead inc dword ptr [ebp-0Ch] stosb jmp @@NextByte @@EndRead: mov byte ptr [edi],0 mov eax,dword ptr [ebp-0Ch] leave ret 3*4 ReadFileString endp ; ××××××××××××× DecChar: Подпрограмма форматирования строки из числа ××××××××××× ; eax=number to digit, edi=offset result string in format 00000000(@n[ebx]) ; ebx=начальный делитель DecChar proc pushad pushfd cld @GetDec: xor edx,edx div ebx add al,'0' stosb push edx mov eax,ebx xor edx,edx mov ebx,10 div ebx mov ebx,eax pop eax test ebx,ebx jnz @GetDec popfd popad ret DecChar endp ; ××××××××××××××××××××××× Str_Len: Get lenght of string ××××××××××××××××××××××××× ; esi=addr of string; result: ecx=length Str_Len proc xor ecx,ecx push eax push esi @@GetStrLen: lodsb test al,al jz @@ExitStrLen inc ecx jmp @@GetStrLen @@ExitStrLen: pop esi pop eax ret Str_Len endp end startОценим работоспособность алгоритма на его способности аппроксимировать Y(X) полиномом вида Y(x) = a2x2 + a1x+a0 (параболой): найти неизвестные коэффициенты a2, a1, a0.
Для примера рассмотрим простую параболу:
Y(x) = 3x2 - 5x -8;на интервале значений x={1..10}:
X Y 1.0 -10.00 2.0 -6.00 3.0 4.00 4.0 20.00 5.0 42.00 6.0 70.00 7.0 104.00 8.0 144.00 9.0 190.00 10.0 242.00Поместим эту зависимость Y(x) в файл входных данных программы (splain.dat), указав предполагаемую степень полинома первой строкой:
Содержимое входного файла splain.dat: 2 1.0 -10.00 2.0 -6.00 ... 10.0 242.00В результате программа должна выдать следующий результат (содержимое MessageBox'а):
Assume that you data was approximated by polinom with power 0002 a[0]= -7.282... a[1]= -5.288... a[2]= 3.0247... Xq = 0.413... for 00010 points.Как видно, программа достаточно точно смогла определить коэффициенты параболы. Заметим, однако, что несмотря на то, что входные данные являлись абсолютно точно параболой с опредленными коэффициентами, они не были определены на 100%. Почему так происходит? Дело в том, что критерий Xq выполняет подбор коэффициентов по величине квадратов отклонения точек экспериментальных данных от теоретических. Очевидно, что для небольших значений Y приближение будет менее точным, нежели чем для больших: см. например, то, что значение a[2] определено гораздо точнее, чем a[1] или тем более a[0].
Оценка предложенного метода численной аппроксимации с помощью Xq.
В результате проведенного теста оказывается, что предложенный способ достаточно точно позволяет определять параметры аппроксимирующей зависимости, по крайней мере позволяет выполнять предварительные оценки свойств экспериментальных зависимостей. Вместе с этим он претендует на некоторую универсальность алгоритма относительно выбора вида и параметров аппроксимирующей зависимости: она может быть практически сколь угодно сложной и иметь широкий спектр параметров; алгоритм не требует сложных математических вычислений. © Chingachguk / HI-TECH
Сплайн. Аппроксимация данных
Дата публикации 21 ноя 2002
