Solve_ODE (ЭВМ для спецгруппы)
Описание файла
Файл "Solve_ODE" внутри архива находится в папке "ЭВМ для спецгруппы". Документ из архива "ЭВМ для спецгруппы", который расположен в категории "". Всё это находится в предмете "практика расчётов на пэвм" из 2 семестр, которые можно найти в файловом архиве МГУ им. Ломоносова. Не смотря на прямую связь этого архива с МГУ им. Ломоносова, его также можно найти и в других разделах. .
Онлайн просмотр документа "Solve_ODE"
Текст из документа "Solve_ODE"
2. Численные методы решения задачи Коши
2.1. Ошибки и устойчивость метода
Для простоты начнем с задачи Коши для одного дифференциального уравнения:
; (1)
.
Запишем решение в виде разложения в ряд Тейлора относительно начальной точки x0:
(2)
Пусть точка x, в которой ищется решение, отстоит от начальной точки на расстояние h:
Если h достаточно мало, то в разложении (2) можно ограничиться линейными членами, внеся ошибку порядка h2:
. (3)
Уравнение (3) соответствует методу Эйлера — простейшему методу численного интегрирования дифференциальных уравнений. Имея начальное условие y(x0)=y0, мы можем с помощью этого уравнения получить приближенное решение в точке x1=x0+h. Принимая найденное решение за начальное условие для следующего шага, получим решение в точке x2=x1+h=x0+2h. Многократно повторяя эту процедуру, будем иметь таблицу значений y(x) для x, меняющегося от x0 до некоторого конечного значения xm с шагом h. Однако для практического применения подобного метода решающее значение имеет поведение ошибки, которая накапливается от шага к шагу.
Существуют два источника ошибок, искажающих конечный результат. Во‑первых, это присущая методу интегрирования ошибка усечения (или дискретизации), которая связана с отбрасыванием членов высших порядков в разложении (2). Величина ошибки усечения зависит от шага интегрирования h; в частности, для метода Эйлера она имеет порядок O(h2). Во‑вторых, в любых реальных вычислениях присутствует ошибка округления, вызванная ограниченной точностью представления чисел. Ошибка округления не является специфической для того или иного метода интегрирования и не зависит от h; более того, ее можно произвольно менять, выполняя вычисления с большим или меньшим количеством значащих цифр. Если речь идет о машинных расчетах, то относительная ошибка округления 
не превышает некоторого фиксированного уровня, определяемого формой представления чисел с плавающей точкой на данной ЭВМ. Так, наиболее распространенные форматы обычной и удвоенной точности (числа размером 4 и 8 байтов, соответственно) характеризуются величинами 1.19·107 и 2.22·1016. Переход к вычислениям с увеличенной точностью позволяет сократить порядок ошибки округления по меньшей мере вдвое, хотя при этом обычно увеличивается время выполнения арифметических операций.
Ошибка, сделанная на одном шаге интегрирования, называется локальной ошибкой. Она характеризует погрешность решения в точке xn+1=xn+h, если в исходной точке xn значение yn считается точным. Как известно, дифференциальное уравнение определяет целое семейство решений, а начальное условие выделяет из этого семейства конкретное решение. Наличие локальной ошибки приводит к тому, что точка (xn+1,yn+1), где yn+1 вычислено по формуле Эйлера
, (4)
уже не принадлежит решению, проходящему через точку (xn,yn). Однако, будучи принятой в качестве начального условия следующего шага, точка (xn+1,yn+1) соответствует некоторому другому члену семейства решений. Иными словами, локальная ошибка на каждом шаге вызывает переход с одного точного решения на другое.
Чтобы оценить фактическую погрешность численного интегрирования, надо сравнить полученный результат с точным решением, проходящим через начальную точку (x0,y0). Пусть yn – точное решение, отвечающее значению независимой переменной xn, т. е. точки (x0,y0) и (xn,yn) принадлежат одному и тому же члену семейства решений дифференциального уравнения. Приближенное решение, полученное на n‑м шаге численного интегрирования, обозначим 
. Отклонение от yn называют глобальной ошибкой:
.
Каждый шаг добавляет к накопленной глобальной ошибке собственную локальную ошибку и дополнительную составляющую, обусловленную искажением начального условия из‑за ошибок на предыдущих шагах. Связь между локальной и глобальной ошибками метода Эйлера показана на рис. 1. Вообще говоря, глобальная ошибка на n‑м шаге может быть как больше, так и меньше суммы локальных ошибок в зависимости от характера решений дифференциального уравнения.
Рассмотрим поведение глобальной ошибки метода Эйлера. Для упрощения анализа будем считать, что ошибки округления отсутствуют, т. е. все арифметические действия выполняются абсолютно точно. Тогда локальная ошибка обусловлена только ошибкой усечения и имеет величину порядка O(h2). На (n+1)‑м шаге приближенное решение 
получается по формуле (4) из результата , достигнутого после n шагов, что можно записать следующим образом:
, (5)
С другой стороны, для точного решения имеет место соотношение
. (6)
Вычитая (6) из (5), получим:
.
П
Рис. 1. Локальная и глобальная ошибки при численном интегрировании дифференциального уравнения методом Эйлера.
ренебрегая членами высших порядков относительно n, можно записать:
Подставляя это выражение в предыдущее уравнение и обозначая 
, получим окончательно:
. (7)
Из (7) следует, что поведение глобальной ошибки зависит от знака и величины шага h. Для простоты рассмотрим случай линейного дифференциального уравнения, т. е. =const. При >0 независимо от выбора h множитель перед n больше единицы, т. е. глобальная ошибка с каждым шагом будет экспоненциально возрастать. Иная ситуация наблюдается в случае <0. Если 
, то n умножается на коэффициент, меньший единицы, так что вклады от ошибок предыдущих шагов затухают. Глобальная ошибка не возрастает и определяется, по существу, локальной ошибкой последнего шага. При увеличении шага выше порогового значения 
множитель (1+h) в (7), будучи отрицательным, по абсолютной величине станет больше единицы. Это значит, что глобальная ошибка начнет неограниченно расти, причем знак ее на каждом шаге будет меняться, то есть возникнет «раскачка» решения.
Тенденция к неограниченному росту глобальной ошибки называется неустойчивостью численного решения (часто также говорят о неустойчивости метода интегрирования). Сам факт неустойчивости еще не означает безусловной неприемлемости решения. Так, в случае >0 глобальная ошибка метода Эйлера экспоненциально растет, причем скорость роста (показатель экспоненты) зависит от h. Точное решение дифференциального уравнения также растет по экспоненциальному закону (простейшее уравнение такого типа имеет вид 
с решением 
). Поэтому при достаточно малом h относительная погрешность решения 
может оставаться постоянной или даже уменьшаться по мере роста n. В таких случаях говорят об относительной устойчивости. Заметим, что при <0 y(x) является экспоненциально убывающей функцией, так что здесь нельзя получить относительно устойчивое решение, если глобальная ошибка начинает расти. В дальнейшем нас будет интересовать именно случай убывающих решений как наиболее характерный для дифференциальных уравнений химической кинетики (неограниченно возрастающие решения несовместимы с физическим смыслом рассматриваемых моделей).
2.2. Неявная схема численного интегрирования
Вывод уравнения метода Эйлера в предыдущем разделе опирался на разложение решения в ряд Тейлора в окрестности точки x0. Однако, с таким же успехом можно построить разложение относительно конечной точки шага x=x0+h:
Ограничиваясь двумя первыми членами разложения, получаем формулу, аналогичную (3):
, (8)
откуда следует схема численного интегрирования
. (9)
Здесь подлежащая определению величина yn+1 входит как в левую, так и в правую часть уравения, в отличие от формулы Эйлера (4), где yn+1 явно выражается через известные величины xn и yn.
Схемы построения численного решения, подобные формуле (9), называют неявными в отличие от явных схем типа (4). В частности, (9) иногда называют неявным методом Эйлера.
В общем случае при произвольной функции F(x,y) вычисления по неявной схеме обычно организуют по принципу итерации. Вначале каким-либо способом получают оценку (начальное приближение) неизвестной величины yn+1, которую подставляют в правую часть уравнения (9) и вычисляют результат, который считают уточненным значением yn+1. В зависимости от качества начального приближения итерация может быть однократной либо повторяться несколько раз до достижения сходимости. На практике обычно комбинируют явную и неявную схемы, получая с помощью первой из них оценку (прогноз) решения в точке (n+1), а затем применяя вторую схему для уточнения (коррекции) результата. Методы такого типа называются методами прогноза и коррекции.
Неявные схемы численного интегрирования отличаются от явных с точки зрения устойчивости. Рассмотрим, например, устойчивость решений, получаемых неявным методом Эйлера. Записав по аналогии с (5) уравнение (9) с учетом присутствующих ошибок,
,
и вычитая соотношение для точного решения
,
а затем представляя 
по формуле Тейлора в линейном приближении
,
получим связь между глобальными ошибками n‑го и (n+1)‑го шагов:
, (10)
где, как и раньше, использовано обозначение 
.
Из соотношения (10) следует, что для дифференциального уравнения с <0 неявная схема интегрирования обеспечивает устойчивость решения при любой величине шага! Действительно, в данном случае знаменатель в (10) всегда больше единицы, так что влияние ошибок, сделанных на предыдущих шагах, постепенно затухает. Устойчивость, которая не нарушается при сколь угодно больших значениях h (в области отрицательных ), получила название A‑устойчивости. Это, однако, не значит, что величина шага интегрирования не имеет значения. Факт устойчивости говорит лишь о том, что погрешность решения приблизительно соответствует текущей локальной ошибке, которая увеличивается с ростом h (слагаемое O(h2) в числителе формулы (10)).
