В данном документе речь пойдет об интерпретации экспериментальных данных.
Под экспериментальными данными понимается некий набор измерений, являющийся результатом опыта.
Так, например вы стоите на выходе из метро и считаете число людей, проходящих мимо вас за каждую минуту. Допустим, в результате у вас получились следующие измерения: 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, троек бит... Все они должны быть равновероятны, т.е. вероятность выпадения в каждом ключе определенного числа нулевых бит должно описывается так называемым биномиальным распределением:
где n - полное число бит в ключе (в нашем случае 8?8 = 64), k - определенное число бит, Ckn - число вариантов размещения k нулевых бит по n свободным местам:
k
C = n! / k! ? (n-k)!;
n
2n - число всех битовых комбинаций в числе из 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 - это функция от одного аргумента:
поскольку множества 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:
число разбиений решетки пусть будет равно 5, тогда в первом приближении (на первой решетке) будет проверены следующие значения a:
при этом шаг решетки равен:
(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 masm32includewindows.inc
include masm32includekernel32.inc
includelib masm32libkernel32.lib
include masm32includeuser32.inc
includelib masm32libuser32.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. Подпрограммы для задания режима округления сопроцессора при операциях конвертирования чисел с плавающей точкой в числа типа integer.
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
wasm.ru