Основы решения дифференциальных уравнений
Дифференциальные уравнения очень часто встречаются при построении моделей динамики объектов исследования. Они описывают, как правило, изменение параметров объекта во времени (хотя могут быть и другие случаи). Результатом решения дифференциальных уравнений являются функции, а не числа, как при решении конечных уравнений, вследствие чего методы решения их более трудоемки. Особенно это касается дифференциальных уравнений в частных производных. В данном пособии такие методы не рассматриваются. Владение методами решения дифференциальных уравнений обязательно при моделировании и очень важно для получения правильного результата.
Концепция решения дифференциальных уравнений
При использовании численных методов решение дифференциальных уравнений dy/dx = f(x,y) или у' = f(x,y) представляется в табличном виде, т.е. получается совокупность значений yi и хi. Решение носит шаговый характер, т.е. по одной или по нескольким начальным точкам (х, у) за один шаг находят следующую точку, затем следующую и т.д. Разница между двумя соседними значениями аргумента h =хi+1 - хi называется шагом.
Наибольшее распространение имеют задачи Коши, в которых заданы начальные условия: при х = x0 у(х0) = у0. Имея их, легко начинать процесс решения, т.е. найти у1 при x1, y2 - при x2 и т.д. Задачи другого типа - краевые задачи (например, с конечными условиями или с условиями в промежуточной точке) - решаются специальными приемами, в том числе нередко сведением к другим эквивалентным задачам с начальными условиями.
Выделяют два класса методов решения: одношаговые и многошаговые. Первый класс методов требует для нахождения следующего значения функции только одной текущей точки, т.е. yi+1 = F[f(xi, yi)], а второй - нескольких, например, уi+1 = F(yi-3, yi-2, yi-1, yi). Поэтому методы второго класса не обладают свойством "самостартования", т.е. ими нельзя начать решение задачи Коши, это всегда делается одношаговыми методами. К недостаткам многошаговых методов относится также и невозможность изменения в процессе решения величины шага (так как они используют предыдущие точки с ранее применяемым шагом, а учет меняющегося шага очень сложен и громоздок), что бывает необходимо для повышения эффективности метода. Заметим, что величина шага существенно влияет на точность и скорость решения, поэтому изменение ее в процессе решения - увеличение при медленно изменяющемся решении и уменьшение при быстро изменяющемся - очень важно для эффективности решения. К достоинствам многошаговых методов относят в основном меньший объем памяти компьютера, требующейся для реализации, и возможность теоретической оценки погрешности решения. Представителем класса многошаговых методов являются методы прогноза и коррекции. К классу одношаговых методов относятся методы Эйлера, Рунге - Кутта и др.
Основная идея получения простейших вычислительных алгоритмов в одношаговых методах сводится к разложению искомого решения у(х) в ряд Тейлора в окрестности текущей точки и усечения его. Количество оставленных членов ряда определяет порядок и, следовательно, точность метода. По полученному разложению, зная значение у в точке разложения уi и производную f(xi, yi), находят значение функции у через шаг h: уi+1 = уi + уi. Если в разложении удерживается большее число членов, то необходимо рассчитывать f(xi, yi) в нескольких точках (таким способом избегают необходимости прямого вычисления высших производных, присутствующих в разложении в ряд Тейлора). Расчетные алгоритмы многошаговых методов базируются на построении интерполяционных или аппроксимирующих функций, от которых берется интеграл.
Численными методами решаются не только отдельные уравнения, но и системы уравнений (чаще всего первого порядка), причем большинство методов решения одного уравнения легко распространяются на решение систем. Дифференциальные уравнения высших порядков вида
y(n) = f(x, y, y', y'', …, y(n-1))
решаются в основном сведением к системе уравнений первого порядка путем замены переменных: у1 = у', у2 = у'', у3 = у''' и т.д. При этом дифференциальное уравнение n-го порядка заменяется системой из n уравнений:
Основные методы
Метод Эйлера и модифицированный метод Эйлера
Метод Эйлера является одним из самых простых методов решения дифференциальных уравнений первого порядка y' = f(x, y). Он в основном используется как учебный, в практических расчетах он дает значительную погрешность. Вычислительный алгоритм представляется следующим образом:
yi+1 = y(xi + h) = yi + D yi = yi + hf(xi, yi),
где h - шаг по х (в общем случае может быть непостоянным). Запускается метод из начальных условий у(х0) = у0.
yi+1 = yi + 0,5h [f(xi, yi) + f(xi+1, yi+1)].
Здесь в формуле используется значение f(xi+1, уi+1) с еще пока неизвестным значением уi+1. Это значение может быть найдено предварительно, например, по методу Эйлера, а затем использовано в алгоритме. Если же выражение f(x, y) несложное, то можно выразить уi+1 из уравнения в явном виде и найти его или решить его относительно уi+1 численными методами.
Точность вычислений обычно контролируют двойным просчетом: сначала вычисляют решение уравнения на каком-то текущем шаге h, т.е. находясь в точке хi и вычисляя значение 
, затем в эту же точку хi+1 приходят за два шага по h/2, получают 
, сравнивают их: если для обоих вариантов различие 
в пределах желаемой погрешности, то решение принимают, а если нет, то опять делят шаг на два и т. д., до тех пор, пока не получится приемлемый результат. Однако следует помнить, что при очень маленьком шаге, получающемся в результате его последовательного деления, может значительной оказаться накапливающаяся вычислительная ошибка.
Метод Рунге-Кутта
Существует целая группа методов Рунге - Кутта (в последние годы начинает в литературе появляться "русский" вариант произношения фамилий авторов метода, в соответствии с которым название звучит как метод Рунге - Кутты), среди которых наибольшее распространение получил метод четвертого порядка. Следовательно, он более точен, чем метод Эйлера, который является методом первого порядка. Для расчета одного значения функции необходимо четыре раза вычислять правую часть дифференциального уравнения, а не два, как в модифицированном методе Эйлера второго порядка. Вычислительный алгоритм записывается следующим образом
где
Здесь также для контроля точности можно применять прием двойного просчета.
Метод Милна
Метод Милна относится к многошаговым методам и представляет один из методов прогноза и коррекции. Решение в следующей точке находится в два этапа. На первом этапе осуществляется по специальной формуле прогноз значения функции, а затем на втором этапе - коррекция полученного значения. Если полученное значение у после коррекции существенно отличается от спрогнозированного, то проводят еще один этап коррекции. Если опять имеет место существенное отличие от предыдущего значения (т.е. от предыдущей коррекции), то проводят еще одну коррекцию и т.д. Однако очень часто ограничиваются одним этапом коррекции.
Блок-схема метода прогноза и коррекции представлена на рисунке.
Метод Милна имеет следующие вычислительные формулы:
а) этап прогноза:
где для компактности записи использовано следующее обозначение fi = f(xi, yi);
б) этап коррекции:
Абсолютная погрешность определяется по формуле 
Метод требует несколько меньшего количества вычислений (например, достаточно только два раза вычислить f(x, y), остальные запомнены с предыдущих этапов), но требует дополнительного "расхода" памяти. Кроме этого, как уже указывалось выше, невозможно "запустить" метод: для этого необходимо предварительно получить одношаговыми методами первые три точки.
