Исследование методов прогнозирования временных рядов с использованием дифференциальных уравнений
 Скачать 2.03 Mb.
|
|
Реферат Работа включает: 34 страницы, 19 рисунков, 16 использованных источников. Объект исследования – временные ряды. Предмет исследования – прогнозирование временных рядов с использованием дифференциальных уравнений. Ключевые слова: временные ряды, дифференциальные уравнения, методы прогнозирования. Цель и задачи исследования – Исследование методов прогнозирования временных рядов с использованием дифференциальных уравнений. В работе исследованы следующие задачи: Изучена литература по теории прогнозирования временных рядов. Рассмотрены классические методы прогнозирования временных рядов. Исследованы методы прогнозирования временных рядов с использованием дифференциальных уравнений. Составлен автоматизированный алгоритм на языке программирования python и в вычислительной среде MATLAB для расчёта прогнозирования временного ряда и построения соответствующих графиков. Содержание Введение ............................................................................................. 2 1. Классические методы прогнозирования .................................................. 3 1.1 Метод экспоненциального сглаживания. .................................... 3 1.2 Метод Хольта. ............................................................................... 6 2. Постановка задачи о прогнозировании временного ряда дифференциальным уравнением ................................................................... 10 3. Метод решения ........................................................................................... 13 4. Результаты .................................................................................................. 16 4.1 Первый метод .............................................................................. 16 4.2 Второй метод .............................................................................. 20 Заключение ....................................................................................................... 25 Список литературы.......................................................................................... 26 ПРИЛОЖЕНИЕ А Листинг кода алгоритма 1.Python ............................... 28 ПРИЛОЖЕНИЕ Б Листинг кода алгоритма 2.Matlab ............................... 33 2 Введение Ряды - дискретные последовательности значений величины, характеризующей объект, измеренные или вычисленные при некоторых значениях независимой переменной. Если независимой переменной является время, ряды называют временными Термин "прогноз" (от греч, prognosis) традиционно обозначает предвидение, предсказание о развитии чего-либо, основанное на определенных данных. Прогнозирование - это оценка, предвидение, предсказание будущего развития избранного объекта управления. Оно должно предшествовать планированию. Главная цель прогнозирования - это формирование научных предпосылок принятия управленческих решений. Данные предпосылки включают: Анализ тенденций и закономерностей изменения объекта управления; Альтернативное предвидение его будущего развития; Оценку последствий активного воздействия на предвидимые процессы в объекте управления. Главные задачи прогнозирования: 1. Научный анализ социально-экономических процессов и явлений, оценка сложившейся ситуации и выявление узловых проблем развития. 2. Оценка действия этих тенденций в будущем, предвидение новых экономических ситуаций, новых проблем, требующих своего разрешения. 3. Выявление возможных альтернатив развития в перспективе. 3 1. Классические методы прогнозирования Существует много моделей прогнозирования: регрессионные модели прогнозирования, авторегрессионные модели прогнозирования, модели экспоненциального сглаживания, модель по выборке максимального подобияи т.д. Мы же в данной работе проведём обзор 2 из них: Метод экспоненциального сглаживания и Метод Хольта. Оба метода оптимальны для прогнозирования на небольшой промежуток времени с малым количеством начальных данных. 1.1 Метод экспоненциального сглаживания. Прогноз по методу экспоненциального сглаживания – оптимальный вариант прогноза, когда данные есть только за несколько периодов (месяцев, дней, недель, кварталов) и еще не понятно - существуют ли тенденции роста или падения Ŷ t+1 =k*Y t +(1-k) * Ŷ t, где Ŷ t+1 – прогноз на следующий период t+1; Y t – данные для прогноза за текущий период t (например, продажи по месяцам); k - коэффициент сглаживания ряда , k задается вами вручную и находится в диапазоне от 0 до 1, 0 < k < 1; Ŷ t – значение прогноза на текущий период t Причем в первый период (месяц, день...) Ŷ 1 =Y 1 , т.е. Ŷ t в первый период равны продажам в этот период. Рассчитаем точность прогноза при k=0.9 (1 случай) и k=0.1 (2 случай). Для каждой модели определим: Ошибку модели - для этого в каждом периоде из фактических продаж вычитаем прогноз продаж на этот период. 4 Квадратическое отклонение - для каждого периода рассчитаем отношение квадрата ошибки модели к квадрату прогноза на этот период. Среднеквадратическое отклонение - среднее значение квадратических отклонений за весь анализируемый период. Точность прогноза = 1 - Среднеквадратическое отклонение. На Рисунке 1 показан пример расчетов прогнозирования продаж. Видно, что точность расчётов при к=0.9 равна 98,8%. Это высокая точность расчётов. Проверим при других значениях k. Рисунок 1.Расчёты по экспоненциальной модели в 1 случае. На Рисунке 2 показаны расчёты 2 случая для меньшего значения коэффициента k.Точность этих расчётов 94.7%. Точность ниже чем в первом случае. Этот прогноз менее точный по сравнению с первым. 5 Рисунок 2.Расчёты по экспоненциальной модели во 2 случае. На рисунке 3 показаны графики начальных данных, 1 модель и 2 модель. По графику можно определить какая модель более точна в расчётах и более совпадает с начальными данными. Данный способ используется для прогнозирования на 1 точку вперёд. Рисунок 3. График всех расчётов 6 1.2 Метод Хольта. Метод Хольта используется для прогнозирования временных рядов, когда есть тенденция к росту или падению значений временного ряда. А также для рядов, когда данные есть не за полный цикл, и сезонность еще не выделить (например, за неполный год для прогноза по месяцам). Если временной ряд имеет тенденцию к росту или падению, то вместе с оценкой текущего уровня ряда стоит выделить тренд. Для управления уровнем и наклоном в модели Хольта вводится 2 коэффициента сглаживания – коэффициент сглаживания ряда и тренда. Рассчитываем экспоненциально-сглаженный ряд: L t =k*Y t + (1-k)*(L t-1 -T t-1 ), где L t – сглаженная величина на текущий период; k – коэффициент сглаживания ряда; Y t – текущие значение ряда (например, объём продаж); L t-1 – сглаженная величина за предыдущий период; T t-1 – значение тренда за предыдущий период; Определяем значение тренда T t =b*(L t - L t-1 ) + (1-b)*T t-1, где T t – значение тренда на текущий период; b – коэффициент сглаживания тренда; L t – экспоненциально сглаженная величина за текущий период; L t-1 – экспоненциально сглаженная величина за предыдущий период; T t-1 – значение тренда за предыдущий период. 7 Делаем прогноз по методу Хольта. Прогноз на p периодов вперед равен Ŷt+p = Lt + p *Tt, где Ŷt+p – прогноз по методу Хольта на p период; Lt – экспоненциально сглаженная величина за последний период; p – порядковый номер периода, на который делаем прогноз; Tt – тренд за последний период. Результаты действий на Рисунке 4. Рисунок 4. Прогнозирование по методу Хольта на 3 точки. Рассчитаем значения прогнозной модели. Рассчитываем прогноз на 1 период вперед для каждого месяца, когда продажи нам известны. Прогноз на 1 период для анализа это есть значение экспоненциально-сглаженного ряда за предыдущий период + значение тренда за предыдущий период (Значение тренда мы не умножаем на p, т.к. прогноз делаем на 1 период, а в этом случае p=1). Рассчитаем ошибку модели = из фактических данных вычитаем прогноз 8 на этот период. Рассчитаем точность прогноза. Для этого определим отклонение ошибки модели от прогнозной модели = отношение ошибки модели в квадрате к фактическому значению в квадрате. Рассчитаем точность прогноза = единица минус среднее значение отклонений ошибок модели от прогнозной модели. Для подбора коэффициентов сглаживания ряда и тренда k и b, при которых прогноз будет максимально точным, нам необходимо последовательно перебрать все значения k и b в диапазоне от 0 до 1 и найти такое сочетание, при котором точность прогноза будет максимальна приближена к 100%. Мы хотим подобрать коэффициенты, которые будут делать максимально точный прогноз на 3 периода вперед (а не на 1, как в примере выше), то на анализируемых данных нам надо оценить модель с прогнозом на 3-й период. Для этого на входящих данных рассчитываем прогноз по следующей формуле: Ŷ t+3 = L t + 3 *T t . Далее определяем ошибку модели для каждого периода, среднеквадратическое отклонение и точность прогноза, как было описано выше. Последовательно перебирая коэффициенты k и b от 0 до 1 с шагом 0,1, находим такое сочетание, которое дает максимальное значение точности прогноза. Результаты действий на Рисунке 5. Рисунок 5. Полный прогнозируемый ряд 9 Расчётный ряд Ŷ t+p с моделью прогноза состоит из прогноза на 1 период и прогноза на 3 период, с добавлением 3 точек прогноза. Результатна Рисунке 6. Рисунок 6. График всех расчётов (Синий- продажи, красный – экспоненциально сглаженный ряд, зелёный прогнозируемый ряд) 10 2. Постановка задачи о прогнозировании временного ряда дифференциальным уравнением Классическим «материалом» для создания динамических (детерминистических) математических моделей являются дифференциальные уравнения. Мы ограничимся «обыкновенными» дифференциальными уравнениями (ОДУ), содержащими производные лишь по одной независимой переменной: dx dt = F(x) ,где x – зависимая (динамическая) переменная,t – независимая переменная. К уравнению первого порядка или эквивалентной ему системе уравнений: 𝒙 𝟏 ̇ = 𝑭 𝟏 (𝒙 𝟏 , … , 𝒙 𝑫 ), 𝒙 𝟐 ̇ = 𝑭 𝟐 (𝒙 𝟏 , … , 𝒙 𝑫 ), ... 𝒙 𝑫 ̇ = 𝑭 𝑫 (𝒙 𝟏 , … , 𝒙 𝑫 ). D- размерность модели. Имеется векторный временной ряд {𝑋 𝑖 } 𝑖=1 𝑁 , т.е. в каждый момент времени наблюдению доступны Dвеличин 𝑥 1 (𝑡 𝑖 ), … , 𝑥 𝐷 (𝑡 𝑖 ), однозначно определяющих состояние объекта(Рисунок 7). Выбор размерности модели D можно осуществить на основе предварительного анализа временного ряда . Модель имеет вид 𝒙 𝟏 ̇ = 𝒙 𝟐 , 𝒙 𝟐 ̇ = 𝒙 𝟑 , ... 𝒙 𝑫 ̇ = 𝑭 𝑫 (𝒙). Мы рассмотрим 2 метода. В методах рассмотрим динамические модели данных - системы с размерностью D равной 1 и исследуем полученные данные. 11 Рисунок 7. Схема процесса моделирования по временному ряду (Синяя линия – заданный временной ряд). Построение модели сводится к тому чтобы: По временному ряду {x i } рассчитать временной ряд { dx(t i ) dt }, Подобрать функцию F i , аппроксимирующие зависимости dx j dt от x, Проверить эффективность построенной модели. Производны вычисляются по формулам конечных разностей. Нахождение производных: dx i dt = x i+1 − x i−1 2∆t , (1.1) dx i1 dt = x 2 − x 1 ∆t , (1.2) 12 где ∆t –шаг по времени. Производная каждого временного ряда аппроксимируется некоторой заданной функцией с неизвестными коэффициентами. В данной работе использованы два метода подбора такой функции. Первый метод F j (x 1 , x 2 , … , x D ) = ∑ C j,l 1 ,l 2 ,…,l D ∏ X k l k , ∑ l k ≤ K j . (1.3) D k=1 D k=1 K j l 1 ,l 2 ,…,l D =0 В этом случае среднеквадратичная ошибка аппроксимации вычисляется как ε j 2 = 1 N train − m + 1 ∑ [ dx j (t i ) dt − F j (x i (t i ), x 2 (t i ), … , x D (t i ))] 2 = min. i (1.4) Во втором случае при подборе функции будем пользоваться результатами статьи [2]. В рассматриваемой статье в качестве аппроксимирующей функции использовалось x’ = g(t, x) = (∑ a i t i M i=0 ) x + b 0 + ∑ b j sin ( 2πjx θ + c j ) N j=1 . (1.5) Примеры других функций, которые можно использовать в аппроксимирующей функции F 𝑗 показаны в Таблице 1. 𝛼 sin 𝑥 + 𝛽 𝛼𝑒 𝑥 + 𝛽 𝛼𝑥 + 𝛽 𝛼𝑥 𝛽 + 𝑥 𝛼𝑥 2 + 𝛽ч + 𝛾 𝛼 ∗ 2 −𝑥/𝛽 𝛼 ln 𝑥 + 𝛽 𝛼𝛽 𝑥 𝛼 ln 𝑥 + 𝛽𝑥 + 𝛾 𝛼 + 𝛽𝛾 𝑥 𝛼 ln 𝑥 + 𝛽 ln 𝑥 + 𝛾 𝑒 𝛼+𝛽𝛾 𝑥 𝛼 sin 𝑥 + 𝛽𝑥 + 𝛾 ln 𝑥 + 𝛿 𝛼𝑒 𝛽𝑥 + 𝛾 Таблица 1. Пример функций для аппроксимации. 13 3. Метод решения Пусть дан временной ряд, значения которого описываются дискретной функцией x(t). Сначала вычислим производную этой функции по формулам (1.1) - (1.2). Затем, будем аппроксимировать полученную производную согласно уравнению (1.3) или уравнению (1.5) с разным набором функций из Таблицы 1. Уравнение (1.3) преобразуем для нашего случая с одним вектором значений. Получаем функцию аппроксимации (1 метод) Уравнение (1.4) преобразуем для нашего случая и выведем формулу для расчёта коэффициентов C. В нашем случае формула (1.4) примет вид здесь m - количество соседних точек. Для поиска коэффициентов C найдем экстремум функции ε 2 Перенося слагаемое без C вправо, приходим к системе уравнений C 0 (∑ x i n ) + … + C K (∑ x i n+K N i=0 ) = ∑(x ′ i x i n ) N i=0 N i=0 , C 0 (∑ x i K ) + … + C K (∑ x i 2K N i=0 ) = ∑(x ′ i x i K ). N i=0 N i=0 F i = ∑ C i ∗ x i K i=0 = C 0 + C 1 x 1 +...+ C K x K (2.1) ε 2 (C 0 , C 1 , … , C i ) = 1 N − m + 1 ∑(x i ′ − (C 0 + C 1 x 1 + ⋯ + C K x K )) 2 → min, N i ε 2 = 1 N−m+1 ∑ ( dx i dt − F i ) 2 → min. N i (2.2) Путём подробного разложения получаем: ∑(x′ i x i n ) N i=0 − ∑( N i=0 C 0 + C 1 x i + ⋯ + C K x K i ) = 0. 14 Введя обозначение A 𝑖𝑗 для выражений в скобках, получим C 0 A 00 + ⋯ + C K A 0K = B 0 , C 0 A K0 + ⋯ + C K A KK = B K ; Или, в матричной форме ( A 00 ⋯ A 0K ⋮ ⋱ ⋮ A K0 ⋯ A KK ) x ( C 0 … C K ) = ( B 0 … B K ) ; A × C = B; Ищем коэффициенты C через обратную матрицуA −1 С = A −1 B. (2.3) Для второго метода возьмём ещё одну аппроксимирующую функцию 𝐹. Пример данной функции описан в статье [2] из формулы (1.5). В нашем случае подбор функции F осуществляется путем анализа производной x(t), пример функции, использованной в расчетах 𝐹 = 𝐶 0 𝑥 3 + 𝐶 1 𝑥 + 𝐶 2 𝑥 2 + 𝐶 3 1 𝑥+1 + 𝐶 4 1 (𝑥+1) 2 + 𝐶 5 𝑠𝑖𝑛 𝑥 + 𝐶 6 𝑐𝑜𝑠 𝑥 + +𝐶 7 𝑙𝑛 𝑥 + 𝐶 8 𝑒 −𝑥 (2.4) После нахождения коэффициентов разложения C прогноз определяется как решение дифференциального уравнения 𝑥 ̇ = 𝐹(𝑥) с начальным условием в последней точке заданного временного интервала. Данное уравнение решается численно классическим методом Рунге — Кутты четвёртого порядка для нахождения приближённых значений k 1 = F(t i , x i ), 15 k 2 = F (t i + ∆t 2 , x i + ∆t 2 k1), k 3 = F (t i + ∆t 2 , x i + ∆t 2 k2), k 4 = F(t i + ∆t, x i + ∆tk3), x i+1 = x i + ∆t 6 (k 1 + 2k 2 + 2k 3 + k 4 ). (2.5) Вычисленные значения x(t) могут как использоваться для прогнозирования, так и для проверки точности алгоритма прогноза. Для этого у заданного временного рядя обрезается несколько значений, затем по усеченному ряду строится прогноз и полученные значения сравниваются с отброшенным хвостом ряда. В дальнейшем при обсуждении точности прогноза будем говорить об абсолютной ошибке (здесь мы понимает абсолютную разность значений прогноза и заданного временного ряда) и относительной (здесь мы будем понимать отношение абсолютной ошибки к значению временного ряда). 16 4. Результаты 4.1 Первый метод В работе [1] показано, что прогнозирование вектора временных рядов дифференциальным уравнением с правой частью в виде многочлена (уравнение (1.3)) хорошо работает для динамических систем. Действительно, в классических динамических системах (например, аттрактор Лоренца) правые части дифференциальных уравнений задаются в виде многочленов (обычно 1-3 порядка). Для аппроксимации более сложных функций, или временного ряда с шумами, встает вопрос о степени многочлена в формуле (1.3), необходимой для получения хорошей точности аппроксимации. Рассмотрим на следующих примерах точность аппроксимации производной x(t). Для расчётов были использованы формулы (2.1) - (2.5) и написана программа вычисления и построения графиков на языке программирования python. В подрисуночных подписях K - степень полинома F j , x – заданный модельный временной ряд. На Рисунках 8-9 показан пример расчетов для x = Sin T i 5 и K = 100 (Рисунок 1) и K = 50 (Рисунок 2). Синей линией показана вычисленная производная. Красной линией показана аппроксимация этой функции. Стоит заметить, что при больших значения К численное решение матричной задачи (2.3) может иметь проблемы. При вычислении обратной матрицы A могут возникать большие значения в коэффициентах матрицы, а также определитель матрицы может быть очень малым. Следовательно, вычисление значений коэффициентов разложения C может быть некорректным. В целом, на точность вычисления влияет не только число K, но и значения производной временного ряда. 17 Рисунок 8. Графики производной, аппроксимирующей функции и абсолютной ошибки (красная линия- аппроксимирующая функция, синяя – производная, зелёная – абсолютная ошибка),K = 100, x = Sin T i 5 Рисунок 9. Графики производной, аппроксимирующей функции и абсолютной ошибки (красная линия- аппроксимирующая функция, синяя – производная, зелёная – абсолютная ошибка),K = 50, x = Sin T i 5 18 На Рисунках 10-12 показан пример расчетов производной временного ряда, аппроксимирующей функции и абсолютной ошибки для случайным образом заданных значений временного ряда. Получено, что для точной аппроксимации необходимо брать многочлен не менее 50 степени. При этом точность аппроксимации может уменьшаться при увеличении количества точек на рассматриваемом интервале времени для заданного временного ряда. Можно заключить, что при использовании данного метода (аппроксимация многочленом) точность прогноза может быть неудовлетворительной в результате большого шанса получения вычислительной ошибки. Для использования этого метода необходимо аккуратно решать матричную задачу (2.3). Рисунок 10. Графики производной, аппроксимирующей функции и абсолютной ошибки (красная линия- аппроксимирующая функция, синяя – производная, зелёная – абсолютная ошибка), x = [4.0, 25.0, 1.0, 4.0, 64.0, 16.0, 81.0, 9.0], К=50. 19 Рисунок 11. Графики производной, аппроксимирующей функции и абсолютной ошибки (красная линия- аппроксимирующая функция, синяя – производная, зелёная – абсолютная ошибка),x = [4.0, 25.0, 1.0, 4.0, 64.0, 16.0, 9.0, 9.0, 1.0, 4.0, 64.0, 16.0, 81.0, 9.0], К=50. Рисунок 12. Графики производной, аппроксимирующей функции и абсолютной ошибки (красная линия- аппроксимирующая функция, синяя – производная, зелёная – абсолютная ошибка),x = [4.0, 25.0, 1.0, 4.0, 64.0, 16.0, 9.0, 9.0, 1.0, 4.0, 64.0, 16.0, 81.0, 9.0], К=10. 20 4.2 Второй метод В сравнении с первым методом, во втором методе аппроксимирующая функция задается в виде комбинаций разных функций. Причем, эти функции не обязательно должны являться многочленами. Для данного способа важным является анализ вычисленной производной от временного ряда. Подбором правильных функций можно добиться хорошей точности с малыми вычислительными затратами. Для реализации второго метода будем использовать вычислительную среду Matlab, в которой нам потребуется функция сглаживания данных. Зададим случайным образом временной ряд и добавим к нему шумы (Рисунок 13). Рисунок 13. Сгенерированный временной ряд. 21 Для расчёта значений используем формулы (1.1), (1.2) и (2.4). Вычисленное значение производной временного ряда показано на Рисунке 14. Заметим большое количество шумов на графике. Для аппроксимации производной предлагается использовать как саму производную, так и ее "сглаженную" форму. Для этого используем стандартное сглаживание в библиотеке Matlab, smoothinggaussian, c коэффициентом сглаживания «7». после сглаживания производная временного ряда примет вид, показанный на Рисунке 15. Рисунок 14. Производная от начальных данных. Используя классический метод Рунге — Кутты четвёртого порядка для поиска приближенных значений по итерационной формуле (2.5) построим прогноз. Результаты прогнозирования показаны на Рисунках 16-17. 22 Рисунок 15. График сглаженной производной. Рисунок 16.Графики исходного временного ряда и прогноза (красная линия – прогноз, синяя - начальные данные, зелёная - сглаженные начальные данные). 23 Рисунок 17. Графики исходного временного ряда и прогноза (красная линия – прогноз, синяя - начальные данные, зелёная - сглаженные начальные данные). Относительная ошибка данных между начальными данными и прогнозными данными не значительна, это видно на Рисунке 18. Значение ошибки в данном случае от 0.4% до 1%. Среднее значение ошибок равно 0.53%. Значение ошибки меньше чем в статье [2], алгоритм которой применяется в данной работе во втором методе. По полученным результатам можем сделать вывод о высокой точности второго метода. При этом, следует заметить, что для второго метода требуется время на анализ производной. 24 Следовательно, данный метод рекомендуется применять к долгосрочному планированию. Рисунок 18. Ошибка между данными (относительное отклонение прогноза и начальных данных). 25 Заключение Исследованы методы прогнозирования временных рядов с использованием дифференциальных уравнений. Временной ряд моделируется дискретной функцией x(t), значения которой равны значениям временного ряда в соответствующих точках по времени. Основная задача состоит в аккуратной аппроксимации производной полученной функции. Рассмотрены разные комбинации функций, аппроксимирующие эту производную. Предложены два метода. В первом методе производная аппроксимируется многочленом k-го порядка от x(t). Данный метод хорошо работает для прогнозирования динамических систем, которые описываются малым числом слагаемых. Для прогнозирования скалярного временного ряда со сложным поведением необходимо учитывать большую степень многочлена, что может приводить к вычислительным ошибкам. Во втором методе вычисленная производная сглаживается и подбирается система функций, которая, в целом, для каждого временного ряда может состоять из разных функций. Рассмотрены примеры аппроксимации на модельных данных. Показано, что второй метод может с хорошей точностью аппроксимировать временные ряды, и может использоваться для прогноза. Составлены и реализован алгоритм прогнозирования временных рядов дифференциальными уравнениями и построены соответствующие графики. 26 Список литературы 1. Б.П. БЕЗРУЧКО, Д.А. СМИРНОВ. Реконструкция обыкновенных дифференциальных уравнений по временным рядам. Учебно-методическое пособие, – Саратов: Издательство ГосУНЦ “Колледж”, 2000 – 46 с. 2. Xue, M., & Lai, C.-H. (2018). From time series analysis to a modified ordinary differential equation. Journal of Algorithms & Computational Technology, 85–90. https://doi.org/10.1177/1748301817751480 3. Подкорытова, О. А. Анализ временных рядов: учебное пособие для бакалавриата и магистратуры / О. А. Подкорытова, М. В. Соколов. — 2-е изд., перераб, и доп. — Москва: Издательство Юрайт, 2018. 4. Галочкин, В. Т. Эконометрика: учебник и практикум для бакалавриата и специалитета / В. Т. Галочкин. — Москва: Издательство Юрайт, 2019. 5. Филиппов А.Ф. Введение в теорию дифференциальных уравнений: Учебник. Изд. 2-е, испр. М., 2007. - 240 с 6. Карташев А.П., Рождественский Б.Л. Обыкновенные дифференциальные уравнения и основы вариационного исчисления. Изд. 2, исп. и доп. 1980. 7. Грибков Д.А., Грибкова В.В., Кравцов Ю.А., Кузнецов Ю.И., Ржанов А.Г. “Восстановление структуры динамической системы по временным рядам”, Радиотехника и электроника, 1994, Т. 39, В. 2, С. 269-277. 8. Kennel M.B., Brown R., Abarbanel H.D.I. “Determining embedding dimension for phase-space reconstruction using a geometrical construction”, Phys.Rev. A, 1992, Vol. 45, № 6, P. 3403-3411. 9. Кравцов Ю.А. “Случайность, детерминированность, предсказуемость”, Успехи физ. наук, 1989, Т. 158, № 1, С. 93-115. 10. Янсон Н.Б., Анищенко В.С. “Моделирование динамических систем по экспериментальным данным”, Изв. ВУЗов. Прикладная нелинейная динамика, 1995, Т. 3, № 3, С.112-121. 27 11. Crutchfield J.P., McNamara B.S. “Equations of motion from a data series”, Complex Systems, 1987, Vol. 1, P. 417-452. 12. Cremers J., Hubler A. “Construction of differential equations from experimental data”, Z.Naturforschung A, 1987, Vol. 42, P. 797-802. 13. Breeden J.L., Hubler A. “Reconstructing equations of motion from experimental data with unobserved variables”, Phys.Rev. A, 1990, Vol. 42, № 10, P. 5817- 5826. 14. Baake E., Baake M., Bock H.J., Briggs K.M. “Fitting ordinary differential equations to chaotic data”, Phys.Rev. A, 1992, Vol. 45, № 8, P. 5524- 5529. 15. Gouesbet G., Maquet J. “Construction of phenomenological models from numerical scalar time series”, Physica D, 1992, Vol. 58, P. 202-215. 16. Brown R., Rulkov N.F., Tracy E.R. “Modeling and synchronizing chaotic systems 44 from time-series data”, Phys.Rev. E, 1994, Vol. 49, № 5, P. 3784-3800. 28 ПРИЛОЖЕНИЕ А Листинг кода алгоритма 1.Python import numpy as np import math import matplotlib.pyplot as plt import numpy as np import xlrd import xlwt K = 50 A = np.zeros((K,K)) B = np.zeros((K,1)) read_book = xlrd.open_workbook('D:\Deeplom\Coopper1.xlsx', on_demand=True) read_sheet = read_book.sheet_by_index(0) T = [] row_number = read_sheet.nrows if row_number>0: for i in range(1,row_number): T.append(read_sheet.cell_value(i,2)) print('T = ',T) print('') x = [] for i in range(len(T)): x.append(T[i]**2/T[i]) print('x = ',x) print('') 29 dx1_dt = x[1]-x[0] dxi_dt = [dx1_dt] for i in range(1,len(x)-1): dxi_dt.append((x[i+1]-x[i-1])/2) i = i+1 dxi_dt.append(x[len(x)-1]-x[len(x)-2]) print('dxi_dt = ',dxi_dt) print('') for i in range(K): n = i for j in range(K): b=0 if j == 0: n += j else: n += 1 for z in range(len(x)): b += x[z]**n if z == len(x)-1: A[i][j] = b print('A = ',A) print('') for i in range(K): for j in range(1): b_2=0 30 for z in range(len(x)): b_2 += dxi_dt[z]*x[z]**i if z == len(x)-1: B[i][j] = b_2 print('B = ',B) print('') C = np.linalg.solve(A, B) print('C = ',C) print('') F = [] for i in range(1): for j in range(len(x)): F1 = 0 for z in range(K): F1 += C[z][i]*x[j]**z if z == (K-1): F.append(F1) print('Функция F =',F) print('') E = [] for i in range(len(T)): E.append(math.fabs(dxi_dt[i]-F[i])) print('Функция E =',E) print('') 31 fig = plt.figure() plt.grid() # включение отображение сетки plt.plot(F,'s-.r', label="F") plt.plot(dxi_dt,'o-.b', label="dx/dt") plt.plot(E,'p-.g', label="Ошибки") plt.legend() k1 = F[len(x)-1] xk2 = [] xk3 = [] xk4 = [] k2 = 0 k3 = 0 k4 = 0 for b in range(len(x)): xk2.append(x[b]+0.5*k1) k22 = len(x)-1 for k22 in range(len(x)): F1 = 0 for z in range(K): F1 += C[z][0]*xk2[k22]**(z+0.5) if z == (K-1): k2 = F1 for b in range(len(x)): xk3.append(x[b]+0.5*k2) for k22 in range(len(x)): 32 F1 = 0 for z in range(K): F1 += C[z][0]*xk3[k22]**(z+0.5) if z == (K-1): k3 = F1 for b in range(len(x)): xk4.append(x[b]+k3) for k22 in range(len(x)): F1 = 0 for z in range(K): F1 += C[z][0]*xk4[k22]**z if z == (K-1): k4 = F1 X_prog = x[len(x)-1]+(1/6)*(k1+2*k2+2*k3+k4) print('k1 =',k1) print('k2 =',k2) print('k3 =',k3) print('k4 =',k4) print('Прогноз =',X_prog) 33 ПРИЛОЖЕНИЕ Б Листинг кода алгоритма 2.Matlab tau = 0.05; N = 100; for i=1:1:N t(i) = (i-1) *tau; end tau=1; dx(1) = (x(2) - x(1)) /tau; dx(N) = (x(N) - x(N-1)) /tau; for i=2:1:N-1 dx(i) = (x(i+1) - x(i-1)) /(2 *tau); end DX1 = smoothdata(dx, 'gaussian' ,7); k = 9; for i = 1:1:k for j = 1:1:k A(i,j) = 0; for q = 1:1:N A(i,j) = A(i,j) + Fcalc(x(q),i,t(q)) *Fcalc(x(q),j,t(q)); end end B(i,1) = 0; for q = 1:1:N B(i) = B(i) + DX1(q) *Fcalc(x(q),i,t(q)); end end C = A\B; for q=1:1:N t(q) = (q-1) *tau; dxc(q) = 0; for i = 1:1:k dxc(q) = dxc(q) + C(i) *Fcalc(x(q),i,t(q)); end end t2=t; x2=x; 34 xx = smoothdata(x, 'gaussian' ,10); for q=N-6:1:N t2(q) = (q-1) *tau; x2(q) = x2(q-1); for i = 1:1:k vv(i) = C(i) *Fcalc(x2(q-1),i,t2(q)) *tau; x2(q) = x2(q) + C(i) *Fcalc(x2(q-1),i,t2(q)) *tau; end err_abs(q) = abs(x2(q) - xx(q)); err_otl(q) = err_abs(q)/xx(q); end plot(t,x,t2,x2) function F = Fcalc(x,i,t) if i == 1 F = x^3; end if i == 2 F = x; end if i ==3 F = x^2; end if i == 4 F = 1/(x+1); end if i == 5 F = 1/(x+1)^2; end if i == 6 F = sin(x); end if i == 7 F = cos(x); end if i == 8 F = log(x); end if i == 9 F = exp(-x); end end |