CC BY
33
УДК 517.912:519.6 ББК В161.6:В192.1
В.Х. ФЕДОТОВ, Н.В. НОВОЖИЛОВА
ВЫСОКОТОЧНЫЙ МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
Ключевые слова: математическое моделирование, численное интегрирование дифференциальных уравнений, высокоточный численный метод решения задачи Коши, динамические процессы в областях науки и техники, численные методы, компьютерные вычисления, сложные технические системы.
Точность численного интегрирования дифференциальных уравнений зависит от малоизученных законов алгебры приближенных вычислений - погрешностей методов, дискретизации, устойчивости разностных схем и др. Как следствие, на каждой итерации решается задача, отличная от предыдущей, что приводит к накоплению ошибки расчета (эффект «сползания»). В работе описан и апробирован простой высокоточный численный метод решения задачи Коши без использования высших производных и с минимальным эффектом «сползания».
Обыкновенные дифференциальные уравнения (ОДУ) часто используются в качестве математических моделей сложных динамических процессов в различных областях науки и техники. При этом, как правило, эти уравнения не решаются аналитически, а исследуются с помощью различных численных методов.
Рассмотрим задачу Коши для системы ОДУ
У = ёуШ = Лу(0, 0, у(к) = Уо, (1)
где = (уь у2,..., ук) - вектор зависимых переменных; у0 - начальное значение; ^ - независимая переменная (время); ¿0 - начальный момент времени (начало эволюции); . [ - известная (заданная) функция.
Задача Коши в общем виде не имеет точного решения в квадратурах [3]. Численные методы решения основаны на различных разностных схемах, представляющих собой конечное приближение бесконечного ряда, аппроксимирующего функцию у(0. Чем выше порядок приближения, тем, как правило, ниже погрешность метода (увеличения порядка приближения не обязательно снижает погрешность). Все применяемые на практике численные методы без используют высшие производные и характеризуются погрешностью порядка ~0(в5), где в << 1.
Целью данной работы является построение простого численного метода решения задачи Коши с минимальной погрешностью без использования высших производных.
Численные методы решения систем ОДУ включают одношаговые (методы Эйлера, методы Рунге - Кутта), многошаговые (методы Адамса), адаптивные с автоматическим выбором шага и др., модифицированные (прогноза и коррекции), а также их комбинации [2]. Предлагаемый ниже метод можно отнести к модификациям одношаговых методов Эйлера или методам прогноза и коррекции.
Одномерные системы. Рассмотрим вначале одномерный вариант. Од-ношаговые методы являются простейшим способом численного интегрирования ОДУ, основанным на аппроксимации непрерывной и бесконечно-
дифференцируемой функции y(t), конечным числом p членов бесконечного ряда Тейлора в малой окрестности начального приближения [3]
y(to+h) *y(to) + y'(to)h + y" (to)h2/2!+У" (to)h3/3! +.... + y(p)(to)h p /p!, (2) где h - шаг (предполагается постоянным).
Это предположение не существенно для целей данной работы и не нарушает общности анализа, так как нас интересует точность вычислений, но не их скорость. Погрешность формулы (2) (остаток порядка p) определяется бесконечным, вообще говоря, рядом
Rph) = y(P+1)(to)hp+1 /(p+1)! + y(p+2)(to)hp+2/(p+2)! + ... - o(hP)- 0(hp+l), (3) где O(hp+1) обозначает бесконечно малую порядка p + 1, т.е. < Chp+1 (C - 1), а o(hp) - бесконечно малую более высокого порядка, чем p, т.е. < shp1 (в << 1).
Основной вклад в эту погрешность в случае общего положения вносит первое слагаемое, зависящее от (p + 1)-й производной в точке to и шага h. Сумма остальных образует погрешность порядка малости p + 2 и т.д.
При разработке алгоритмов численных вычислений аналитические соотношения представляют в виде дискретных итерационных формул. Подставляя в (2) вместо у' ее известное значение f из (1) с учетом (3) при p = 1, получим обобщенную разностную формулу первого порядка:
Уп+1 = Уп + hy'(tn) + R1(h) *Уп + h fy„, t„)+ RUn, n = o, 1, 2,.. , (4) где yn - значение на предыдущей итерации; yn+1 - значение на следующей итерации; R1n - локальная ошибка на текущей итерации; o <h <<1 - шаг разностной схемы.
В этом случае основной вклад в погрешность вносит вторая производная в точке tn. При R1(h) = o из (4) следует классический метод Эйлера первого порядка, для которого локальная ошибка максимальна R1n = R^^a - O(h2).
На практике погрешность итерационных формул вида (4) зависит от малоизученных законов алгебры приближенных вычислений (нарушение коммутативности и др.) [5] - локальной ошибки метода (на данной итерации), ошибок дискретизации (округлений, отсечений), устойчивости разностных схем и др. Их следствием является то, что вычисление yn+1 для каждой итерации фактически означает решение задачи Коши, отличной от предыдущей, так как представляет собой линию семейства кривых y(t, С) с началом в точке yn. Этот эффект «сползания» может привести к накоплению больших глобальных ошибок R(h) = YRp(h) на всем интервале интегрирования. Например, глобальная ошибка метода Эйлера на порядок выше локальной R(h) - O(h) [5]. Бесконечный остаток Rp(h) в соотношении (3) является ошибкой метода, а его конечное приближение R1n в (4) является ошибкой дискретизации. Они близки, но не равны и обе «немного» отличаются от неизвестной (истинной, точной) локальной ошибки En решения исходной задачи Коши (1), так как не известно ее точное решение y(tn). Критерий высокой точности численного решения задачи Коши запишем в виде задачи нелинейной минимизации локальной ошибки на интервале [tm+1, tn]:
R1,n * En = lyn+1-y(tn) l ^ min. (5)
Соотношения (4), (5) за счет выбора разных способов оценки остатка позволяют конструировать алгоритмы численного интегрирования различной точности. Например, оценка остатка в форме Шлёмильха - Роша использует значение
второй производной внутри отрезка (t^t). Ее частными случаями являются оценки Лагранжа - Коши [4]
Rln * ДЛагранж = /'(6„^2/2 * ^Коши = у"^^)^ - O(h\ (6)
где e„ = t„+a„h, 0 < a„ < 1.
Позднее Пеано получил оценку Дпеано = P«h2/2 - o(h2), где Р„ = pn(h) ^ 0 при h ^ 0 - неизвестная функция [4. С. 295], которая еще неопределеннее. Примерно в это же время Рунге [1] предложил, по-видимому, наиболее простую оценку
Rln * Дунге = С^у p-1(e„) - O(hP), (7)
где С - константа.
Эти оценки используют значения высших производных в неопределенных точках и являются разновидностями теорем о среднем [3], которые недостаточно конструктивны. Если высшие производные можно найти из условий задачи (1) через полные дифференциалы у"=ff+ft, у"'=fvyf2 + f2 + ftt, то an , en , P„ остаются неизвестными1.
Найдем более точные и конструктивные формы оценки R1,n, не использующие высшие производные. Рассмотрим поведение у() вблизи точки tn. Локально любую непрерывную и дифференцируемую функцию можно считать монотонной (не обязательно строго). Возможные случаи:
1) у(1) выпукла вниз (лежит над касательной);
2) у(1) выпукла вверх (лежит под касательной);
3) у(^ имеет перегиб (лежит по разные стороны от касательной). Рассмотрим их подробнее.
Случай 1. Пусть функция у() на интервале [tn, tn+1] лежит над касательной (выпукла вниз, у'(tn) Ф 0, y"(tn) > 0). Тогда метод Эйлера на каждой итерации занижает значение yn+1 по сравнению с точным значением y(tn) на неизвестную величину E = |yn+1 - y(tn)| * R1n. Оценим ее через первые производные, значения которых у' (tn) = f(tn) для задачи Коши известны в любой предыдущей точке. Очевидно, что грубой оценкой является разность между двумя последовательными приближениями R1n < |yn+1 - yn| = h\y' (tn)|. Модули учитывают знак производной sign у' (tn) Ф 0, который может быть любым. Улучшим оценку за счет неизвестного значения первой производной на правом конце tn+1 интервала R1,n < h\y' (tn)| < h\y' (tn+1)|. Эту оценку можно улучшить за счет другого значения первой производной в точке en = tn + an(t„+1 — tn) = tn + anh = antn+1 + (1 - an)tn, 0 < an < 1, расположенной немного левее an = 1 - в правого конца интервала R1,n< h\f(en)\. Понятно, что улучшение можно продолжать, сдвигаясь еще левее. При этом внутри интервала должно существовать неизвестное оптимальное значение en*, обеспечивающее минимум ошибки En (разновидность теорем о среднем) [3]
h min|f(tn+an*h)| < Rоптим * h |у'(e„*)\ < h max|f(tn+an*h)|, (8)
где en* = tn + an*h - оптимальное значение; 0 < an* < 1 - коэффициент оптимальной точности.
1 Примечание. Теорема о среднем [1, 4]. Если функция/,$) непрерывна и дифференцируема на (а, Ь), то существует такое число 0, что/Ь) -/а)=/ '(0)(Ь - а), где 0 = а + а(Ь - а), 0 < а < 1.
Оценка (8) в отличие от (6) использует только первую производную. Наиболее близкой к (8) является оценка Рунге (7), использующая производные на единицу меньшего порядка, включая и первый. Отличие (8) от (7) состоит в том, что она вообще не зависит от производных высших порядков.
Чтобы сделать оценку (8) конструктивной, зададим алгоритм вычисления у'(0„*). Проще всего это сделать с помощью линейной у'(9„*) * ап*у'(¿„) + у„, квадратичной или кубической интерполяции [3]. Для линейной интерполяции при у„ = 0 соотношение (4) с учетом (8) запишется в виде оптимизированного алгоритма Эйлера
Уп+1 * Уп + Ьуп' + коптим = Уп + Ьуп' + ап-Нуп' = (9)
= уп + к(1 + ап*)уп', п = 0, 1, 2,.....
Критерием точности выберем разность между интерполированным значением и нужным числом членов остатка. Для определения нужного числа членов остатка Др(к) зададим желаемую точность ~0(Нк). Тогда критерий (5) запишется в виде полиномиального уравнения, корни которого дают оптимальные значения ап*:
Еп* Доптим - ад *к[ап*уп'| - ^у(Р+1)п кР/(Р+1)!] = 0,р = 1, 2, 3, ..., к-1. (5') Для линейной интерполяции, соответственно, получим единственное оптимальное значение
ап* = Ту(р+Г)„ЬР/[(р+1)!уп 1], Р = 1, 2, 3, ..., к-1. (10)
Чтобы обеспечить точность ~0(к2), полагаем к = 2, т.е. р = 1, и из (10) получим ап* = у"п к/(2|уп'|). Чтобы обеспечить точность ~0(к3), полагаем к = 3 и из (10) получим ап* = [у"пк/2 + у'"п к2/6]/|уп'| и т.д. Отметим, что критерий точности (5) в отличие от самого метода (9) использует производные высших порядков. Однако практика показывает, что достаточно применить его один раз на первом (самом грубом) шаге. В численных вычислениях критерием точности обычно считают абсолютную величину разности двух последовательных приближений Дп = |уп+г - уп|, которая легко вычисляется, но отличается от истинной погрешности Еп и может неверно оценивать реальные ошибки вычислений. Критерий (5' ) в этом смысле, конечно, объективнее.
Возможны и другие способы определения оптимального значения ап*. Если разложить в ряд не само решение, а производную у'(^+к) * у'(¿0) + + у"(^о)к + у"'(¿0)к2/2! + ... и вновь применить линейное приближение, то получим рекуррентное соотношение /п+\ *./п + ./' п к для вычисления первой производной решения на конце интервала. С его помощью получаем уравнение касательной в конце интервала. Затем параллельно переносим ее в начало интервала. Полученная линия представляет собой секущую, проходящую через начало интервала, но выше его конца. Точка ее пересечения (¿*, у*) с вертикалью t = tn+1 дает уточнение верхней оценки (8) Лоптим * * ку' (9п*)| < т*)1
Простейшими методами определения ап* являются перебор и различные эвристики. Например, для простейшей эвристики «среднего» асредн = 1/2 и оценка (8) запишется
Д1,п * Дсреднее = 0,5к|у 'п I . (11)
Для эвристики «золотого» сечения 1/а = а/(1 - а) существуют два симметричных значения азолот ~ 0,382, 1 - азолот ~ 0,618 и два варианта записи оценки (8) Лзолото1 = 0,382/у '„!, ^золото2 = 0,618/у '„!• Сравним найденные модификации метода Эйлера с точными значениями, оценкой Лагранжа и методами Эйлера первого (Эйлер 1) и второго (Эйлер2) порядка.
Пример 1. Решим задачу Коши для уравнения у' = ay при t0 = 0, у0 = 1 на отрезке [0, /г] для разных /. Ее точным решением является функция у(^ = eat, которая выпукла вниз (лежит над касательной) и зависит от числа а. При а < 0 решение и ошибки затухают, при а > 0 - неограниченно возрастают. Пусть а = -10 (устойчивый случай). Расчет по методу Эйлера1 дает сильно заниженное значение уЭйлер = 0,0. Расчет по методу Эйлера2 требует вычисления второй производной у" =/у[+= а2у, но дает завышенное значение. Оценка Лагранжа также требует расчета второй производной и неизвестного значения ап, например ап = 1/2. Оптимизированные модификации метода Эйлера проще и точнее (табл. 1-4).
Таблица 1
Точность модификаций метода Эйлера для уравнения у' = ау при Н = 0.1, а = -10 (у0 = 1)
г о X & я а ю 1 я а ю •а я а ю X ю 1 ю 2 ¡2 ю -а 1 =13 ю
о н £0 13 О и * 13 О £0 13 О а О 13 О 13 О 13 О 5 а X О 13 О
0,00 1,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,000000
0,10 0,37 0,00 -0,37 0,25 -0,12 0,50 0,13 0,50 0,13 0,38 0,01 0,62 0,25 0,37 0,000021
0,20 0,14 0,00 -0,14 0,06 -0,07 0,25 0,11 0,25 0,11 0,15 0,01 0,38 0,25 0,14 0,000015
0,30 0,05 0,00 -0,05 0,02 -0,03 0,13 0,08 0,13 0,08 0,06 0,01 0,24 0,19 0,05 0,000008
0,40 0,02 0,00 -0,02 0,00 -0,01 0,06 0,04 0,06 0,04 0,02 0,00 0,15 0,13 0,02 0,000004
0,50 0,01 0,00 -0,01 0,00 -0,01 0,03 0,02 0,03 0,02 0,01 0,00 0,09 0,08 0,01 0,000002
0,60 0,00 0,00 0,00 0,00 0,00 0,02 0,01 0,02 0,01 0,00 0,00 0,06 0,05 0,00 0,000001
Глобальная ошибка -0,58 -0,25 0,40 0,40 0,04 0,95 0,00005
По строке 2 видно, что наименьшую локальную ошибку на первом шаге дают методы золотого сечения £золото1 = 0,01 ~ 0(/ ) и оптимальный Ешт™ = 0,00002 ~ 0(/ ), полученные при значении аоптим ~ 0,3679, найденном бисекцией. Оптимальное значение по формуле (10) для р = 1 равно ап* = -а//2 = 1/2, что совпадает с модификацией среднего.
Для р = 2
ап* = -а//2 + а2/2/6 = 1/2 -1/6 = 1/3.
Для р = 3
ап* = -а//2 + а2/2/6 - а3/3/24 = 1/3 + 1/24 = 9/24 « 0,375.
Для р = 4
ап* = -а//2 + а2/2/6 - а3/3/24 + а4/4/120 = 1/3 + 1/24 = = 9/24 - 1/120 = 44/120 « 0,366667.
Для р = 5
ап* = -а//2 + а2/2/6 - а3/3/24 + а4/4/120 - а5/5/720 = = 44/120 + 1/720 =265/720 = 0,368056.
Для р = 6
а„* = -аН/2 + а2Н216 - а3Н3/24 + а4Н4/120 - а5Н5/720 + а6Н6/5040 = = 265/720 - 1/5040 = 0,367857 «а,™. Это значение обеспечивает нужную точность.
Ошибка метода «Среднее» £среднее = 0,4 ~ О(Н), который также не использует вторую производную, сравнима с ошибками по Лагранжу и методу Эйлера второго порядка. Наименьшую глобальную ошибку дал также оптимальный метод, даже без пересчета оптимального значения на каждом шаге. Если на каждом шаге пересчитывать значения аоптим, то результаты будут еще точнее. Для а = 10 (неустойчивый случай) оценки приведены в табл. 2.
Таблица 2
Точность модификаций метода Эйлера для уравнения У = ау при Н = 0.1, а = +10 (у0 = 1)
о X ю и ю •а ю = ю 1 ю 2 ¡2 ю -а п ?« ю
о н ;= £0 13 О и я Ч Е О ;= £0 13 О & О Е О о эт 13 О о эт 13 О ¡е а с О 13 О
0,00 1,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,000000
0,10 2,72 2,00 -0,72 2,25 -0,47 2,50 -0,22 2,50 -0,22 2,38 -0,34 2,62 -0,10 2,72 0,000718
0,20 7,39 4,00 -3,39 5,06 -2,33 6,25 -1,14 6,25 -1,14 5,67 -1,72 6,85 -0,54 7,39 0,003905
0,30 20,09 8,00 -12,09 11,39 -8,69 15,63 -4,46 15,63 -4,46 13,52 -6,57 17,94 -2,14 20,10 0,015926
0,40 54,60 16,00 -38,60 25,63 -28,97 39,06 -15,54 39,06 -15,54 32,19 -22,40 46,98 -7,62 54,66 0,057730
0,50 148,41 32,00 -116,41 57,67 -90,75 97,66 -50,76 97,66 -50,76 76,68 -71,73 122,98 -25,43 148,61 0,196185
0,60 403,43 64,00 -339,43 129,75 -273,68 244,14 -159,29 244,14 -159,29 182,66 -220,77 321,97 -81,46 404,07 0,640028
Глобальная ошибка -510,63 -404,89 -231,40 -231,40 -323,52 -117,29 0,915
По строке 2 видно, что наименьшую ошибку теперь дают вторая модификация метода золотого сечения Езолото2 = -0,1 ~ О(Н) и оптимальный метод Дшт™ = 0,0007 ~ О(Н4) при значении аоптим « 0,7190001. Графические иллюстрации приведены на рис. 1.
1,2
1,0 0,8 0,6 0,4 0,2 0,0 -0,2
\
\
0,0
0,2
0,4
0,6
Эйлер
Эйлер2
■ Оптим
500 400 300 200 100
0
0,0
0,2
0,4
0,6
Эйлер
Эйлер2
Оптим
а б
Рис. 1. Точность модификаций метода Эйлера: а - для уравнения у' = а1 при а = -10, Н = 0.1, у0 = 1; б - для уравнения у' = а1 при а = +10, Н = 0.1, у0 = 1
Расчеты для меньших значений шага приведены ниже.
Таблица 3
Точность модификаций метода Эйлера для уравнения У = ау при Н = 0.01, а = -10 (у0 = 1)
о X ю i ю X •а ю = ю 1 ю 1 ю ¡X 5 i ю
о н ;= £0 13 О U я ч э о ;= £0 13 о & о 13 О СП 13 О СП 13 О с = О § s 13 О
0,00 1,00 1,00 0,00 1,00 0,00 1,00 0,00000 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,000000
0,01 0,90 0,90 0,00 0,90 0,00 0,91 0,00016 0,95 0,05 0,94 0,03 0,96 0,06 0,90 0,000063
0,10 0,37 0,35 -0,02 0,36 -0,01 0,37 0,00066 0,60 0,23 0,53 0,16 0,68 0,31 0,37 0,000255
0,20 0,14 0,12 -0,01 0,13 -0,01 0,14 0,00049 0,36 0,22 0,28 0,14 0,46 0,32 0,14 0,000187
0,30 0,05 0,04 -0,01 0,05 0,00 0,05 0,00027 0,21 0,16 0,15 0,10 0,31 0,26 0,05 0,000103
0,40 0,02 0,01 0,00 0,02 0,00 0,02 0,00013 0,13 0,11 0,08 0,06 0,21 0,19 0,02 0,000051
0,50 0,01 0,01 0,00 0,01 0,00 0,01 0,00006 0,08 0,07 0,04 0,03 0,14 0,14 0,01 0,000023
0,51 0,01 0,00 0,00 0,01 0,00 0,01 0,00006 0,07 0,07 0,04 0,03 0,14 0,13 0,01 0,000022
0,52 0,01 0,00 0,00 0,00 0,00 0,01 0,00005 0,07 0,06 0,04 0,03 0,13 0,13 0,01 0,000020
0,53 0,00 0,00 0,00 0,00 0,00 0,01 0,00005 0,07 0,06 0,03 0,03 0,13 0,12 0,01 0,000018
0,54 0,00 0,00 0,00 0,00 0,00 0,00 0,00004 0,06 0,06 0,03 0,03 0,12 0,12 0,00 0,000017
Глобальная ошибка -0,01 -0,01 0,0004 0,47 0,23 0,92 0,000153
Оптимальное значение aQm™ ~ 0,049 обеспечивает на втором шаге Дшт™ = 0,00006 ~ O(h3) и глобальную точность Еоптим = 0,0001 ~ O(h2).
Таблица 4
Точность модификаций метода Эйлера для уравнения у' = ау при Н = 0.01, а = +10 (у0 = 1)
о X ю 1 ю •а ю = ю 1 ю ю -а Ц >х ю
о н ;= £0 13 О U я Ч 13 О ;= £0 13 О & О 13 О о СП 13 О о СП 13 О z а X О 13 О
0,00 1,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,00 1,00 0,000000
0,01 1,11 1,10 -0,01 1,10 0,00 1,10 0,00 1,15 0,04 1,14 0,03 1,16 0,06 1,11 0,000001
0,10 2,72 2,59 -0,12 2,65 -0,06 2,65 -0,06 4,05 1,33 3,65 0,93 4,48 1,76 2,72 0,000027
0,20 7,39 6,73 -0,66 7,04 -0,35 7,04 -0,35 16,37 8,98 13,32 5,93 20,07 12,68 7,39 0,000145
0,30 20,09 17,45 -2,64 18,68 -1,41 18,68 -1,41 66,21 46,13 48,59 28,51 89,94 69,85 20,09 0,000590
0,40 54,60 45,26 -9,34 49,56 -5,04 49,56 -5,04 267,86 213,27 177,31 122,72 402,95 348,36 54,60 0,002138
0,50 148,41 117,39 -31,02 131,50 -16,91 131,50 -16,91 1083,66 935,24 647,04 498,62 1805,38 1656,97 148,42 0,007265
0,51 164,02 129,13 -34,89 144,98 -19,04 144,98 -19,04 1246,21 1082,18 736,46 572,44 2097,49 1933,47 164,03 0,008189
0,52 181,27 142,04 -39,23 159,84 -21,43 159,84 -21,43 1433,14 1251,86 838,24 656,96 2436,86 2255,59 181,28 0,009228
0,53 200,34 156,25 -44,09 176,22 -24,11 176,22 -24,11 1648,11 1447,77 954,08 753,74 2831,15 2630,81 200,35 0,010395
0,54 221,41 171,87 -49,53 194,29 -27,12 194,29 -27,12 1895,32 1673,92 1085,93 864,53 3289,23 3067,82 221,42 0,011705
Глобальная ошибка -250,83 -136,96 -136,96 7896,89 4158,42 14180,35 0,058934
Оптимальное значение аоптим « 0,05172 обеспечивает на втором шаге Еоптим = 0,000001 ~ O(h3) и высокую глобальную точность.
Случай 2. Если функция у(0 лежит под касательной (выпукла вверх, y'(tn) Ф 0, y"(tn) < 0), то все происходит с точностью до «наоборот». На каждой итерации метод Эйлера завышает значение yn+1 по сравнению с точным y(tn) на ту же величину, и в алгоритме (9) нужно поменять знак перед an*
Уп+1 « Уп + h(1 - an*)yn', n = 0, 1, 2,..........(9 ')
Для удобства записи соотношения (9) и (9 ' ) можно записать в форме одного общего алгоритма
Уп+1 « Уп + h(1 + sign(yn")an*)yn', n = 0, 1, 2,..........(12)
Заметим, что в этой форме он также использует вторую производную. Этого, однако, тоже можно избежать, если заменить точные признаки выпуклости на менее точные эвристические алгоритмы.
Пример 2. Рассмотрим то же уравнение у' = М с отрицательным начальным условием. Тогда функция, представляющая собой точное решение, выпукла вверх (лежит под касательной), табл. 5-6.
Таблица 5
Точность модификаций метода Эйлера для уравнения у' = ау при Н = 0.1, а = -10 (у0 = -1)
г о X ю и ю •а ю = ю 1 ю 1 ю ¡X £ £ ю
о н ;= £0 13 о и я Ч 13 о ;= £0 13 О & о 13 О о 13 О о 13 О О § 13 О
0,00 1,00 -1,00 -1,00 -1,00 -1,00 -1,00 -1,00 -1,00
0,10 -0,37 0,00 0,37 -0,25 0,12 -0,50 -0,13 -0,50 -0,13 -0,38 -0,01 -0,62 -0,25 -0,37 -0,000021
0,20 -0,14 0,00 0,14 -0,06 0,07 -0,25 -0,11 -0,25 -0,11 -0,15 -0,01 -0,38 -0,25 -0,14 -0,000015
0,30 -0,05 0,00 0,05 -0,02 0,03 -0,13 -0,08 -0,13 -0,08 -0,06 -0,01 -0,24 -0,19 -0,05 -0,000008
0,40 -0,02 0,00 0,02 0,00 0,01 -0,06 -0,04 -0,06 -0,04 -0,02 0,00 -0,15 -0,13 -0,02 -0,000004
0,50 -0,01 0,00 0,01 0,00 0,01 -0,03 -0,02 -0,03 -0,02 -0,01 0,00 -0,09 -0,08 -0,01 -0,000002
0,60 0,00 0,00 0,00 0,00 0,00 -0,02 -0,01 -0,02 -0,01 0,00 0,00 -0,06 -0,05 0,00 -0,000001
Глобальная ошибка 0,58 0,25 -0,40 -0,40 -0,04 -0,95 -0,00005
Как видно, сравнительные характеристики точности методов не изменились. Наилучшие результаты, по-прежнему, показывает оптимальный метод Дшт™ = 0,00002 ~ 0(И5) при значении аоптим « 0,3679, что совпадает с результатами табл. 1, полученными для такого же шага.
Таблица 6
Точность модификаций метода Эйлера для уравнения у' = -ау при Н = 0.1, а = +10 (у0 = 1)
о X я а ю 1 я а ю •а я а ю X я а ю 1 я а ю 2 ¡2 я а ю ¡X £ £ ю
о н ;= £0 Е О и * 13 О ;= £0 13 О а О 13 О 15 о 13 О 15 о 13 О о 1 13 О
0,00 1,00 -1,00 -1,00 -1,00 -1,00 -1,00 -1,00 -1,00
0,10 -2,72 -2,00 0,72 -2,25 0,47 -2,50 0,22 -2,50 0,22 -2,38 0,34 -2,62 0,10 -2,72 -0,000718
0,20 -7,39 -4,00 3,39 -5,06 2,33 -6,25 1,14 -6,25 1,14 -5,67 1,72 -6,85 0,54 -7,39 -0,003905
0,30 -20,09 -8,00 12,09 -11,39 8,69 -15,63 4,46 -15,63 4,46 -13,52 6,57 -17,94 2,14 -20,10 -0,015924
0,40 -54,60 -16,00 38,60 -25,63 28,97 -39,06 15,54 -39,06 15,54 -32,19 22,40 -46,98 7,62 -54,66 -0,057722
0,50 -148,41 -32,00 116,41 -57,67 90,75 -97,66 50,76 -97,66 50,76 -76,68 71,73 -122,98 25,43 -148,61 -0,196158
0,60 -403,43 -64,00 339,43 -129,75 273,68 -244,14 159,29 -244,14 159,29 -182,66 220,77 -321,97 81,46 -404,07 -0,639939
Глобальная ошибка 510,63 404,89 231,40 231,40 323,52 117,29 -0,914366
Как видно, наилучшие результаты по-прежнему показывает оптимальный метод Еоптим = 0,0007 ~ 0(И4) при значении аоптим ~ 0,7190001, что совпадает с данными табл. 2.
Случай 3. Пусть функция у(0 имеет перегиб (лежит по разные стороны от касательной, у'(¿„) Ф 0, у"(^) = 0). Обращение второй производной в нуль -необходимое, но не достаточное условие перегиба [4. С. 341]. Тогда справа и слева от точки перегиба решение лежит над касательной (для четных функций) и под касательной (для нечетных функций). Для прямого времени (эволюция вперед) определяющим является поведение справа (для обратного
времени - слева). Для этого вырожденного, но не типичного случая алгоритм (12) можно скорректировать так:
Уп+1 « Уп + h(1 + sign(y"(yn + h/2))an*)yn ', п = 0, 1, 2,..........(12 ')
Анализ показал, что точность найденных модификаций метода Эйлера (8)-(12) превосходит и высокоточные методы, например Рунге - Кутта, имеющего, как известно, порядок точности ~ O(h5).
Многомерные системы. Все приведенные рассуждения обобщаются на многомерный случай, т.е. на системы ОДУ, аналогично классическим многомерным методам Эйлера. При этом скалярное соотношении (13) следует рассматривать в векторной форме
Уп+1, i«Уп, i + h(1+sign(yn, i")an*)yn, i ', n = 0, 1, 2, ....; i = 1, ..., k. (13)
Рассмотрим примеры применения модификации метода Эйлера (13) для многомерных систем.
Пример 3. Модель эволюции 1) X3 = X1, 2) X1 = X2, 3) X2 + 2X3 = 3X3 описывается системой двух дифференциальных уравнений (простейший осциллятор [6])
X ' 1 = Ю1(1- Х1 - Х2) - Ю-1Х1 - Ю2Х1 + Ю-2Х2, (14)
X'2 = Ю2Х1 - Ю -2X2 - «3X2(1 - Х1 - Х2)2. ( )
Анализ показал, что динамика этой системы может быть монотонной, колебательной (затухающей и незатухающей), т.е. сочетает участки выпуклости вверх (лежат под касательной) и вниз (лежат над касательной). Равновесных состояний может быть три или одно. Если равновесие единственное и неустойчивое, то в системе возникают незатухающие колебания (автоколебания). Например, при ю1 = 2.89, ю-1 = 0.01, ю2 = 3/89, ю-2 = 0.1, ю3 = 2000 возникают автоколебания вокруг неустойчивого фокуса. Результаты численного интегрирования системы (14) для различных модификаций метода Эйлера (Эйлер1 - первого порядка, Эйлер2 - второго порядка, Эйлер3 - третьего порядка, Точно -максимально возможная точность) и нашего метода (Оптимальный) при начальных значениях Х10 = 0.7, Х20 = 0.2 приведены в табл. 7. Максимальная точность определена по порядку ненулевых высших производных
Х ' \ = -ю1 - Ю-1 - Ю2,
Х"2 = -Ю-2 - Ю3(1 - Х1 - Х2)2 + 2«3Х2(1 - Х1 - Х2),
Х ' '' 1 = 0,
х" ' 2 = 4ю3(1 - Х1 - Х2) - 2ю3х2
Х"" 2 = -6Ю3,
х ' ""2 = 0.
Как видно, наилучшие результаты показывает оптимизированный метод, среднеквадратическая ошибка (корень из суммы квадратов, деленный на число точек) минимальна Еоптим = 0,0002 даже при не полностью оптимальном значении a1 = a2 = 0,5, которое соответствует эвристике «среднего», описываемой формулой (11). Графическая иллюстрация приведена на рис. 2.
Таким образом, простой численный метод решения задачи Коши позволяет вычислять решения систем обыкновенных дифференциальных уравнений с высокой точностью и минимальным эффектом «сползания». Этот метод может быть применен для уточнения решений обыкновенных дифференциальных уравнений без использования производных высших порядков.
Таблица 7
Точность модификаций метода Эйлера для системы (14) при Н = 0.01
г Эйлер1 Эйлер1 я а ю 5 •а Э •а л Э я а ю 5 ЭйлерЗ Эйлер3 я а ю 5 о X г о X г я а ю X й X А й X -а я а ю X
О О о н н О |5 X о 15 X О О
0,00 0,70 0,20 0,70 0,20 0,70 0,20 0,70 0,20 0,70 0,20
0,01 0,70 0,16 0,0030 0,70 0,16 0,0000 0,70 0,16 0,0000 0,70 0,16 0,00 0,70 0,17 0,00525307
0,02 0,71 0,10 0,0070 0,71 0,11 0,0001 0,71 0,11 0,0000 0,71 0,11 0,00 0,70 0,11 0,00491431
0,03 0,71 0,03 0,0079 0,71 0,03 0,0002 0,71 0,03 0,0000 0,71 0,03 0,00 0,71 0,03 0,00700015
0,04 0,72 -0,01 0,0052 0,72 -0,01 0,0003 0,72 -0,01 0,0000 0,72 -0,01 0,00 0,71 -0,02 0,00472926
0,05 0,73 0,01 0,0048 0,73 0,00 0,0001 0,73 0,00 0,0000 0,73 0,00 0,00 0,72 0,00 0,00467941
0,10 0,76 0,00 0,0050 0,76 -0,01 0,0003 0,76 0,00 0,0000 0,76 0,00 0,00 0,76 0,00 0,00359888
0,20 0,82 0,00 0,0050 0,82 0,00 0,0004 0,82 0,00 0,0000 0,82 0,00 0,00 0,82 0,00 0,00266897
0,30 0,86 0,00 0,0049 0,86 0,00 0,0005 0,86 0,00 0,0000 0,86 0,00 0,00 0,86 0,00 0,00197777
0,40 0,89 0,00 0,0049 0,89 0,00 0,0006 0,89 0,00 0,0000 0,89 0,00 0,00 0,89 0,00 0,00146269
0,50 0,92 0,00 0,0049 0,92 0,00 0,0008 0,92 0,00 0,0000 0,92 0,00 0,00 0,92 0,00 0,00107800
1,00 0,97 0,01 0,0047 0,97 0,00 0,0019 0,97 0,00 0,0002 0,97 0,00 0,00 0,97 0,00 0,00021155
25,00 0,91 0,09 0,0302 0,67 0,32 0,3075 0,88 0,11 0,0132 0,88 0,10 0,00 0,88 0,10 0,00012234
25,10 0,91 0,09 0,0257 0,67 0,32 0,3104 0,88 0,11 0,0150 0,88 0,10 0,00 0,88 0,10 0,00015355
25,20 0,91 0,09 0,0212 0,67 0,32 0,3140 0,88 0,11 0,0173 0,89 0,09 0,00 0,89 0,09 0,00019311
25,30 0,90 0,09 0,0172 0,67 0,32 0,3184 0,88 0,11 0,0200 0,89 0,09 0,00 0,89 0,09 0,00024412
25,40 0,90 0,09 0,0151 0,67 0,33 0,3238 0,88 0,11 0,0235 0,89 0,09 0,00 0,89 0,09 0,00031081
25,49 0,90 0,09 0,0171 0,67 0,33 0,3298 0,88 0,10 0,0274 0,89 0,08 0,00 0,89 0,08 0,00038827
Ошибка 0,0563 0,1907 0,0063 0,0000 0,0002
......... Эйлер1 — • — Эйлер2------ЭйлерЗ -Оптим
Рис.2. Зависимость х1(Г) для системы (14) при ю1 = 2.89, ю_1 = 0.01, ю2 = 3/89, ю_2 = 0.1, ю3 = 2000, х10 = 0.7 и х20 = 0.2
Представленный в работе новый метод позволяет уменьшить погрешность численного исследования моделей динамических процессов в сложных технических системах.
Литература
1. БерезинИ.С., ЖидковН.П. Методы вычислений: в 2 т. М.: Физматлит, 1962. Т. 1. С. 307.
2. Дьяконов В., Круглов В. Математические пакеты расширения Ма1ЬаЪ. СПб.: Питер,
2001.
3. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1974. 832 с.
4. Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления: в 3 т. М.; Л.: Гостехтеоретиздат, 1949. Т. 1. 690 с.
5. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980, 307 с.
6. Fedotov V.Kh., Alekseev B.V., Koltsov N.I. Self-oscillations in three step catalytic reactions. React. Kinet. Catal. Lett., 1983, vol. 23, no. 3-4, pp. 301-306.
ФЕДОТОВ ВЛАДИСЛАВ ХАРИТОНОВИЧ - кандидат химических наук, доцент кафедры информационных систем, Чувашский государственный университет, Россия, Чебоксары ([email protected]).
НОВОЖИЛОВА НИНА ВАСИЛЬЕВНА - кандидат экономических наук, доцент кафедры информационных систем, Чувашский государственный университет, Россия, Чебоксары ([email protected]).
V. FEDOTOV, N. NOVOZHILOVA
HIGH-ACCURACY METHOD FOR NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS
Key words: mathematical modeling, numerical integration of differential equations, highly accurate numerical method for solving the Cauchy problem, dynamic processes in the fields of science and engineering, numerical methods, computer calculations of a complex technical system.
The accuracy of the numerical integration of differential equations depends on a little-studied laws of algebra of approximate computation - errors of the methods, sampling, stability of difference schemes, etc. As a consequence, at each iteration a problem, different from the previous is solved, which leads to accumulation of calculation errors (the effect of «slipping»). In the work a simple high-precision numerical method for solving the Cauchy problem without using higher derivatives and with minimum effect of "slipping" is described and tested.
References
1. Berezin I.S., Zhidkov N.P. Metody vychislenii: v 2 t. [Methods of calculations. 2 vols]. Moscow, Fizmatlit Publ., 1962, vol. 1, p. 307.
2. D'yakonov V., Kruglov V. Matematicheskie pakety rasshireniya MatLab [Mathematical expansion packs MATLAB]. St. Petersburg, Piter Publ., 2001.
3. Korn G., Korn T. Spravochnikpo matematike dlya nauchnykh rabotnikov i inzhenerov [Mathematical Handbook for scientists and engineers]. Moscow, Nauka Publ., 1974.
4. Fikhtengol'ts G.M. Kurs differentsial'nogo i integral'nogo ischisleniya: v 3 t. [A course of differential and integral calculus. 3 vols]. Moscow, Leningrad, Gostekhteoretizdat Publ., 1949, vol. 1, 690 p.
5. Forsait Dzh., Mal'kol'm M., Mouler K. Mashinnye metody matematicheskikh vychislenii [Machine methods of mathematical calculations]. Moscow, Mir Publ., 1980.
6. Fedotov V.Kh., Alekseev B.V., Koltsov N.I. Self-oscillations in three step catalytic reactions. React. Kinet. Catal. Lett., 1983, vol. 23, no. 3-4, pp. 301-306.
FEDOTOV VLADISLAV - Candidate of Chemical Sciences, Associate Professor of Information Systems Department, Chuvash State University, Russia, Cheboksary ([email protected]).
NOVOZHILOVA NINA - Candidate of Economics Sciences, Associate Professor of Information Systems Department, Chuvash State University, Russia, Cheboksary ([email protected]).
Ссылка на статью: Федотов В.Х., Новожилова Н.В. Высокоточный метод численного решения дифференциальных уравнений // Вестник Чувашского университета. - 2016. - № 1. -С. 174-184.
