Реферат
Решение многихтехнических, химических, а также биологических задач требует решения задачиКоши. Эту задачу можно решать разными способами, как аналитическими, так ичисленными, применяя ЭВМ. Очень часто бывает важно получить результат в сжатыесроки. В этом случае предпочтение отдается численным методам. Кроме того,встречаются такие сложные дифференциальные уравнения, найти аналитическоерешение которых либо вообще не представляется возможным, либо для этоготребуются очень большие затраты времени и сил.
В работе детально рассматриваетсяметод Рунге-Кутты четвертого порядка с автоматическим выбором длины шагаинтегрирования (это обеспечивает гораздо более высокую точность вычислений посравнению с методом, использующим шаг постоянной длины), приводится необходимаятеоретическая сводка, описание метода, а также программа для ЭВМ, результаты еевыполнения и иллюстрации.
Ключевые слова: дифференциальное уравнение, методРунге-Кутты, метод Эйлера, порядок метода Рунге-Кутты, задача Коши, ряд Тейлора,отрезок, коэффициенты, шаг интегрирования, интегральная кривая.
Работа содержит 36 листов, включая 8 графиков, 4иллюстрации и 12 таблиц.
Содержание
Введение
1. Теоретическая часть
1.1 Постановка задачи
1.2 Метод Эйлера
1.3 Общая формулировка методовРунге-Кутты
1.4 Обсуждение методов порядка 4
1.5 «Оптимальные» формулы
1.6 Условия порядков для методовРунге-Кутты
1.7 Оценка погрешности и сходимостьметодов Рунге-Кутты
1.7.1 Строгие оценки погрешности
1.7.2 Главный член погрешности
1.7.3 Оценка глобальной погрешности
1.8 Оптимальный выбор шага
2. Практическая часть
2.1 Описание программы «Ilya RK-4 версия 1.43» Заключение Список использованных источников Приложение А. Графики функций Приложение Б. Пример таблицы значений функции y(x) Приложение В. Листинг программы «Ilya RK-4 версия 1.43»
Введение
Ввиду того, что дляметодов Рунге-Кутты не нужно вычислять дополнительные начальные значения, этиметоды занимают особое место среди методов классического типа. Ниже будутрассмотрены их свойства, а также некоторые ограничения, присущие этим методам.
С увеличением числа этапов длябольших задач, решаемых этими методами, возникли бы трудности с памятью ЭВМ,кроме того (и это важнее), для больших задач, как правило, всегда великиконстанты Липшица. В общем случае это делает методы Рунге-Кутты высокогопорядка не пригодными для таких задач. Во всяком случае, другие методы обычноэффективнее и им следует отдавать предпочтение. Однако методы Рунге-Куттычетвертого порядка являются достаточно легко реализуемыми на ЭВМ, а наличиеавтоматического выбора шага дает возможность производить вычисления с хорошейточностью. Поэтому их целесообразно применять для довольно широкого множествазадач.
Методы Рунге-Кутты имеютнесколько весомых достоинств, определивших их популярность среди значительногочисла исследователей. Эти методы легко программируются, обладают достаточнымидля широкого круга задач свойствами точности и устойчивости. Эти методы, как ивсе одношаговые методы, являются самостартующими и позволяют на любом этапевычислений легко изменять шаг интегрирования.
В работе основноевнимание сконцентрировано на вопросах точности и эффективности решения задачтого типа, для которых методы Рунге-Кутты приемлемы.
Программная реализацияметодов Рунге-Кутты четвертого порядка с автоматическим выбором шагапредставлена в виде программы, написанной на языке высокого уровня BorlandC++ 3.1. Программу можно запускать в среде MS-DOS или Windows® 95/98/Me/2k/XP. В качестве выхода программа пишеттаблицу значений в файл на диск и рисует график на экране ЭВМ.
Для проверки результатовработы созданной программы одни и те же дифференциальные уравнения решались вматематическом пакете WaterlooMaple9.01 и припомощи созданного приложения (версия 1.43), проводился анализ таблиц значений играфиков решений.
1. Теоретическая часть1.1 Постановка задачи
Дано дифференциальное уравнение и начальное условие, то есть поставлена задача Коши:
/> (2.1.1)
Требуется отыскать интегральную кривую,удовлетворяющую поставленной задаче Коши с помощью метода Рунге-Куттычетвертого порядка с автоматическим выбором шага на отрезке />. Задачу можно решить аналитически, найдярешение дифференциального уравнения и подставив в него начальное условие, темсамым, отыскав требуемую интегральную кривую. Но для нас интерес представляетрешение данной задачи с применением численного метода, а конкретнее – методаРунге-Кутты 4-го порядка с автоматическим выбором шага, то есть численноерешение. Автоматический выбор шага – необходимое условие адекватного поведенияпрограммы при резко изменяющихся функциях, задающих интегральную кривую,позволяющее отразить все моменты в поведении интегральной кривой и добитьсявысокой точности.
1.2 Метод Эйлера
Метод Эйлера для решения начальной задачи (2.1.1) былописан Эйлером в 1768 году. Этот метод весьма прост. Его глобальная погрешностьимеет вид />, где /> –постоянная, зависящая от задачи, и /> – максимальная длина шага.Если желательно, скажем, получить 6 точных десятичных знаков, то требуется,следовательно, порядка миллиона шагов, что не слишком удовлетворительно. Сдругой стороны, еще со времен Ньютона известно, что можно найти гораздо болееточные методы, если /> не зависит от />, то есть если мы имеем задачу (2.1.1),решаемую квадратурой
/>. (2.2.1)
В качестве примера можно рассмотреть первуюквадратурную формулу Гаусса, также называемую «правилом средней точки»:
/> (2.2.2)
где /> и /> –граничные точки подинтервалов, на которые разбит интервал интегрирования.Известно, что оценка глобальной погрешности этой формулы /> имеет вид />. Такимобразом, если желаемая точность составляет 6 десятичных знаков, ее обычно можнополучить приблизительно за 1000 шагов, то есть этот метод в тысячу раз быстрее.Поэтому Рунге поставил следующий вопрос: нельзя ли распространить этот метод наисходную задачу Коши? Первый шаг длины /> должениметь вид
/>. (2.2.3)
Но какое значение взять для />?За неимение лучшего естественно использовать один малый шаг метода Эйлера длины/>. Тогда из предыдущей формулы получим:
/> (2.2.4)
Решающим обстоятельством здесь является умножение /> в третьем выражении на />, в результате чего влияние погрешностистановится менее существенным. Точнее, вычислим для /> разложениеТейлора по степеням />:
/> (2.2.5)
Его можно сравнить с рядом Тейлора для точногорешения, который получается из того, что /> путемповторного дифференцирования с заменой /> на /> каждый раз, когда оно появляется:
/> (2.2.6)
Вычитая из последнего равенства предыдущее, получимдля погрешности первого шага выражение
/> (2.2.7)
Таким образом, если все частные производные /> второго порядка ограничены, то
/>.
Чтобы получить приближенное значение решения исходнойзадачи в конечной точке />, будем применять формулы(2.2.4) последовательно к интервалам />. Приведенные вышеформулы являются усовершенствованным методом Эйлера. Для вычислений с высокойточностью, однако, следует пользоваться другими методами, одним из которых какраз является метод Рунге-Кутты.
1.3 Общая формулировка методов Рунге-Кутты
Рунге и Хойн построили новые методы, включив вуказанные формулы один или два добавочных шага по Эйлеру. Но именно Куттасформулировал общую схему того, что теперь называется методом Рунге-Кутты.
Пусть /> – целое положительноечисло (число стадий, этапов) и /> – вещественныекоэффициенты. Тогда метод
/> (2.3.1)
называется />-стадийным явным методомРунге-Кутты для исходной задачи Коши (2.1.1)
Обычно коэффициенты /> удовлетворяютусловиям
/>. (2.3.2)
Эти условия были приняты Куттом без каких-либокомментариев. Смысл их заключается в том, что все точки, в которых вычисляется />, являются приближениями первого порядка крешению. Эти условия сильно упрощают вывод условий, определяющих порядокаппроксимации для методов высокого порядка. Однако для методов низких порядковэти предположения необходимыми не являются.
Метод Рунге-Кутты имеет порядок />, если для достаточно гладких задач (2.1.1)справедливо неравенство
/>, (2.3.3)
то есть ряды Тейлора для точного решения /> и для /> совпадаютдо члена /> включительно.
После статьи Бутчера вошло в обычай символическипредставлять метод (2.3.1) по средствам следующей таблицы:
/>
/>
/>
1.4 Обсуждение методов порядка 4
Подойдем теперь вплотную к определению 4-стадийныхметодов Рунге-Кутты (2.3.1) с таким расчетом, чтобы они имели порядок 4. Дляэтого необходимо вычислить производные порядков 1, 2, 3 и 4 от /> при /> и сравнитьих с производными точного решения. Теоретически при известных правилахдифференциального исчисления это совершенно тривиальная задача. Однако сиспользованием (2.3.2) получаются следующие условия:
/>
Эти вычисления очень утомительны и емки. Ихгромоздкость очень быстро растет для более высоких порядков.
Лемма 1.
Если
/> (2.4.2)
то уравнения d), g) и h)являются следствием остальных.
Доказательство.
Покажем это для g). Cпомощью уравнений c) и e) получим:
/>
Для уравнений d) и h)процедура аналогична.
Покажем, что в нашем случае условие
/>
является и необходимым.
Лемма 2.
При /> (2.4.2) следует изуравнений (2.4.1) и уравнений (2.3.2).
Для доказательства потребуется следующая лемма 3.
Лемма 3.
Пусть /> и /> суть3x3-матрицы, такие что
/>, (2.4.3)
тогда либо />, либо />, где />.
Доказательство.
Если />, то из /> следует />. Если же />, то существует вектор />, такой, что />, ипоэтому />. Но тогда из (2.4.3) следует, что /> должен быть пропорционален вектору />.
Докажем теперь предыдущую лемму. Введем величины /> для />. Итак,надо доказать, что />. Введем теперь матрицы
/> (2.4.4)
Перемножение этих матриц с использованием условий(2.4.1) дает
/>, (2.4.5)
причем
/>
Далее последний столбец /> не можетбыть нулевым, так как из того, что />, следует
/>
в силу условия h). Такимобразом, из последней леммы следует, что />. Последнеетождество /> вытекает из равенства />, которое является следствием условий a) и b).
Теорема.
Если выполнены предположения />,то уравнения (2.4.1) эквивалентны следующим:
/> (2.4.6)
Доказательство.
Из j) и h) следует, что
/>. (2.4.7)
Отсюда, в частности, вытекает, что в силу k) />.
Решение уравнений (2.4.6). Уравнения a)-e) и k)выражают тот факт, что коэффициенты /> и /> являютсявесами и узлами квадратурной формулы четвертого порядка при /> и />. В силу (2.4.7) возможныследующие четыре случая:
1) />. (2.4.8)
Тогда уравнения a)-e)образуют невырожденную линейную систему для определения />.Эта система имеет решение:
/>
Остальные три случая с двойными узлами основаны направиле Симпсона:
2) />;
3) />;
4) />.
После того, как выбраны /> и />, получаем /> изуравнения j), и тогда два уравнения f) и i)образуют линейную систему для определения /> и />.
Определитель этой системы
/>,
согласно (2.4.7) не равен нулю. Наконец, из того, что /> находим />, /> и />.
Особенно популярными стали два варианта, которыевыбрал Кутта в 1901 году. Это случай 3) при /> и случай1) при />. Оба метода обобщают классическиеквадратурные формулы, сохраняя их порядок. Первый из них более популярен,однако второй более точен.
Правило 3/8
/>
/>
/>
Классический метод Рунге-Кутты
/>
/>
/>
1.5 «Оптимальные» формулы
Предпринималось много исследования, чтобы выбратьвозможно «лучшие» из множества различных формул Рунге-Кутты 4-го порядка.
Первой попыткой в этом направлении был оченьпопулярный метод, который в 1951 году предложил Гилл. Он преследовал цельуменьшить на сколько возможно количество требуемой машинной памяти(«регистров»). Этот метод широко использовался на первых компьютерах впятидесятых годах и представляет поэтому исторический интерес. Гилл заметил,что больше всего машинной памяти нужно при вычислении />,когда «требуются регистры для хранения в какой-либо форме» величин
/>.
Ясно, что для третьей стадии будет достаточно трехрегистров, если подлежащие хранению величины линейно зависимы, то есть если
/>.
Гилл заметил, что это условие удовлетворяется дляметодов типа 3), если />. Получающийся метод можно в такомслучае переформулировать следующим образом:
/>
В настоящее время производительность компьютеров оченьсильно возросла по сравнению с машинами 50-х годов, что привело к прекращениюиспользования данного метода на практике при расчетах. Справедливости радистоит отметить, что этот метод все-таки рационально употреблять в случаях,когда требуется решать системы дифференциальных уравнений очень высокой размерностисо сложными функциями.
1.6 Условия порядков для методов Рунге-Кутты
Рассмотрим структуру условий, определяющих порядокметода, или условий порядка, как их называют для краткости. Способ выводаусловий порядка прошел большую эволюцию. Он совершенствовался главным образомпод влиянием работ Бутчера.
Так как явные методы Рунге-Кутты являются частнымслучаем неявных, то можем выписать условия, при которых метод имеет заданныйпорядок.
Метод
/>
/>
/>
(где на свободных местах должны стоять нули) имеетпорядок />, если удовлетворяется уравнение
/> (2.6.1)
для каждого дерева /> с корнем ине более чем с /> разветвлениями[1].
При /> эти условия,обеспечивающие порядок 4, и соответствующие деревья имеют следующий вид:
/>
/> (2.6.2)
/>
/> (2.6.3)
/>
/> (2.6.4)
/>
/> (2.6.5)
/>
/> (2.6.6)
/>
/> (2.6.7)
/>
/> (2.6.8)
/>
/> (2.6.9)
Заметим, что для меньших значений /> мы берем соответствующее подмножество этихусловий, а для меньших /> оставляем лишь некоторые из указанныхчленов.
Из (2.9) видим, что действительно необходимо 4 этапа,так как если бы их было меньше, то был бы опущен единственный член в левойчасти этого уравнения. Для явных методов в общем случае выполняется неравенство/>. Фактически (для тех значений, для которыхэто известно) минимальное значение /> для данного /> указано в следующей таблице:
/> 1 2 3 4 5 6 7
/> 1 2 3 4 6 7 9
Общие классы методов с этими значениями /> и /> легко найти в случае />.
Для />: 1
Это известный метод Эйлера.
Для />:
/>
/>
/>
/>
Это однопараметрическое семейство имеет требуемыйпорядок для любого ненулевого значения />.
Для /> имеется три семейства, изкоторых первые два таковы:
/>
/> />
/>
/>
/> />
/>
Каждое из них имеет один параметр />. Третье семейство имеет в качестве параметров/> и />, причем
/>.
Вывод методов с /> болеесложен, но его можно упростить, положив
/> (2.6.10)
(что влечет равенство />), так какэто позволяет опустить уравнения (2.6.3), (2.6.5), (2.6.8) и (2.6.9). Интереснотакже, что (2.6.10) является следствием (2.6.2) – (2.6.9).
План вывода конкретного метода этого порядка можновыполнить при условии, что не возникает несовместных систем.
Шаг 1.Выбираем значения />, /> иполагаем />.
Шаг 2. Из(2.6.2), (2.6.3), (2.6.4) и (2.6.6) находим />.
Шаг 3. Изуравнения /> (это уравнение есть разность уравнений(2.6.5) и (2.6.7)) находим />.
Шаг 4. Из(2.6.10) находим />.
Шаг 5.Вычисляем />.
В случае /> шаг 2 приводит к выбору /> и /> при условии, что />, />. В частности, имеемизвестный метод:
/>
/>
/>
1.7 Оценка погрешности и сходимость методовРунге-Кутты
Со времен работы Лагранжа и особенно Коши всякийустановленный численно результат принято сопровождать надежной оценкойпогрешности. Лагранж дал известные оценки погрешности многочленов Тейлора, аКоши вывел оценки для погрешности метода ломаных Эйлера. Через несколько летпосле первых успехов методов Рунге-Кутты также пришел к заключению, что дляэтих методов нужны оценки погрешностей[2].
1.7.1 Строгие оценки погрешности
Способ, которым Рунге получил оценку погрешности,делаемой на одном шаге («локальной погрешности»), может быть описан следующимобразом. Для метода порядка /> рассмотрим локальнуюпогрешность
/> (2.7.1)
и воспользуемся ее тейлоровским разложением:
/>, (2.7.2)
где /> и />.Явное вычисление /> дает выражение вида
/>, (2.7.3)
где /> и /> содержатчастные производные /> до порядков /> и /> соответственно. Далеепоскольку />, имеем />. Такимобразом, если ограничены все частные производные /> допорядка /> включительно, имеем /> и/>. Следовательно, существует постоянная /> такая, что /> и
/>. (2.7.4)
Бибербах использовал несколько иной подход. Запишем
/> (2.7.5)
и воспользуемся тейлоровскими разложениями
/> (2.7.6)
Для векторных функций эти формулы справедливыпокомпонентно (возможно, с различным />). В силу условий порядкапервые члены разложения (2.6.5) по степеням /> обращаютсяв нуль. Таким образом, справедлива следующая теорема.
Теорема.
Если метод Рунге-Кутты (2.3.1) имеет порядок /> и если все частные производные /> до порядка /> включительносуществуют и непрерывны, то локальная погрешность метода (2.3.1) допускаетследующую строгую оценку:
/>, (2.7.7)
или
/>. (2.7.8)
Продемонстрируем этот результат, применяя к скалярномудифференциальному уравнению первый метод Рунге-Кутты (2.2.4), который имеетпорядок />. Дифференцируя (2.1.1), получим
/>. (2.7.9)
Вторая производная величины /> имеетвид
/>
Если условия теоремы выполнены, то легко видеть, чтовыражения (2.7.9) и (2.7.10) ограничены постоянной, которая не зависит от />, что и дает оценку (2.7.8).
1.7.2 Главный член погрешности
Для методов высших порядков строгие оценкипогрешностей, подобные (2.7.7), становятся очень непрактичными. Поэтому гораздоболее реалистично рассматривать первый ненулевой член в тейлоровским разложениипогрешности.
Теорема.
Если метод Рунге-Кутты имеет порядок /> и если /> непрерывнодифференцируема /> раз, то для главного членапогрешности имеем:
/>. (2.7.11)
/> (2.7.12)
1.7.3 Оценка глобальной погрешности
Глобальной (накопленной) погрешностью[3]называется погрешность численного решения после выполнения нескольких шагов. Пустьмы имеем некоторый одношаговый метод, с помощью которого при заданных начальныхданных /> и длине шага /> мыопределяем численное решение />, аппроксимирующее />. Воспользуемся обозначениями Хенричи дляэтого процесса:
/>, (2.7.13)
и назовем /> функцией приращения дляданного метода.
/>
Оценивание глобальной погрешности методами a) и b)
Тогда численное решение в точке /> получается с помощью пошаговой процедуры
, (2.7.14)
и наша задача состоит в оценке глобальной погрешности
/> (2.7.15)
Эта оценка находится простым способом: локальныепогрешности переносятся в конечную точку и затем складываются. Этот «переноспогрешностей» можно выполнить двумя разными способами:
a)перенося погрешность вдоль кривых точных решений; этот способ может датьхорошие результаты, если известны хорошие оценки распространения погрешностидля точных решений.
b)перенося погрешность />-го шага посредством выполнения /> шагов численного метода; этот способиспользовали в своих доказательствах Коши (1824) и Рунге (1905), он легкообобщается на многошаговые методы.
В обоих случаях оценим сначала локальные погрешности:
/>. (2.7.16)
Займемся теперь оценкой перенесенных погрешностей />.
a) Теорема.
Обозначим /> окрестность точки />, где /> – точноерешение уравнения
/>.
Пусть в /> справедливы оценкилокальных погрешностей (2.7.16) и выполнено одно из условий:
/> или />. (2.7.17)
Тогда имеет место следующая оценка глобальнойпогрешности (2.7.15):
/>, (2.7.18)
где />,
/>
и /> достаточно мало для того,чтобы численное решение оставалось в />.
Доказательство.
При /> оценка (2.7.18) переходитв />.
/>. (2.7.19)
Подставляя в неравенство
/>
выражение (2.7.18) с учетом (2.7.16) и принимая вовнимание, что />, приходим к такому неравенству:
/>.
Выражение в квадратных скобках мажорируется следующимиинтегралами:
/>, (2.7.20)
/>. (2.7.21)
Отсюда вытекает справедливость оценки (2.7.18).
b)При втором способе переноса погрешностей рассмотрим кроме (2.7.14) еще одночисленное решение, значения которого в соседних узлах связаны равенством
/>.
Оценим норму разности /> через />. Для /> формулыметода Рунге-Кутты запишем в следующих обозначениях:
/>
Вычитая из этих формул соответствующие формулы(2.3.1), получим для норм разностей такие оценки:
/>
/>
Оценивание римановых сумм методом a) и b)
Пусть /> – постоянная Липшица дляфункции /> и пусть />. Тогдафункция приращения /> для метода (2.3.1) удовлетворяетнеравенству
/>, (2.7.22)
где
/>. (2.7.23)
Из (2.7.22) получаем искомую оценку:
/>, (2.7.24)
и с её помощью оценку перенесенных погрешностей вместооценки (2.7.19).
Предположим, что для начальных значений, лежащих наточном решении, локальная погрешность удовлетворяет оценке
/> (2.7.25)
и что в окрестности решения функция приращения /> удовлетворяет неравенству
/>. (2.7.26)
Тогда для глобальной погрешности (2.7.15) справедливаследующая оценка:
/>, (2.7.27)
где />.
1.8 Оптимальный выбор шага
Предположим, что при интегрировании от точки /> до точки /> с шагом /> погрешность приближенно равна />. Так как это соответствует росту погрешностисо скоростью, приблизительно равной />, то />, где /> – функция,определяющая шаг. Положим /> и получим оценкуинтеграла, который приближенно равен полной погрешности:
/>
С другой стороны, затраты будут пропорциональны числушагов, которое приближенно равно
/>
Методами вариационного исчисления можно показать, чтоесли мы хотим минимизировать затраты /> при некоторомфиксированном значении погрешности />, то следует сохранятьпостоянной величину />. Это означает, что окончательнаяпогрешность должна быть одинаковой на каждом шаге.
В современных программах[4],реализующих методы Рунге-Кутты, обязательно используется некоторый алгоритмавтоматического изменения шага интегрирования. Интуитивно ясно, что на участкахплавного изменения решения счет можно вести с достаточно крупным шагом. В то жевремя, на тех участках, где происходят резкие изменения поведения решения,необходимо выбирать мелкий шаг интегрирования. Обычно начальное значение шагазадет пользователь или оно определено в программе. Далее шаг интегрирования изменяетсяв соответствии с величиной, получаемой в ходе вычисления оценки локальнойпогрешности.
Существует достаточно много способов оценки локальнойпогрешности, среди которых так называемое правило Рунге. Однако в моейпрограмме я использовал самый простой и в то же время эффективный способ оценкилокальной погрешности, который описан в разделе 3.1. «Описание программы Ilya RK-4 версия 1.43». Этот метод базируется на удвоении илиделении пополам длины шага в зависимости от отношения локальной погрешности имаксимально локальной допустимой погрешности />.
2. Практическая часть
2.1 Описание программы«IlyaRK-4 версия 1.43»Программа для нахождения интегральной кривой,удовлетворяющей поставленной задаче Коши написана на языке высокого уровня BorlandC++ 3.1. Программа состоит из четырех функций. При помощидиректив препроцессора #define определены максимальный шаг и величина локальноймаксимальной погрешности, а также номер версии программы. Рассмотрим подробнееработу программы в комплексе.
Функция title() предназначена для печати на экране названияпрограммы.
Функция do_step() совершает один шаг Рунге-Кутты и возвращаетполученное значение. В качестве входных параметров в нее передается текущееположение, значение искомой функции, вычисленное на предыдущем шаге и величинашага, с которым требуется произвести шаг.
Функция f() задает правую часть дифференциального уравнения,левая часть дифференциального уравнения равна />. Вкачестве аргументов функции передается /> и />.
Функция main() – основная функция программы. У пользователязапрашивается точка, начиная с которой необходимо отобразить решение задачиКоши, точка, правая граница интегрирования и значение в левой точке, черезкоторое обязана проходить искомая интегральная кривая. После этого программаначинает вычислительный процесс, выводя полученные значения на экран в видесписка и в текстовый файл “rk4.txt” на диск.После того, как будет достигнута правая граница интегрирования, процессостановится и пользователю будет предложено нажать на любую клавишу для того,чтобы построить график. Для построения графика программа переключается вграфический режим. График масштабируется с учетом того, чтобы он всегда былвиден на экране вне зависимости от того, как высоко или низко от оси абсцисс онлежит. Кроме того, программа постарается сделать его максимально наглядным. Дляэтого будут проведены пунктирные линии, соответствующие минимальному имаксимальному значению интегральной кривой, левому и правому концаминтегрирования, а также значению интегральной кривой в указанной точке />. Для того, чтобы пользователь мог легкоориентироваться на графике, рядом с пунктирными линиями пишутся координатныезначения с точностью до двух десятичных знаков. Как показали многочисленныетесты, проведенные на компьютере на базе процессора Intel Pentium 4B стактовой частотой 2.4 ГГц, построение графика происходит значительно быстрее,чем первичный расчет с выводом на экран и записью в файл. В этом легкоубедиться[5], если задать довольнобольшой отрезок интегрирования, например [-100,100].
Программа применяетследующий метод Рунге-Кутты четвертого порядка.
/>Для получения более точных результатов и возможно дажеувеличения скорости работы в определенных ситуациях программа используетавтоматическое управление длиной шага. Производятся одна итерация с шагом />, а затем две итерации подряд с шагом />, начиная с позиции x_cur.Разность между первым и последним упомянутым здесь результатами, взятая помодулю, считается локальной погрешностью. Если она больше, чем наперед заданнаявеличина />, то шаг делится пополам и тело циклаповторяется.
Для того, чтобы программамогла «разгоняться» после уменьшения шага, предусмотрено условие увеличениядлины шага. Оно состоит в том, что если погрешность между вторым из двухзначений, вычисленных с шагом /> и значением, вычисленнымс шагом />, не превосходит />, то шагувеличивается вдвое и тело цикла повторяется. Отметим, что величина шага неможет превосходить значения MAXSTEP,которое определяется директивой препроцессора #define. Если ни одно из двух описанных выше условий невыполняется, это означает, что шаг оптимален и программа производит вычислениезначения функции с записью его в файл и отображением на экране.
Программа снабженамеханизмом защиты от сбоев – в случае, если интегрируемая функция терпит разрыв(ее нельзя интегрировать на данном участке), программа останавливается ивыдается сообщение о невозможности продолжать. Работоспособность этогомеханизма проверена на некоторых разрывных функциях, таких как тангенс и др.
Заключение
В работе детально рассмотрен метод Рунге-Куттычетвертого порядка с автоматическим выбором длины шага, приведены необходимыетеоретические сведения, освещены альтернативные методы и их эффективность.
Был разработан алгоритм программного модуля,позволяющий автоматически менять величину шага интегрирования при решениизадачи Коши в зависимости от требуемой точности, что является непременнымтребованием, предъявляемым ко всем хорошим современным программам данногокласса, написано приложение, решены примеры.
Список использованных источников
[1].Амоносов А.А., Дубинский Ю.А., Копченова Н.В. «Вычислительные методы дляинженеров», М., Высшая школа, 1994, 544с.
[2].Хайрер Э., Нёрсетт С., Ваннер Г. «Решение обыкновенных дифференциальныхуравнений. Нежесткие задачи», М., Мир, 1990, 512с.
[3].Холл Д., Уатт Д. «Современные численные методы решения обыкновенныхдифференциальных уравнений», М., Мир, 1979, 312с.
Приложение А. Графики функций
В данном приложениирассмотрены три дифференциальных уравнения первого порядка. К каждому уравнениюприлагается по два графика – первый из них построен созданным приложением, авторой создан в пакете Maple9.01.
/>
/>
Интегральная кривая,построенная приложением «IlyaRK-4 версия 1.43»
/>
Интегральная кривая,построенная математическим пакетом WaterlooMaple9.01/>
/>
Интегральная кривая,построенная приложением «IlyaRK-4 версия 1.43»
/>
Интегральная кривая,построенная математическим пакетом WaterlooMaple9.01
/>
/>
Интегральная кривая,построенная приложением «IlyaRK-4 версия 1.43»
/>
Интегральная кривая,построенная математическим пакетом WaterlooMaple9.01
/>
/>
Интегральная кривая,построенная приложением «IlyaRK-4 версия 1.43»
/>
Интегральная кривая,построенная математическим пакетом WaterlooMaple9.01
Приложение Б. Пример таблицы значений функции y(x)
y(-7)=100, h=0.4
y(-6.6)=100.045,h=0.4
y(-6.2)=100.112,h=0.4
y(-5.8)=100.212,h=0.4
y(-5.4)=100.361,h=0.4
y(-5)=100.585,h=0.4
y(-4.6)=100.919,h=0.4
y(-4.2)=101.419,h=0.4
y(-3.8)=102.17,h=0.4
y(-3.4)=103.301,h=0.2
y(-3.2)=104.067,h=0.2
y(-3)=105.011,h=0.2
y(-2.8)=106.175,h=0.2
y(-2.6)=107.615,h=0.2
y(-2.4)=109.4,h=0.2
y(-2.2)=111.62,h=0.2
y(-2)=114.392,h=0.2
y(-1.8)=117.873,h=0.2
y(-1.6)=122.267,h=0.2
y(-1.4)=127.857,h=0.2
y(-1.2)=135.033,h=0.2
y(-1)=144.346,h=0.2
y(-0.8)=156.596,h=0.2
y(-0.6)=172.977,h=0.1
y(-0.5)=183.256,h=0.1
y(-0.4)=195.327,h=0.1
y(-0.3)=209.595,h=0.1
y(-0.2)=226.578,h=0.1
y(-0.1)=246.953,h=0.05
y(-0.05)=258.68,h=0.05
y(2.96985e-15)=271.608,h=0.05
y(0.05)=285.897,h=0.05
y(0.1)=301.73,h=0.05
y(0.15)=319.32,h=0.05
y(0.2)=338.919,h=0.05
y(0.25)=360.821,h=0.05
y(0.3)=385.374,h=0.05
y(0.35)=412.989,h=0.05
y(0.4)=444.156,h=0.05
y(0.45)=479.459,h=0.025
y(0.475)=498.877,h=0.025
y(0.5)=519.603,h=0.025
y(0.525)=541.747,h=0.025
y(0.55)=565.433,h=0.025
y(0.575)=590.794,h=0.025
y(0.6)=617.978,h=0.025
y(0.625)=647.149,h=0.025
y(0.65)=678.489,h=0.025
y(0.675)=712.199,h=0.025
y(0.7)=748.501,h=0.025
y(0.725)=787.645,h=0.025
y(0.75)=829.906,h=0.025
y(0.775)=875.592,h=0.025
y(0.8)=925.047,h=0.025
y(0.825)=978.656,h=0.025
y(0.85)=1036.85,h=0.025
y(0.875)=1100.11,h=0.025
y(0.9)=1168.98,h=0.025
y(0.925)=1244.07,h=0.0125
y(0.9375)=1284.16,h=0.0125
y(0.95)=1326.08,h=0.0125
y(0.9625)=1369.91,h=0.0125
y(0.975)=1415.77,h=0.0125
y(0.9875)=1463.78,h=0.0125
y(1)=1514.04,h=0.0125
y(1.0125)=1566.7,h=0.0125
y(1.025)=1621.89,h=0.0125
y(1.0375)=1679.75,h=0.0125
y(1.05)=1740.44,h=0.0125
y(1.0625)=1804.13,h=0.0125
y(1.075)=1871,h=0.0125
y(1.0875)=1941.23,h=0.0125
y(1.1)=2015.04,h=0.0125
y(1.1125)=2092.63,h=0.0125
y(1.125)=2174.24,h=0.0125
y(1.1375)=2260.12,h=0.0125
y(1.15)=2350.54,h=0.0125
y(1.1625)=2445.79,h=0.0125
y(1.175)=2546.16,h=0.0125
y(1.1875)=2652,h=0.0125
y(1.2)=2763.65,h=0.0125
y(1.2125)=2881.49,h=0.0125
y(1.225)=3005.94,h=0.0125
y(1.2375)=3137.43,h=0.0125
y(1.25)=3276.43,h=0.0125
y(1.2625)=3423.46,h=0.0125
y(1.275)=3579.07,h=0.0125
y(1.2875)=3743.83,h=0.0125
y(1.3)=3918.41,h=0.0125
y(1.3125)=4103.47,h=0.0125
y(1.325)=4299.77,h=0.0125
y(1.3375)=4508.11,h=0.0125
y(1.35)=4729.35,h=0.00625
y(1.35625)=4845.11,h=0.00625
y(1.3625)=4964.45,h=0.00625
y(1.36875)=5087.5,h=0.00625
y(1.375)=5214.41,h=0.00625
y(1.38125)=5345.3,h=0.00625
y(1.3875)=5480.34,h=0.00625
y(1.39375)=5619.66,h=0.00625
y(1.4)=5763.44,h=0.00625
y(1.40625)=5911.83,h=0.00625
y(1.4125)=6065,h=0.00625
y(1.41875)=6223.14,h=0.00625
y(1.425)=6386.44,h=0.00625
y(1.43125)=6555.09,h=0.00625
y(1.4375)=6729.28,h=0.00625
y(1.44375)=6909.25,h=0.00625
y(1.45)=7095.2,h=0.00625
y(1.45625)=7287.36,h=0.00625
y(1.4625)=7485.99,h=0.00625
y(1.46875)=7691.33,h=0.00625
y(1.475)=7903.64,h=0.00625
y(1.48125)=8123.19,h=0.00625
y(1.4875)=8350.28,h=0.00625
y(1.49375)=8585.2,h=0.00625
y(1.5)=8828.27, h=0.00625
Приложение В. Листинг программы «Ilya RK-4 версия 1.43»
//----------------------------------------------------------------------- //
#include
#include
#include
#include
#include
#include
#define EPSILON 0.00001
#define MAXSTEP 1
#define VERSION 1.43
//----------------------------------------------------------------------- //
double f(double x, double y);
double do_step(double h, double x_cur, doubley_cur);
void title(void);
void main(void);
//----------------------------------------------------------------------- //
double f(double x, double y)
{
//Правая часть ДУ f(x,y)
return (pow(2.718,x)*y);
}
//----------------------------------------------------------------------- //
void main(void)
{
inti; // Вспомогательный счетчик
intmetka; // Метка на осях
intflag = 0; // Флаг правого конца интегрирования
intmetka1, metka2; // Переменные меток на осях координат
double err= 0; // Погрешность
doublex0, y0; // Координаты точки начального условия
double big2_step_res, super_step_res; // Результатыдлинныхшагов
doublek = 1; // Коэффициент коррекции
doublezoom = 1; // Масштаб на графике
double big_step_res, small_step_res; // Результатышаговинтегрирования
double a,b; // Границы
doubletemp; // Переменная для служебных нужд
doublex_cur, y_cur; // Переменные метода РК
doubleh; // Шаг интегрирования
doublef_max = 0, f_min = 0; // Максимальное и минимальное значение кривой
doublenorma = 0; // Норма (для корректного масштабирования графика)int c = 8; //Переменная цвета разделительных линий
FILE*myfile; // Указатель на текстовый файл с таблицей значений
//Инициализируем графический адаптер
int gdriver= DETECT, gmode,errorcode;
initgraph(&gdriver, &gmode, "");
errorcode = graphresult();
if (errorcode != grOk)
{
printf(«Ошибкаинициализацииграфики:%s\n», grapherrormsg(errorcode));
getch();
}
textcolor(0);
setbkcolor(0);
title();
printf(«y'=f(x,y), y(x0)=y(a)=y0, [a,b] — отрезокинтегрирования\n»);
label1: printf("\na=");
scanf("%lg", &a);
printf(«b=»);
scanf("%lg", &b);
//Авто смена границ при необходимости
if (a> b)
{
temp = a;
a = b;
b = temp;
}
if (a== b)
{
printf(«Началоотрезка интегрирования совпадает с его концом, повторите ввод!\n»);
gotolabel1;
}
printf(«y(%lg)=», a);
scanf("%lg", &y0);
title();
printf("[%lg,%lg] — границыинтегрирования,y(%lg)=%lg — начальное условие.\n",a, b, a, y0);
//Инициализация
h= fabs(b — a) / 10;
if (h > 0.1) h = 0.1;
x_cur = a;
y_cur = y0;
f_max = y_cur;
f_min = y_cur;
myfile = fopen(«rk4.txt», «w»);
fprintf(myfile, «Program: Ilya RK4 Version%g\n», VERSION);
fprintf(myfile, «Method: Runge-Kutta\n»);
fprintf(myfile, «The order of method:4\n»);
fprintf(myfile, «Automatic integration stepselect: Enabled\n»);
fprintf(myfile, "[a,b]=[%lg,%lg],y(%lg)=%lg\n", a, b, a, y0);
while (x_cur
{
if (flag > 1) break;
big_step_res = do_step(h, x_cur, y_cur);
temp = do_step(h / 2, x_cur, y_cur);
small_step_res = do_step(h / 2, x_cur + h / 2,temp);
err = fabs(big_step_res — small_step_res);
// Уменьшение длины шага
if(err > EPSILON)
{
h= h / 2;
continue;
}
//Увеличение длины шага
big2_step_res= do_step(h, x_cur + h, big_step_res);
super_step_res = do_step(2 * h, x_cur, y_cur);
if (fabs(big2_step_res — super_step_res)
{
h *= 2;
continue;
}
if (h > MAXSTEP) h = MAXSTEP;
// Защита от сбоев
if(h
{
printf(«Ошибка!Возможно, функция разрывна.\nПроинтегрировать на данном интервале невозможно.Скорее всего, g(%lg)=», x_cur);
fprintf(myfile,«Ошибка! Возможно, функция разрывна.\nПроинтегрировать на данном интерваленевозможно. Скорее всего,g(%lg)=», x_cur);
if (y_cur
{
printf("-oo.\n");
fprintf(myfile, "-oo.\n");
}
else
{
printf("+oo.\n");
fprintf(myfile, "+oo.\n");
}
getch();
fclose(myfile);
exit(1);
}
printf(«y(%lg)=%lg, err=%lg, h=%lg\n»,x_cur, y_cur, err, h);
if (y_cur
if (y_cur > f_max) f_max = y_cur;
fprintf(myfile, «y(%lg)=%lg, h=%lg\n»,x_cur, y_cur, h);
if (x_cur + h > b) h = fabs(b — x_cur);
x_cur += h;
y_cur = big_step_res;
if (x_cur >= b) flag++;
}
fclose(myfile);
printf("\nТаблицазначенийзаписанавфайлrk4.txt.\n");
printf("\nНажмителюбую клавишу для построения графика...");
flag = 0;
getch();
// Построениеграфика
cleardevice(); clrscr();
if (fabs(a) > fabs(b)) zoom = fabs(getmaxx() / 2/ a);
else zoom = fabs(getmaxx() / 2 / b);
// Рисуемграницы
for (i = 0; i
{
if (c == 8) c = 0;
else c = 8;
setcolor(c);
line(a * zoom + getmaxx() / 2, i, a * zoom +getmaxx() / 2, i + 5);
line(b * zoom + getmaxx() / 2 — 1, i, b * zoom +getmaxx() / 2 — 1, i + 5);
}
if (fabs(f_min) > fabs(f_max)) norma =fabs(f_min) * zoom;
else norma = fabs(f_max) * zoom;
//Определение коэффициента коррекции
k = (getmaxy()/ 2) / norma;
//Предотвращение чрезмерного масштабирования
if (k
if (k > 10000) k = 10000;
for (i = 0; i
{
if (c == 8) c = 0;
else c = 8;
setcolor(c);
line(i, -y0 * zoom * k + getmaxy() / 2, i + 5, -y0* zoom * k + getmaxy() / 2);
line(i, -f_min * zoom * k + getmaxy() / 2, i + 5,-f_min * zoom * k + getmaxy() / 2);
line(i, -f_max * zoom * k + getmaxy() / 2, i + 5,-f_max * zoom * k + getmaxy() / 2);
}
metka = ceil((-y0 * zoom * k + getmaxy() / 2) / 16);
if (metka
if (metka == 15) metka = 16;
if (metka > 25) metka = 25;
gotoxy(1, metka);
printf(«Y=%.2g», y0, metka);
metka = ceil((-f_max * zoom * k + getmaxy() / 2) /16);
if (metka
if (metka == 15) metka = 16;
if (metka > 25) metka = 25;
gotoxy(1, metka);
printf(«Y=%.2lg», f_max, metka);
metka = ceil((-f_min * zoom * k + getmaxy() / 2) /16);
if (metka
if (metka == 15) metka = 16;
if (metka > 25) metka = 25;
gotoxy(1, metka);
printf(«Y=%.2lg», f_min, metka);
//Пишем границы, делаем отметки на осях координат
metka1 = ceil((a * zoom + getmaxx() / 2) / 8);
if (metka1
if (metka1 > 75) metka1 = 75;
if (metka == 17) metka = 18;
gotoxy(metka1, 15);
if (a != 0) printf("%.2lg", a);
metka2 = ceil((b * zoom + getmaxx() / 2 — 1) / 8);
if (metka2 — metka1
if (metka2
if (metka2 > 75) metka2 = 75;
gotoxy(metka2, 15);
printf("%.2lg", b);
gotoxy(80, 17);
printf(«X»);
gotoxy(42,1);
printf(«Y»);
gotoxy(39, 15);
printf(«0»);
// Рисуемсистемукоординат
setcolor(15);
line(0, getmaxy() / 2, getmaxx(), getmaxy() / 2);
line(getmaxx() / 2, 0, getmaxx() / 2, getmaxy());
line(getmaxx() / 2, 0, getmaxx() / 2 — 5, 10);
line(getmaxx() / 2, 0, getmaxx() / 2 + 5, 10);
line(getmaxx(), getmaxy() / 2, getmaxx() — 10,getmaxy() / 2 + 5);
line(getmaxx(), getmaxy() / 2, getmaxx() — 10,getmaxy() / 2 — 5);
setcolor(10);
h = fabs(b — a) / 10;
if (h > 0.1) h = 0.1;
y_cur = y0;
x_cur = a;
f_max = y_cur;
f_min = y_cur;
x0 = zoom * a + getmaxx() / 2;
y0 = (zoom * (-y_cur)) * k + getmaxy() / 2;
while (x_cur
{
if (flag > 1) break;
big_step_res = do_step(h, x_cur, y_cur);
temp = do_step(h / 2, x_cur, y_cur);
small_step_res = do_step(h / 2, x_cur + h / 2,temp);
err = fabs(big_step_res — small_step_res);
if (err > EPSILON)
{
h = h / 2;
continue;
}
big2_step_res = do_step(h, x_cur + h,big_step_res);
super_step_res = do_step(2 * h, x_cur, y_cur);
if (fabs(big2_step_res — super_step_res)
{
h *= 2;
continue;
}
if (h > MAXSTEP) h = MAXSTEP;
line (x0, y0, zoom * x_cur + getmaxx() / 2, zoom *(-y_cur) * k + getmaxy() / 2);
x0 = zoom * (x_cur) + getmaxx() / 2;
y0 = (zoom * (-y_cur)) * k + getmaxy() / 2;
if (x_cur + h > b) h = fabs(b — x_cur);
x_cur += h;
y_cur = big_step_res;
if (x_cur >= b) flag++;
}
while (getch() != 0);
}
//----------------------------------------------------------------------- //
void title(void)
{
//Печать заголовка программы
cleardevice();clrscr();
printf("Решение дифференциальных уравнений методом Рунге-Кутты 4-го порядка\n");
printf("с автоматическим выбором длины шага\n");
printf("Разработал Щербаков Илья, гр. 520212, версия %g\n", VERSION);
printf("____________________________________________________\n");
}
//----------------------------------------------------------------------- //
double do_step(double h, double x_cur, double y_cur)
{
double k1, k2, k3, k4, delta_y_cur;
k1 = f(x_cur, y_cur);
k2 = f(x_cur + (h / 2), y_cur + (h / 2) * k1);
k3 = f(x_cur + (h / 2), y_cur + (h / 2) * k2);
k4 = f(x_cur + h, y_cur + h * k3);
delta_y_cur = (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
return(y_cur + delta_y_cur);
}
//----------------------------------------------------------------------- //