Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 10
Текст из файла (страница 10)
Разлагая сначала по переменной х, получаем ф =С2~+ Сз1ЯХ, У +С1 ЙД+ С1Ц~(х, 7 + С1ЬГ1+ = (сз +сз)Г+ сз сз 1з Пу +Ух) + 2 + с1 сзп (Х Уу»+20ху Ухх)+0(л )- 2 (2.2.25) 1 аУ 1 , ~'1 з.+ й + 1,з +0(йз) ~+ й(~~. +~ )+ 2 Ых 6 Ихз 2 йз + (У .Гу» + 2Их» +Ухх +Ух Уу +Пу ) + 0(Ь. ) 6 (2.2.26) Объединяя теперь (2.2.25) и (2.2.26), получаем 1 — [у(х+Ь) — у(х)1 — р(х„у(х)) = й /1 (1 сз сз) 1+ 1~ с1сз (Оу +Хх) + (2.2.27) 2 2 йз + — ~ — — Сз~ сз) (3'~уу + 2Д~ху + ~хх) + (~х ~у + Д~", ) + 0(йз) Если мы потребуем, чтобы сз +сз =1, с,сз =1/2, (2.2.28) то первые два члена в. (22,27) обратятся в нуль для любой функции ~, Однако член ~хХ» + Цу в общем случае не будет тождественно равен нулю; так что, как бы мы ни выбирали константы, самое большое, что мы можем иметь 1 — [ У(х + Ь) — у(х)~ — ~о (х, У(х)) = 0(йз), Ь (2.2.29) что и будет выполняться всякий раз, когда удовлетворены условия (2.2.28) и все необходимые производные ограничены.
Таким образом, Х (Ь) = = 0(йз) и все методы, выделенные условием (2.2.28), имеют второй порядок. Если положить с, = 7/2 и выразить сз и сз из двух уравнений (2.2.28) через 7, то мы получим функцию у, всегда удовлетворяющую (2.2.28). Таким образом, метод 1~ 1 / 7 7л У, =У, + Ь 1 — — К(хр„У~) + ~~ха + л Уа + Лх(с.Уа) 7 7 имеет второй порядок точности для любого значения 7Ф О. Специальный выбор 7 = 2 дает метод Рунге — Кутта второго порядка (2„2.19) . Построе- С другой стороны, для точного решения дифференциального уравнения у(х) имеем 1 1 ,(х+ ь) У(х)[ -У~(х)+ У~(х)ь+ У~ (х)йз + 0(йз) л 2 6 ние методов Рунге — Кутта более высокого порядка и, в частности, метода четвертого порядка (2.2.20) может быть проведено совершенно аналогично, но требует более громоздких выкладок. Теперь покажем, как описанные выше методы, можно использовать для систем уравнений.
Рассмотрим систему (2.2.1), которую запишем в векторной форме: у (х) =Д(х, у(х)). (2.2.30) Здесь у(х) — вектор с координатами у, (х),...,у„(х), а à — вектор с л координатами ~~,...,~„. Через у обозначим вектор, составленный из начальных значений (2.2.2) . Тогда метод Эйлера (2.2,5) для системы (2.2.30) может быть записан как л уе =у, ук+1 =ук+ЙД(хк,ук), й=0,1,..., (2.2.31) где ус,уг,... — векторы приближений к решению у в соответствующих точках. Мы могли бы, конечно, выписать соотношения (2.2.31) покоординатно; при п = 2 это бы выглядело так: л л У1о =У1 Уао =Уз. У1,к+1 У1к + М1 (хк У~ к У~ к ) 1 й=0,1, Ут,к+1 =Угк +Из(хк У1к.У1к) ~ Очевидно, что краткая векторная запись (2.2.31) гораздо более удобна. Аналогично, метод Хьюна (2.2.19) для системы (2.2.30) может быть записан в векторной форме следующим образом: 6 Ук+1=ук+ — У(хк Ук)+йхк+1 Ук+Ю(хк,ук))1.
(2.2.32) В упражнении 2.2.7 предлагается записать в векторной форме метод Рунге— Кутта четвертого порядка (2.2,20) . Теперь кратко обсудим влияние ошибок округления, возникающих при реализации на ЭВМ методов этого раздела. Рассмотрим сначала метод Эйлера. Здесь имеется два источника ошибок окрутления. Во-первых, это ошибка, возникающая при вычислении ~(хк,ук). Обозначим этуошибку через ек. Во-вторых, ошибка пк, возникающая при счете по самой формуле Эйлера. Таким образом, фактически вычисленные приближения ук удовлетворяют соотношениям Ук+1=ук+й ~Яхк,ук)+ек]+ ск 1=0,1, (2.2.33) Можно получить оценку влияния этих ошибок через оценки величин е» и лк.
Мы, однако, ограничимся следующим интуитивным рассуждением. Как мы видели, глобальная ошибка дискретизации метода Эйлера стремится к нулю при стремлении Ь к нулю. Следовательно, за счет уменьшения шага й мы можем сделать ошибку дискретизации сколь угодно малой. Однако чем меньше сс, тем больше потребуется шагов по методу Эйлера и, вообще говоря, тем больше скажутся на полученном решении ошибки округления. На практике, когда прн выполнении арифметических операций используются слова фиксированной длины, всегда существует такая величина шага й, меньше которой вклад ошибок округления начинает доминировать в суммарной ошибке. Эта ситуация схематически изображена иа рис.
2.6, где йе определяет минимальный шаг, который можно использовать в практических расчетах. Эту минимальную величину шага очень трудно установить зарнее, но в задачах, где не требуется слишком высокая точность, необходимый шаг обычно будет значительно больше, чем этот минимум, и основной вклад в полную ошибку будет вносить ошибка дискретизации.
Такое поведение характерно и для других методов, хотя минимальный размер шага будет меняться от метода к методу и от задачи к задаче. 4$ Ю и Ю к Ю Рис. 2.б. Поведение ошибки в методе Эйлере и других методах В заключение этого раздела приведем результаты некоторых простых расчетов для уравнений модели хищник — жертва, полученных в разделе 2.1. Напомним, что это были уравнения ду — = уу+бху, г> О, Н Ых — = ах + 13ху, с~г (2.2.34) с начальными условиями х(О) =х, у(О) =у . (2.2.35) В дополнительных замечаниях к разделу 2,1 было показано, что точка х,= — 7~б, у = — а/Р является стационарной для системы (2.2.34) и что в окрестности точки (х„у,) траектории решений (х(г), у(г)) при г > 0 близки к эллипсам.
В иллюстративных целях мы выбрали начальные значения хе и уа вблизи стационарной точки. Мы использовали также следующие. значения пара. метров: а = 0,25, Р = — 0,01, 7 = — 1,00 и б = 0,01. При этих параметрах стационарной является точка (х„у,) = (100,25). Во всех случаях использовались начальные значения хе = 80 и уе = 30. На рис. 2.7 — 2.9 представлены приближенные решения задачи (2.2.34)— (2.2.35), полученные в результате использования некоторых из изложенных в этом разделе численных методов.
Во всех слу иях мы представили х (численность жертв) и у (численность хищников) как функции времени г в виде параметрических кривых на плоскости (х, у) . При возрастании г движение происходит по часовой стрелке. На рис. 2.7 показаны три возможности, которые обычно доступны пользователям графических систем. На рис. 2.7,а нанесены только дискретные значения (хр,у~) (~ = = О, 1,... ), полученные численным методом, что подчеркивает дискретную природу этих методов. На рис. 2.7,б эти точки соединены отрезками прямых, придающими приближению многоутольную форму. Это обычно саьаая простая возможность„реализуемая графопостроителем. На рис.
2.7,в точки соединены гладкими кривыми. Эта возможность требует специального программного обеспечения, основанного на методах аппроксимации, подобных тем, которые будут рассмотрены в разделе 5.2, На рис. 2.8 показана зависимость приближенного решения от величины шага Л, Представленные на этом рисунке результаты получены по методу Эйлера (2.2.5). При этом использовались шаги 1, 0,5 и 0,25.
Видно, что при уменьшении шага в два раза ошибка также уменьшается примерно в два раза, что соответствует сходимости со скоростью 0(л). Яоно, что ошибка весьма велика даже при )1 = 0,25. Использованное для 12 12 10 80 ВО ) ) ) 60 О 20 40 0 20 40 О 20 4О а б б Рис. 2.7. Возможности графического вывода, Решение задачи (2.2.34) — (2.2.35): а) только данные; б) данные, соединенные отрезками; в) данные, интерполированные функцией х(г 160 12 10 О 20 40 ВО ф~Ц О 20 40 У(Г) Рис. 2.8. Метод Эйлера для задачи (2.2.34) — (2.2.35). Использованы три различные величины шага Рис.
2.9. Метод Рунге — Кутта второго порядка (раскручивающаяся спираль) при шаге Ь = 1 в сравнении с "точным" решением задачи (2.2.34) — (2.2. 35) 42 сравнения "трчное" решение получено по методу Рунге-Кутта высокого порядка. В целях сравнения с методами более низкого порядка так полученное решение вполне можно рассматривать как точное. На рис. 2.9 продемонстрирован эффект от использования метода второго порядка, а не метода Эйлера, имеющего первый порядок. Здесь ошибка при шаге л = 1 меньше, чем ошибка метода Эйлера при шаге й = 0,25. Обратите внимание, по, как и в случае метода Эйлера, приближенное решение представляет собой раскручивающуюся спираль, удаляющуюся от точного решения, Здесь, вероятно, уместно сделать несколько поясняющих замечаний.
Конечно, задача (2.2.34) — (2.2.35) не моделирует абсолютно точно физическую ситуацию. Мы не учли многие факторы, такие,как ограниченность пищи для жертв, вмешательство других хищников, впияние погоды и т.д. Такие упрощения были, в частности, необходимы и для возможности дать описание модели в этой главе. Действительной проверкой модели является сравнение полученных по модели результатов с конкретной физической ситуацией. Вычисленные решения задачи (2.2,34) — (2,2.35) сравнивались с фактическими данными о популяциях видов, При этом оказалось, что форма кривой, представленной на рис. 2.7, в, весьма правдоподобно аппроксимирует реальные данные, хотя, как и следовало ожидать, результаты расчетов ведут себя значительно более гладко. Как и в случае реальной ситуации, хищники не полностью истребляют свои жертвы, При уменьшении численности жертв популяция хищников также уменьшается, поэтому шансы жертв на выживание возрастают.
Поведение популяции хищников все время следует за изменением популяции жертв, и цикл повторяется. Дополнительные замечания и ссылки 2.2 Вероятно, наиболее простой с концептуальной точки зрения ладход к построению одношаговых методов высокого порядка заключается в разложении решения в ряд Тейлора. Рассмотрим следующий '*метод": ьР ух+1 =уа+Ьу'(ху,)+ — Ь'у" (ха)+... + — у~р)(хя), (2.2.36) гдеу — точное решение.
Легко видеть, что такой метод имеет порядок р. Старшие производные решения можно, в принципе, определить из самого дифференциального уравнения. Действительна, у'(х) =)'(х, у(х)), с у" (х) = — 1(х, у(х)) =~х(х, у(х))+~ (х, у(х))у'(х). (2.2.37) Затем мы заменяем у' (ха) на у (хь, уь) и аналогично поступаем со старшими производными. Таким образом, этот метод при р = 1 является просто методом Эйлера, а при р = 2 принимает вид й' ур,~1 = уя + Ь~(хх, у~,) + — ( ~х (х~„у~,) + ~у (х~„уа)~(х~„у>,)] 2 и является методом второго порядка. Можно продолжить дифференцирование (2.2.37) и выразить более высокие производные решения у через старшие частные производные 7 (х,у), но получаемые при этом методы становятся чрезвычайно громоздкими. При определении выражений для этих производных в известной мере полезной оказалась техника символьных вычислений, Дэлънейшее развитие методов на основе разложения в ряд Тейлора можно найти в книге 112] .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
