Теоретические сведения по численным методам решения задач аппроксимации

Министерство науки и высшего образования Российской Федерации Федеральное государственное бюджетное образовательное учреждение высшего образования «Пермский национальный исследовательский политехнический университет» Электротехнический факультет Кафедра «Микропроцессорных средств автоматизации» Н. В. Андриевская ЧИСЛЕННЫЕ МЕТОДЫ Для студентов заочной формы обучения Пермь 2020 УДК 62-52 Рекомендовано к изданию редакционно-издательским советом Пермского национального исследовательского политехнического университета Рецензенты кандидат технических наук, С. В. Березняков (ОАО «СТАР») доктор технических наук, профессор Южаков А. А. (Пермский национальный исследовательский политехнический университет) \ Н. В. Андриевская Учебное пособие – конспект лекций «Численные методы» /Н. В. Андриевская. – Пермь: Изд-во Перм. нац. исслед. политехн. ун-та, 2020 - 31 с. Содержит основные теоретические ведения по численным методам решения задач аппроксимации, обработки экспериментальных данных и решения обыкновенных дифференциальных уравнений. Для студентов специальностей 15.03.04 «Автоматизация технологических процессов и производств в машиностроении и энергетике», 13.03.02 «Электропривод и автоматика», 13.03.02 «Электроснабжение», 27.03.04 «Управление в технических системах» УДК 62-52 © Н. В. Андриевская, 2020 © Пермский национальный исследовательский политехнический университет, 2020 2 Содержание 1. Интерполяция функций ...................................................................................... 5 1.1. Постановка задачи аппроксимации ............................................................. 5 1.2. Интерполяционные формулы для равноотстоящих точек. ....................... 6 1.2.1. Конечные разности .................................................................................. 7 1.2.2. Первая интерполяционная формула Ньютона ..................................... 8 1.2.3. Вторая интерполяционная формула Ньютона...................................... 9 1.2.4. Применение интерполяционных формул Ньютона ........................... 10 1.3. Интерполяционные формулы для неравноотстоящих точек. Формула Лагранжа.............................................................................................................. 11 1.4. Погрешности интерполяционных формул ............................................... 12 1.5. Обратная интерполяция. ............................................................................. 13 2. МЕТОДЫ ОБРАБОТКИ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ. ............... 16 2.1. Постановка задачи обработки экспериментальных данных ................... 16 2.2. Выбор критерия согласия ........................................................................... 18 2.2.1. Среднеквадратичный критерий ........................................................... 18 2.2.2. Минимаксный критерий, или критерий Чебышева ........................... 19 2.2.3. Вероятностно-зональные критерии ..................................................... 19 2.3. Выбор аппроксимирующей функции ........................................................ 19 2.4. Выбор узловых точек .................................................................................. 20 2.5. Метод наименьших квадратов ................................................................... 21 2.5.1. Применение МНК для определения параметров линейной функции ........................................................................................................................... 22 2.5.2. Применение МНК для определения параметров квадратичной функции ............................................................................................................ 24 2.5.3. Применение МНК для определения параметров степенной функции ........................................................................................................................... 25 2.5.4. Применение МНК для определения параметров показательной функции ............................................................................................................ 27 2.5.5. Применение МНК для определения параметров логарифмической функции ............................................................................................................ 28 2.5.6. Применение МНК для определения параметров гиперболической функции ............................................................................................................ 29 2.5.7. Применение МНК для определения параметров дробно-линейной функции ............................................................................................................ 30 3 2.5.8. Применение МНК для определения параметров дробнорациональной функции ................................................................................... 31 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ. ......................................................... 33 3.1. Постановка задачи решения обыкновенных дифференциальных уравнений ............................................................................................................ 33 3.2. Методы решения обыкновенных дифференциальных уравнений ........ 35 3.3. Численные методы решения обыкновенных дифференциальных уравнений ............................................................................................................ 36 3.3.1. Метод Эйлера. ........................................................................................ 37 3.3.2. Первая улучшенная формула Эйлера. ................................................. 39 3.3.3. Вторая улучшенная формула Эйлера. ................................................. 42 3.3.4. Метод Рунге-Кутта. ............................................................................... 45 4 1. Интерполяция функций 1.1. Постановка задачи аппроксимации На практике при решении многих прикладных задач возникает проблема замены одной функции другой по разным причинам: исходная функция сложна для исследования, исходная функция задана таблично, исходная сложна для реализациях в моделях и другое. Замена одной функции другой называется аппроксимация. Методы, осуществляющие аппроксимацию, называются аппроксимационными. Постановка задачи аппроксимации определяется следующим образом: Пусть дана таблично определенная функция: № 1 … n x x x 1 … x n y y y 1 … y n Требуется построить функцию F ( x ) , которая в точках x принимала бы i значения F ( x )  y . Такая задача называется задачей аппроксимации. i i Таблично заданные точки называются узловыми точками Рис. 1.1. Аппроксимация точечно (таблично) заданной функции. 5 Аппроксимация очень востребована в различных задачах автоматизации технологических процессов, в теории автоматического управления, в задачах прогнозирования и других. Задачи аппроксимации подразделяются на задачи интерполяции и экстраполяции. В задачах интерполяции аппроксимирующая функция используется для исследования функции внутри интервала наблюдения [ x , x ] : 0 n (1.1) F ( x ), x  [ x , x ], x  x i  0, n 0 n i В задачах экстраполяции построенная аппроксимирующая функция используется для исследования функции вне интервала наблюдения [ x , x ] : 0 n (1.2) F ( x), x [ x , x ], x  x  x  x 0 n n Наиболее распространенной и удобной для исследования является формула полинома: F ( x)  a  a x  a x2  a xn (1.3) 0 1 2 n Поэтому на практике наибольшее распространение получили методы построения аппроксимирующей функции в виде полинома (многочлена) (1.3) Изначально задача аппроксимации с помощью полиномов применялась для определения значений функции внутри рассматриваемого интервала исследования, поэтому и разработанные алгоритмы (формулы) аппроксимации носили и носят название интерполяционные формулы В зависимости от расположения узловых точек внутри интервала наблюдения методы построения полиномов разделяются на две группы:  методы (формулы) для равноотстоящих точек: (1.4) x  x  const  h i 1 i  методы (формулы) для неравноотстоящих точек: (1.5) x  x  const i 1 i 1.2. Интерполяционные формулы для равноотстоящих точек. К формулам для равноотстоящих точек относятся: первая и вторая интерполяционные формы Ньютона, центральные интерполяционные формулы Гаусса, Бесселя, Стирлинга. Данные интерполяционные формулы для равноотстоящих точек использую конечные разности. 6 1.2.1. Конечные разности x [ x , x ] с 0 n постоянным шагом – приращением аргумента x  h , тогда конечная разность 1-го порядка определяется как: Пусть таблично задана функция y на интервале (1.6) y  y  y i  0, n  1 i i 1 i Конечная разность 2-го порядка имеет вид:  2 y  y  y i  0, n  2 (1.7) i i 1 i Конечная разность n-го порядка имеет вид: n y  n  1 y  n  1 y i  0 (1.8) i i 1 i Нетрудно заметить, что конечная разность максимального порядка – одна. Конечные разности с точки зрения вычисления и дальнейшего применения удобнее всего оформлять в виде горизонтальной таблицы конечных разностей (табл. 1.1): Табл. 1.1. x y y x y y x 1 y 1 y 1 x n y n 2 y 2 y 2 y 1 n y n y ПРИМЕР Пусть задана функция в виде таблицы x 2 4 6 8 10 y 1 2 8 16 24 30 Тогда в соответствии с (1.6-1.8) горизонтальная таблица конечных разностей будет иметь вид: 7 № x y y 1 2 3 4 5 2 4 6 8 10 1 2 8 16 25 36 1 6 8 9 11 2 y 5 2 1 2 3 y -3 -1 1 4 y 2 2 5 y Следует отметить следующие особенности, как конечных разностей, так и горизонтальной таблицы конечных разностей: 1. Максимально «возможная» (расчетная) конечная разность имеет порядок n  1 , где n - количество узловых точек. 2. Таблица имеет диагональный вид, так количество конечных разностей с каждым порядком уменьшается. Рассмотрим основные интерполяционные формулы 1.2.2. Первая интерполяционная формула Ньютона Первая интерполяционная формула Ньютона имеет вид: q( q  1) 2 q( q  1)( q  2) 3  y   y  2! 3! q( q  1) ( q  n  1) n   y n! P ( x )  y  qy  n , (1.9) где n - максимальный порядок полинома; i y - конечные разности i -го порядка; q - параметр первой интерполяционной формулы Ньютона (1ИФН), который определяется по формуле: xx 0, q (1.10) h Нетрудно понять, что: 8 xx x  ( x  h) x  x 0 1  1  h h h xx x  ( x  2) x  x 0 2 2 q2  h h h q 1  q  n 1  , (1.11) xx x  ( x  (n  1)h ) x  x 0  ( n  1)  n  h h h Тогда с учетом (1.11) 1ИФН можно представить в следующем виде: (x  x ) ( x  x )( x  x ) 2 1  y  P ( x)  y  y  n h 2! h 2 ( x  x )( x  x ) ( x  x ) 1 n  1 n y  n !hn , (1.12) Следует отметить, что порядок аппроксимирующего полинома не превышает величины n  1 , где n - количество узловых точек. 1.2.3. Вторая интерполяционная формула Ньютона Вторая интерполяционная формула Ньютона (2ИФН) имеет вид: q(q  1) 2 q(q  1)(q  2) 3  y   y  n2 n3 2! 3! q(q  1) (q  n  1) n   y n! P ( x)  y  qy  n n n 1 , (1.13) где n - максимальный порядок полинома; i y - конечные разности i -го порядка; j q - параметр второй интерполяционной формулы Ньютона (2ИФН), который определяется по формуле: xx n, (1.14) q h По аналогии с (1.11) можно получить: 9 xx x  ( x  h) x  x n 1  n n 1  h h h xx n2 q2 , h q 1  q  n 1  (1.15) xx h Тогда с учетом (1.14) и (1.15)2ИФН можно представить в следующем виде: P ( x)  y  n n (x  x ) (x  x )( x  x ) n  1 y n 1 n  2 2 y   n 1 n2 2 h 2!h (x  x ) (x  x ) n n 1 0  y  n !h n , (1.16) 1.2.4. Применение интерполяционных формул Ньютона Обе интерполяционные формулы могут использоваться как для получения общего вида полинома (1.3), так и для решения частной задачи: расчета конкретного значения функции в точке, не являющейся узловой. В этом случае под точкой ( x , y ) в 1ИФН принимается ближайшая 0 0 узловая точка. Во 2ИФН соответственно ближайшая узловая точка – это ( x , y ) n n ПРИМЕР Пусть дана таблично-заданная функция. Рассчитать значение в точке x5 y № x 1 2 3 2 4 6 8 2 7 15 4 10 31 10 Построим горизонтальную таблицу конечных разностей № x y y 2 y 3 y 4 y 2 2 7 -3 11 1 4 2 6 1 8 2 6 8 7 9 3 4 8 10 15 31 16 Ближайшая узловая точка x  4 y  2 . Параметр xx 0  5  4  0.5 q h 2 Тогда значение в точке x  5 по 1ИФН рассчитывается следующим образом: q( q  1) 2 q( q  1)( q  2) 3 P ( x  5)  y  qy   y   y  4 2! 3! 0.5( 0.5) 0.5( 0.5)( 1.5)  2  0.5  6  1 8  5.375 2! 3! Нетрудно заметить, что чем больше слагаемых применяется в 1ИФН, тем точнее аппроксимация. Целесообразность 2ИФН становится очевидной, если возникнет задача расчета значения в точке, находящейся в конце таблицы. 1.3. Интерполяционные формулы для неравноотстоящих точек. Формула Лагранжа Если узловые точки по аргументу функции расположены неравномерно, то интерполяционные формулы Ньютона применять нельзя. В этом случае применяются интерполяционные формулы для неравноотстоящих точек, например формула Лагранжа: ( x  x )( x  x ) ( x  x )( x  x ) (x  x ) 1 i 1 i2 n , L ( x)   y n i ( x  x )( x  x ) ( x  x )( x  x ) ( x  x ) i 0 i 0 i 1 i i 1 i i  2 i n (1.17) n 11 ПРИМЕР Пусть задана функция в виде таблицы Пример: № 1 2 3 x 1 2 5 y 2 3 12 147 ( x  x )( x  x )( x  x ) ( x  x )( x  x )( x  x ) 1 2 3 2 3  L ( x)  y  y  3 0 ( x  x )( x  x )( x  x ) 1 ( x  x )( x  x )( x  x ) 0 1 0 2 0 3 1 0 1 2 1 3 ( x  x )( x  x )( x  x ) ( x  x )( x  x )( x  x ) 1 3 y  1 2  y  2 ( x  x )( x  x )( x  x ) 3 ( x  x )( x  x )( x  x ) 2 0 2 1 2 3 3 1 3 2 3 3 ( x  1)( x  2)( x  5) ( x  0)( x  2)( x  5) ( x  0)( x  1)( x  5) 3  12  (0  1)(0  2)(0  5) (1  0)(1  2)(1  5) (2  1)(2  1)(0  5) ( x  0)( x  1)( x  2) 147  (5  0)(5  1)(5  2) 2  x3  x 2  x  2 1.4. Погрешности интерполяционных формул Естественно, любая интерполяционная формула дает погрешность, оценка которой осуществляется по остаточному члену: R ( x)  f ( x)  P ( x) n n (1.18) Формула остаточного члена для 1ИФН имеет вид: R ( x)  n q( q  1)( q  2) ( q  n) n  1  y (n  1)! (1.19) или R ( x)  n ( x  x )( x  x ) ( x  x ) n  1 1 n  y n  1 h ( n  1)! 12 (1.20) Для 2-й интерполяционной формулы соответственно: R ( x)  n q( q  1)( q  2) ( q  n ) n  1  y ( n  1)! (1.21) или R ( x)  n ( x  x )( x  x ) (x  x ) n 1 n n 1 0  y n  1 h ( n  1)! (1.22) Если вид функции, которая аппроксимируется, известен, то погрешность 1ИФН имеет вид q( q  1)( q  2) ( q  n ) [n  1] R ( x)  hn  1 f ( ) , n (n  1)! (1.23) где f [n  1] ( ) - ( n  1) -я производная исходной функции в точке  , обеспечивающей наибольшее значение данной производной на интервале [ x , x ]. 0 n После подстановки (1.10), (1.11) формула (1.23) для 1ИФН имеет вид: ( x  x )( x  x ) ( x  x ) [n  1] 1 n f (1.24) R ( x)  ( ) n (n  1)! Для 2-й интерполяционной формулы: соответственно q( q  1)( q  2) ( q  n ) [n  1] R ( x)  hn  1 f ( ) , n ( n  1)! R ( x)  n ( x  x )( x  x ) n n 1 (n  1)! или ( x  x ) [n  1] 0 f ( ) (1.25) (1.26) Погрешность формулы Лагранжа оценивается по формуле: max f (n  1) ( x ) x xx n R ( x)  0 ( x  x )( x  x ) n 1 (n  1)! (x  x ) n 1.5. Обратная интерполяция. Пусть дана таблично определенная функция: № x x 13 y y (1.27) 1 … n x 1 … x n y 1 … y n Задача обратного интерполирования заключается в определении значения аргумента функции x по заданному значению функции y . Предположим, что заданная табличная функция f ( x ) монотонна и значение y содержится между узловыми точками. Тогда заменяя y интерполяционным полиномом Ньютона на основе первой интерполяционной формулы: y0  2 y0  n y0 y  f ( x)  y0  q q(q  1)  ...  q(q  1)...(q  n  1) (1.28) 1! 2! n! Значение аргумента функции вычисляется по формуле: x  x0  qh (1.29) Определить q можно из формулы (1.28): yy 2 y n y 0 0 q (q  1)  ...  0 q(q  1)...(q  n  1) q (1.30) y 2!y n!y Формула (1.30) является основой для итерационной процедуры: yy 2 y n y 0 q (q  1)...(q  n  1) (1.31) q   q (q  1)  ...  i 1 i i i y 2!y n!y i i За начальное приближение q принимается: y y (1.32) q  y Расчет q выполняется до тех пор, пока не выполнится условие: i (1.33) q q  i 1 i При выполнении условия (1.32) q  qi 1 и расчет значения аргумента определяется по формуле (1.29) Известно, что первая интерполяционная формула дает хорошие результаты для значений функции, находящихся в начале таблицы. Для значений функции, находящихся в конце таблицы применяется вторая интерполяционная формула Ньютона. Формула для обратного интерполирования рассчитывается аналогично и имеет вид: 14 y y q  i  1 y 2 y n y n  n q ( q  1)...  n q ( q  1)...( q  n  1) (1.34) i i i i i 2! y n ! y n 1 n 1 n 1 Следует обратить особое внимание, что пара точек ( x , y ) (1ИФН) и 0 0 ( x , y ) (2ИФН) – это ближайшие к исходной точке значения. n n Если аргументы функции расположены неравномерно, то применяется интерполяционная формула Лагранжа. Для задачи обратного интерполирования она имеет вид: n ( y  y )...( y  y )( y  y )...( y  y ) i 1 i 1 n (1.35.8) x  x i ( y  y )...( y  y )( y  y )...( y  y ) i0 i i i 1 i i 1 i n 15 2. МЕТОДЫ ОБРАБОТКИ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ. 2.1. Постановка задачи обработки экспериментальных данных Пусть в результате измерений в процессе опыта получена таблица некоторой зависимости: x y x 1 y 1 x y n n Необходимо найти формулу, выражающую эту зависимость аналитически. Необходимо отметить, что такая постановка задачи характерна и для задач интерполяции. Однако, следует отметить, что различие методов интерполяции и методов обработки экспериментальных данных заключается в том, что методы интерполяции предполагают построение полинома определенного порядка, который проходит через узловые точки (рис 2.1), а методы обработки экспериментальных данных предполагают построение функции заданного вида, график которой проходит как можно ближе к узловым точкам (рис. 2.2). Отсутствие требования точного совпадения функции, построенной методом экспериментальных данных, объясняется тем, что полученные в результате эксперимента данные обладают погрешностью (погрешностью измерения, погрешностью округления и т.д.), поэтому требование точного прохождения функции через узловые точки является некорректным. 16 Рис. 2.1. Интерполяция функций y x Рис. 2.2. Обработка экспериментальных данных. При решении задачи обработки экспериментальных данных необходимо ответить на следующие вопросы: 1. Выбор метода на основе реализации критерия согласия: 2. Выбор узловых точек 3. Выбор функции 17 2.2. Выбор критерия согласия Это самый важный вопрос. Решение этого вопроса дает ответ решения задачи наилучшего приближения. Что означает математически наилучшее приближение? Это означает выбрать критерий согласия, который является функцией невязки узловых точек и значениями аппроксимирующей функции: J  J (F ( x )  y ) , i i (2.1) где x , y - узловые точки; i i F ( x ) - значение аппроксимирующей функции в узловой точке x ; i i J - функция невязки (ошибки) Выбор наилучшей функции осуществляется по минимуму этого критерия: J  J ( F ( x )  y )  min , i i (2.2) Существуют три наиболее широко распространенных критерия согласия:  -среднеквадратичный;  - минимаксный или Чебышева;  - вероятностно-зональный. 2.2.1. Среднеквадратичный критерий Среднеквадратичный критерий предполагает минимизацию суммы квадратов ошибки в узловых точках: n J   ( F ( x )  y )2 , i i i 1 (2.3) Среднеквадратичный критерий позволяет получить сглаживание кривой, то есть позволяет отфильтровать зашумленные данные, не требуя никакой дополнительной информации о шумовых характеристиках помех. 18 2.2.2. Минимаксный критерий, или критерий Чебышева Минимаксный критерий Чебышева определяет функцию невязки как: J  max F ( x )  y i i (2.4) Если применение среднеквадратичного критерия уменьшает среднеквадратичную ошибку, при этом допуская отдельные большие ошибки, то чебышевское приближение - минимаксное - уменьшает экстремальную наибольшую ошибку. Поэтому этот критерий используется, когда необходимо при аппроксимации избежать больших ошибок. Минимаксный критерий также не использует дополнительную информацию о шумовых характеристиках помех. 2.2.3. Вероятностно-зональные критерии К данным критериям относится целая группа критериев. Данные критерии используют (требуют) дополнительную информацию о шумовых характеристиках объекта:  обобщенный метод наименьших квадратов - ковариационные матрицы шума;  максимальное правдоподобие - распределение вероятностей и т.д. Поэтому область применения данных критериев ограничена требованиями вышеприведенной априорной информацией 2.3. Выбор аппроксимирующей функции Практически вид приближающей функции F можно определить следующим образом. Строится точечный график функции, заданной таблично, а затем проводиться плавная кривая, по возможности 19 наилучшим образом отражающая характер расположения точек. По полученной таким образом кривой устанавливается вид приближающей функции (обычно из числа простых по виду аналитических функций). В качестве функции приближения могут быть использованы следующие элементарные функции: 1. линейная: F ( x)  ax  b 2. (2.5) квадратичная: F ( x)  ax2  bx  c 3. (2.6) степенная: F ( x)  axb 4. (2.7) показательная: F ( x)  aebx 5. (2.8) логарифмическая: F ( x)  a ln x  b 6. гиперболическая: F ( x)  7. a b x (2.10) дробно-линейная: F ( x)  8. (2.9) 1 ax  b (2.11) дробно-рациональная: F ( x)  x ax  b (2.12) 2.4. Выбор узловых точек В принципе это вопрос статистики, а именно той области, которая называется "планированием эксперимента". 20 На практике узловые точки заданы внешними обстоятельствами или используются равноотстоящие точки. Если же существует возможность выбора точек, то выбор осуществляется по формуле Чебышева: ba ba  , 2 2 i 2i  1  , i  0, n . где    cos i 2n  2 Узлы интерполяции, выбранные по формуле Чебышева, x  i (2.5) располагаются не равномерно, а группируются концов отрезка [a, b] На практике наибольшее распространение получил метод обработки экспериментальных данных с использованием среднеквадратичного критерия согласия. Из-за вида функции невязки данный критерий получил название Метод наименьших квадратов. 2.5. Метод наименьших квадратов Близость функции к таблично заданной функции осуществляется методом наименьшего квадрата (МНК), который предполагает минимизацию суммы квадратов ошибки в узловых точках: n J   ( F ( x )  y )2 i i i 1 где y - значение исходной функции в точке x (табличное); i i F ( x ) -значение аппроксимирующей (приближающей) функции. i Нетрудно заметить, что в аппроксимирующих функциях неизвестными являются параметры (2.5-2.12) a, b, c . Рассмотрим задачу в общем виде. 21 Приближающая функция имеет общий вид F ( x, a, b, c) , тогда сумма квадратов может быть определена как функция неизвестных параметров a, b, c : n Ф(a, b, c)   ( y  F ( x , a, b, c)) 2 i i i 1 (2.13) Функция Ф(a, b, c) примет минимальное значение при выполнении следующих условий – условий экстремума:  дФ(a, b, c) 0  дa   дФ(a, b, c) 0  дb   дФ(a, b, c) 0  дc  (2.14). С учетом (2.13) и правилами дифференцирования сложных функций система уравнений (2.14) примет вид:  n '   ( yi  F ( xi , a, b, c)  Fa ( xi , a, b, c)  0 i  1  n  '   ( yi  F ( xi , a, b, c)  Fb ( xi , a, b, c)  0 i  1  n   ( y  F ( x , a, b, c)  F ' ( x , a, b, c)  0 i c i i  1 i  (2.15) Решение системы (2.15) относительно неизвестных параметров a, b, c и определяет параметры соответствующей F ( x, a, b, c) 2.5.1. Применение МНК для определения параметров линейной функции Рассмотрим подробно применение метода МНК для линейной функции. Система уравнений (2.15) для линейной функции будет иметь вид: 22  n   ( yi  axi  b)  xi  0 i  1   n   ( yi  axi  b)  0  i 1 (2.16) Или после преобразований:  x  y  a  x 2  b x  0 i i  i i i i i    yi  a  xi  nb  0  i i (2.17) Разделив каждое уравнение на n – количество измерений или узловых точек, получаем:  1 2 b 1 x  1 x  y a  x     n i ni i ni i i  i  1 1  a x b  y  ni i ni i (2.18) Введем обозначения: 1 x  M x ni i 1 y M y ni i (2.19) 1 y x M xy ni i i 1 2 x  M ni i x2 С учетом обозначений (2.19) система уравнений (2.18) примет вид: a  M  b  M  M x xy  x2  aM b  M  x y  (2.20) Разрешив данную систему уравнений относительно неизвестных параметров а и b, определим искомую функцию F ( x, a, b)  ax  b 23 Результаты вычислений удобнее и нагляднее оформить в виде таблицы (2.1) Табл. 2.1 y x x 1 x x y 1 y 2 y n M x 2 n M y xy x2 F ( x) x y 1 1 x2 1 F(x ) 1 x y 2 2 x2 2 F(x ) 2 x y n n x2 n F(x ) n M M xy e  y  F ( x) e 1 e 2 e n x2 e2 e2 1 e2 2 e2 n E Значение E определяет близость аппроксимирующей функции к исходной. E определяется по следующей формуле: n 2  e i i  1 Е n (2.21) Естественно, чем меньше E, тем аппроксимирующая функция ближе к исходной функции. 2.5.2. Применение МНК для определения параметров квадратичной функции Система уравнений (2.15) для квадратичной функции будет иметь вид: 24  ( y  a  x 2  b  x  c)  x 2  0 i i i i   2   ( yi  a  xi  b  xi  c)  xi  0    ( y  a  x 2  b  x  c)  0 i i i  (2.22) После преобразований и подставки обозначений аналогичных (2.19) система уравнений примет вид:  M 4  a  M 3  b  M 2  c  M 2 x x x y  x   M 3  a  M 2  b  M x  c  M xy x x   M a  M b  c  M x y  x2 (2.23) Разрешив данную систему уравнений относительно неизвестных параметров a, b, с а, b и с, определим искомую функцию F ( x, a, b, с)  ax2  bx  c . Результаты вычислений для квадратичной функции оформляются в виде таблицы (2.2): Табл. 2.2 x y M M x y xy x2 x3 x4 M M M M xy x2 x3 x4 x2 y M F ( x) e x2 y 2.5.3. Применение МНК для определения параметров степенной функции Степенная функция имеет вид: 25 e2 E F ( x, a, b)  axb Если прологарифмировать данное уравнение, то можно получить следующее: ln( F ( x))  ln(a)  b  ln( x) Введя новые переменные: Ab B  ln(a ) u  ln( x)  ( x)  ln( F ( x)) (2.24) можно получить следующее линейное уравнение:  (u)  Au  B (2.25) Тогда система уравнений для новых параметров A, B будет иметь вид: A M  B  M  M u u  u2  A M  B  M  u   (2.26) Определив параметры А и В линейной функции (2.25), можно определить параметры степенной функции: b A a  eB (2.27) Результаты вычислений оформляются в виде таблицы (2.3) Табл. 2.3 x y u  ln(x)   ln( y) M M u  26 u u2 M M u u2 F ( x) e e2 E Следует заметить, что данный алгоритм справедлив только для положительных значений x и y. 2.5.4. Применение МНК для определения параметров показательной функции Показательная функция имеет вид: F ( x, a, b)  aebx Если прологарифмировать данное уравнение, то можно получить следующее: ln( F ( x))  ln(a)  bx Введя новые переменные: Ab B  ln(a )  ( x)  ln( F ( x)) (2.28) можно получить следующее линейное уравнение:  ( x)  Ax  B (2.29) Тогда система уравнений для новых параметров A, B будет иметь вид: A M  B  M  M x x  x2  A M  B  M  x   (2.30) Определив параметры А и В линейной функции (2.30), можно определить параметры степенной функции: b A a  eB (2.31) Результаты вычислений оформляются в виде таблицы (2.4) 27 Табл. 2.4 y x M x   ln( y) M  x x2 M M x F ( x) x2 e e2 E Следует заметить, что данный алгоритм справедлив только для положительных значений y. 2.5.5. Применение МНК для определения параметров логарифмической функции Показательная функция имеет вид: F ( x, a, b)  aln( x)  b Если ввести новые переменные: u  ln( x) (2.32) можно получить следующее линейное уравнение: F ( x)  au  b (2.33) Тогда система уравнений (2.20) будет иметь вид: aM  bM  M u uy  u2   aM u  b  M y  (2.34) Решая систему уравнений (2.34) относительно a, b , можно определить параметры логарифмической функции: Результаты вычислений оформляются в виде таблицы (2.5) 28 Табл. 2.5 y x M u  ln( x) M x M u u2 uy uy M F ( x) u2 e e2 E Следует заметить, что данный алгоритм справедлив только для положительных значений x . 2.5.6. Применение МНК для определения параметров гиперболической функции Гиперболическая функция имеет вид: F ( x, a , b )  a b x Если ввести новые переменные: u 1 x (2.35) можно получить следующее линейное уравнение: F ( x)  au  b (2.36) Тогда система уравнений (2.20) будет иметь вид: aM  bM  M u uy  u2   aM u  b  M y  (2.37) Решая систему уравнений (2.37) относительно a, b , можно определить параметры логарифмической функции: Результаты вычислений оформляются в виде таблицы (2.6) 29 Табл. 2.6 y x M x 1 u x M M u u2 uy uy M F ( x) u2 e e2 E Следует заметить, что данный алгоритм справедлив только для x  0 . 2.5.7. Применение МНК для определения параметров дробно-линейной функции Дробно-линейная функция имеет вид: F ( x, a , b )  1 ax  b Если ввести новые переменные:  ( x)  ln( F ( x)) (2.38) можно получить следующее линейное уравнение:  ( x)  ax  b (2.39) Тогда система уравнений (2.20) будет иметь вид: aM  bM  M x x  x2   aM x  b  M  (2.40) Решая систему уравнений (2.40) относительно a, b , можно определить параметры логарифмической функции: Результаты вычислений оформляются в виде таблицы (2.7) 30 Табл. 2.7 y x M  1 y M x  x x2 M M x F ( x) x2 e e2 E Следует заметить, что данный алгоритм справедлив только для y  0 . 2.5.8. Применение МНК для определения параметров дробнорациональной функции Дробно-рациональная функция имеет вид: F ( x, a , b )  x ax  b Введя новые переменные: Ab Ba 1 u x  (2.41) 1 F ( x) можно получить следующее линейное уравнение:  (u)  Au  B (2.42) Тогда система уравнений для новых параметров A, B будет иметь вид: 31 A M  B  M  M u u  u2  A M  B  M  u   (2.43) Определив параметры А и В линейной функции (2.43), можно определить параметры дробно-рациональной функции: b A aВ (2.44) Результаты вычислений оформляются в виде таблицы (2.8) Табл. 2.8 x y u M u 1 x  M 1 y  u u2 M M u F ( x) e u2 Следует заметить, что данный алгоритм справедлив только для положительных значений x  0, y0 32 e2 E 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ. 3.1. Постановка задачи решения обыкновенных дифференциальных уравнений Обыкновенные дифференциальные уравнения (ОДУ) встречаются при описании движения системы взаимодействующих материальных носителей (механике, химической кинетике, электрических цепях, сопротивлении материалов и многих других явлениях жизни). Ряд важнейших задач сводится к задачам с обыкновенными дифференциальными уравнениями, а именно теория автоматического управления базируется на ОДУ. Простейшим ОДУ является уравнение 1-го порядка: dy  y  f ( x, y ) dx (3.1) ОДУ n -го порядка имеет вид: dny  y n  f ( x, y, y, y,..., y (n  1) dx n (3.2) Однако любое уравнение n -го порядка можно свести к системе из n уравнений 1-го порядка:  dy ( x)  dx  y1( x)   dy1( x)  y2 ( x )   dx ...   dyn  1( x)  f ( x, y0 , y1,..., yn  1)  dx (3.3) Таким образом, решение дифференциального уравнения n-го порядка сводится к решению системы дифференциальных уравнений 1-го порядка из n уравнений. 33 Можно систему уравнений представить в матричном виде. Для этого вводятся обозначения:  y1    y  Y   2 ...  y   n  f1    f  F   2 ...  f   n Тогда систему дифференциальных уравнений (3.3) можно переписать в векторно-матричном виде: dY ( x)  F ( x, Y ( x)) dx (3.4) Известно, что система n -го порядка имеет множество решений, которые в общем случае зависят от постоянных С, общее количество которых -n. Общее решение системы ОДУ может быть записано в виде: Y  Y ( x, C) , (3.5)  c1    c  где C   2  - некоторые константы. ...  c   n Для определения констант необходимо учитывать дополнительные условия. В зависимости от выбора дополнительных условий определены следующие типы задач для обыкновенных дифференциальных уравнений:  задача Коши;  краевая задача;  задача на собственные значения. В задаче Коши дополнительными условиями являются Y ( X 0 )  Y0 (3.6) Поэтому задачу Коши называют задачей с начальными условиями. В краевой задаче кроме начальных условий (3.6) заданы и конечные условия 34 Y ( X n )  Yn (3.7) В задаче на собственные решение Y зависит не только от X , но и от параметров k k  1, l : Y ( x)  F ( x,Y , 1, 2 ,..., l ) (3.8) Для разрешения системы уравнений (3.8) число дополнительных условий должно быть соответственно n  l . Этими условиями являются, как начальные условия (3.6), так и дополнительные l условий для параметров k : i  i* (3.9) Функции Y ( x),  , где k  1, n и k r r  1, l , удовлетворяющие всем условиям (3.6), (3.9) называются собственными дифференциальными или собственными значениями задачи. На практике, а особенно в задачах автоматизации и управления наибольшее распространение получили задача с начальными значениями, или задача Коши, поэтому далее будут представлены методы решения задач данного класса. 3.2. Методы решения обыкновенных дифференциальных уравнений Обыкновенные дифференциальные уравнения могут быть решены следующими методами:  аналитическими  приближенными  графическими  численными Аналитические методы дают решение в виде аналитического выражения (3.1). 35 Приближенные методы также дают аналитическое решение, но при этом функция f ( x, y )  P( x, y ) заменяется на другую функцию, более удобную для интегрирования. Графические методы дают приближенное решение в виде графика. Численные методы дают частное решение для определенного интервала x   x0 , xn    3.3. Численные методы решения обыкновенных дифференциальных уравнений Численные методы применяются только к корректно поставленным (хорошо обусловленным) задачам. к таким, у которых малое изменение начальных значений, приводит к малому изменению интегральных кривых. Рассмотрим пример плохо обусловленной задачи. Дано обыкновенное дифференциальное уравнение: dy  y  x. dx Необходимо решить задачу Коши для 0  x  100 при начальных условиях y(0)  1 Общее решение имеет вид: y( x, C )  1  x  Ce x При y(0)  1, значение C  0 и решение y(100)  101. Однако при малом изменении начальных условий y(0)  1,000001  C  106 решение в точке x  100 : y(100)  2.7 1037 , то есть имеет место плохо обусловленная задача. На практике наибольшее распространение получили следующие численные методы решения дифференциальных уравнений:  метод Эйлера;  первый улучшенный метод Эйлера; 36  второй улучшенный метод Эйлера;  метод Рунге-Кутта. 3.3.1. Метод Эйлера. В основе метода Эйлера (метод ломаных) лежит идея графического построения решения дифференциальных уравнений, однако этот метод дает одновременно и способ нахождения искомой функции в численной (табличной) форме. Пусть дано дифференциальное уравнение: dy  f ( x, y ) dx (3.10) с начальными условиями: y ( x0 )  y0 (3.11) Выбрав достаточно малый шаг h, строится система равноотстоящих точек: xi  x0  ih (i  1, n) (3.12) В методе Эйлера приближенные значения y ( xi )  yi вычисляются последовательно по формулам: yi  1  yi  hf ( xi , yi ) (3.13) При этом искомая интегральная кривая y  y( x) , проходящая через точки М i ( xi , yi ) , заменяется ломанной М 0 М1 ... М n с вершинами М i ( xi , yi ) ; каждое звено этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной y '  f ( x, y) , которая проходит через точку M i : 37 кривой уравнения y M0 M1 M2 M3 x x0 x1 x2 x3 Рис. 3.1. Построение ломаной интегральной линии Эйлера Несмотря на простоту в реализации, метод редко применяется на практике, поскольку обладает малой точностью. Погрешность каждого нового шага систематически возрастает. Существует модификации метода, повышающие его точность, - методы Эйлера-Коши – первая и вторая улучшенные формулы. ПРИМЕР. Решить ОДУ dy  3e2 x  y с начальными условиями dx y(0)  1 на интервале [0;1] с шагом интегрирования h  0.2 с помощью формулы Эйлера (ФЭ). 38 i 1 x1  x0  h  0  0.2  0.2 y1  y0  hf ( x0 , y0 )  1  0.2(3e 20  1)  1.400 i2 x2  x1  h  0.2  0.2  0.4 y2  y1  hf ( x1, y1)  1.400  0.2(3e 20.2  1.400)  2.015 Результаты приведены в таблице ФЭ y x 0.0 1.000 0.2 1.400 0.4 2.015 0.6 2.947 0.8 4.350 1.0 6.452 3.3.2. Первая улучшенная формула Эйлера. Для исходных данных (3.10, 3.11) решение в каждой точке xi определяется по формуле: yi  1  yi  hf ( x1  1/ 2 , yi  1/ 2 ) , (3.14) h где xi  1 / 2  xi  2 h yi  1 / 2  yi  f ( xi , yi ) . 2 Геометрически это означает, что отрезок ломанной между точками M i , M i  1 заменяется на два отрезка  M i , M i  1 / 2  ,  M i  1 / 2 , M i  1  . Направление первого отрезка совпадает с направлением интегральной 39 кривой в точке ( xi , yi ) , а направление второго отрезка определяется направлением, интегральной кривой в вспомогательной точке ( xi  1 / 2 , yi  1 / 2 ) . (рис. 3.2) y yiЭ M M i 1 M i 1 2 f ( x, y ) i xi 1/2 xi x xi 1 Рис. 3.2. Построение ломаной интегральной линии по первой улучшенной формуле Эйлера При применении первой улучшенной формулы Эйлера уменьшение погрешности достигалось за счет промежуточной точки ( xi  1 / 2 , yi  1 / 2 ) . Однако вспомогательную точку можно выбрать другим способом, как это сделано во второй улучшенной формуле Эйлера. 40 ПРИМЕР. Решить ОДУ dy  3e2 x  y с начальными условиями dx y(0)  1 на интервале [0;1] с шагом интегрирования h  0.2 с помощью первой улучшенной формулы Эйлера (1УФЭ). x0  0 i 1 y0  1 h 0.2 x0  1 2  x0   0   0.1 2 2 h 0.2 20 y0  1 2  y0  f ( x0 , y0 )  1  (3e  1)  1.200 2 2 x1  x0  h  0  0.2  0.2 y1  y0  hf ( x0  1 2 , y0  1 2 )  1  0.2(3e 20.1  1.200)  1.493 i2 h 0.2 x1  1 2  x1   0.2   0.3 2 2 h 0.2 20.2 y1  1 2  y1  hf ( x1, y1)  1.493  (3e  1.493)  1.791 2 2 x2  x1  h  0.2  0.2  0.4 y2  y1  hf ( x1  1 2 , y1  1 2 )  1.493  0.2(3e 20.3  1.791)  2.228 Результаты приведены в таблице 41 ФЭ 1 УФЭ y y x 0.0 1.000 1.000 0.2 1.400 1.493 0.4 2.015 2.228 0.6 2.947 3.324 0.8 4.350 4.960 1.0 6.452 7.400 3.3.3. Вторая улучшенная формула Эйлера. Для исходных данных (3.10, 3.11) решение в каждой точке xi определяется по формуле: yi  1  yi  h f ( xi , yi )  f ( xi* 1, yi* 1) 2 , (3.15) где xi* 1  xi  h yi* 1  yi  hf ( xi , yi ) Геометрически это означает, что определяется направление интегральной кривой в исходной точке ( xi , yi ) и во вспомогательной точке ( xi  1, yi  1) , а в качестве окончательного направления ломаной берется среднее этих направлений (рис. 3.3) 42 y yiЭ M M i 1 f ( x, y ) yi2УЭ i x xi xi 1 Рис. 3.3. Построение ломаной интегральной линии по второй улучшенной формуле Эйлера Первая и вторая улучшенные формулы Эйлера имеют одинаковый класс точности, так как направление перехода из точки ( xi , yi ) в точку ( xi  1, yi  1) определяется с помощью одной вспомогательной точки. Очевидно, что чем больше вспомогательных точек тем, тем более точное решение дает метод. Поэтому метод Рунге-Кутта, основанный на четырех и более вспомогательных точках, имеет более высокий класс точности. ПРИМЕР. Решить ОДУ dy  3e2 x  y с начальными условиями dx y(0)  1 на интервале [0;1] с шагом интегрирования h  0.2 с помощью второй улучшенной формулы Эйлера (2УФЭ). 43 x0  0 i 1 y0  1 x*  x0  h  0  0.2  0.2 1 y*  y0  hf ( x0 , y0 )  1  0.2(3e 20  1)  1.400 1 x1  x0  h  0  0.2  0.2 h 0.2 y1  y0  ( f ( x0 , y0 )  f ( x* , y* ))  1  ((3e 20  1)  (3e 20.2  1.400))  1.508 2 2 1 1 i2 x*  x1  h  0.2  0.2  0.4 2 y*  y1  hf ( x1, y1)  1.508  0.2(3e 20.2  1.508)  2.101 2 x2  x1  h  0.2  0.2  0.4 h 0.2 y2  y1  ( f ( x1, y1)  f ( x* , y* ))  1.508  ((3e 20.2  1.508)  (3e 20.4  2.101)  2 2 2 2  2.262 Результаты приведены в таблице ФЭ 1 УФЭ 2УФЭ y y 0.0 1.000 1.000 1.000 0.2 1.400 1.493 1.508 0.4 2.015 2.228 2.262 0.6 2.947 3.324 3.358 0.8 4.350 4.960 5.058 1.0 6.452 7.400 7.553 x y 44 3.3.4. Метод Рунге-Кутта. Строго говоря, и метод Эйлера, и первая, и вторая улучшенная формулы Эйлера относятся к семейству методов Рунге-Кутта. Метод Рунге-Кутта q -го порядка имеет вид: q y ( x  h)  y ( x)   pi ki (h) , i 1 где при фиксированных (3.16) значениях 2 , 3, ..., q ; p1, p2, p3,..., pq ; ij некоторых (0  j  i  q) параметров последовательно вычисляются: k1(h)  hf ( x, y ) k2 (h)  hf ( x   2h, y   21k1(h)) ... kq (h)  hf ( x   q h, y   q1k1(h)   q 2k2 (h)  ...   qq kq (h)) (3.17) Наибольшее применение на практике получил метод Руте-Кутта 4-го порядка: 1 yi  1  yi  (k1  2(k2  k3 )  k4 ) , 6 (3.18) k1  hf ( x, y ) k h k2  hf ( x  , y  1 ) 2 2 k h k3  hf ( x  , y  2 ) 2 2 k4  hf ( x  h, y  k3 ) (3.19) где Порядок метода определяется количеством вспомогательных точек, определяющих направление перехода из точки ( x , y ) в точку ( x . i i i  1, yi  1) В методе Рунге-Кутта 4-го порядка направление перехода определяется с помощью 4-х вспомогательных точек. 45 Метод Рунге-Кутта имеет ряд важнейших достоинств:  высокая точность;  явная схема вычислений yn1 за конечное количество шагов;  возможность применения переменного шага интегрирования; ПРИМЕР. Решить ОДУ dy  3e2 x  y с начальными условиями dx y(0)  1 на интервале [0;1] с шагом интегрирования h  0.2 с помощью формулы РунгеКутта 4-го порядка (ФРК) 46 x0  0 i 1 y0  1 k1  hf ( x0 , y0 )  0.2(3e 20  1)  0.400 k h 0.400 k2  hf ( x0  , y0  1 )  0.2(3e 2(00.1)  (1  )  0.493 2 2 2 k h 0.493 k3  hf ( x0  , y0  2 )  0.2(3e 2(00.1)  (1  )  0.484 2 2 2 k4  hf ( x0  h, y0  k3 )  0.2(3e 2(00.2)  (1  0.484)  0.598 x1  x0  h  0  0.2  0.2 k  2(k2  k3 )  k4 0.400  2(0.493  0.484)  0.598 y1  y0  h 1  1  0.2  1.492 6 6 i2 k1  hf ( x1, y1)  0.2(3e 20.2  1.492)  0.597 k h 0.597 k2  hf ( x1  , y1  1 )  0.2(3e 2(0.20.1)  (1.492  )  0.735 2 2 2 k h 0.735 k3  hf ( x1  , y1  2 )  0.2(3e 2( 00.1)  (1.492  )  0.721 2 2 2 k4  hf ( x1  h, y1  k3 )  0.2(3e 2(00.2)  (1.492  0.722)  0.911 x2  x1  h  0.2  0.2  0.4 k  2(k2  k3 )  k4 0.597  2(0.735  0.721)  0.911 y2  y1  h 1  1.492  0.2  2.229 6 6 47 Результаты приведены в таблице ФЭ 1 УФЭ 2УФЭ ФРК y x y y 0.0 1.000 1.000 1.000 1.000 0.2 1.400 1.493 1.508 1.492 0.4 2.015 2.228 2.262 2.229 0.6 2.947 3.324 3.358 3.330 0.8 4.350 4.960 5.058 4.973 1.0 6.452 7.400 7.553 7.427 Можно сравнить результаты аналитическим y решением решений численными методами с данного дифференциального уравнение. Аналитическое решение имеет вид: y  e2x ФЭ 1 УФЭ 2УФЭ ФРК Аналитическое решение x y y y y y 0.0 1.000 1.000 1.000 1.000 1.00 0.2 1.400 1.493 1.508 1.492 1.492 0.4 2.015 2.228 2.262 2.229 2.226 0.6 2.947 3.324 3.358 3.330 3.320 0.8 4.350 4.960 5.058 4.973 4.953 1.0 6.452 7.400 7.553 7.427 7.389 Определить близость полученного решения к точному значению с помощью оценок: 48 n 2  (YTi  YMi ) E  i 1 , n 2  YTi 1 (3.20) где y - точное решение ОДУ в точке x ; Тi i yМi - численное решение уравнение соответствующим методом в точке xi ; n - количество точек на интервале [ x ; x ] n ФЭ 1УФЭ 2УФЭ ФРК Аналитическое решение x y y y y y 0.0 1.000 1.000 1.000 1.000 1.00 0.2 1.400 1.493 1.508 1.492 1.492 0.4 2.015 2.228 2.262 2.229 2.226 0.6 2.947 3.324 3.358 3.330 3.320 0.8 4.350 4.960 5.058 4.973 4.953 1.0 6.452 7.400 7.553 7.427 7.389 0.265 0.003 0.046 0.010 E 49