Численные методы решения задачи Коши
Обыкновенным дифференциальным уравнением называется уравнение, связывающее независимую переменную [math]x[/math] , неизвестную функцию [math]y(x)[/math] этой независимой переменной и ее производные [math]y'(x),y''(x),\ldots,y^(x)\colon[/math]
где [math]F(x,y, y',\ldots,y^)[/math] — функция указанных аргументов, заданная в некоторой области их изменения.
Порядком дифференциального уравнения называется порядок старшей производной, входящей в уравнение.
Если в соотношении (6.1) функция [math]F[/math] такова, что его можно представить в виде
то уравнение называется обыкновенным дифференциальным уравнением л-го порядка, разрешенным относительно старшей производной.
Уравнение называется линейным, если функция [math]f[/math] линейна относительно искомой функции и ее производных, т.е. если уравнение может быть записано в виде
где [math]a_n(x),a_(x),\ldots,a_0(x),f(x)[/math] — известные в общем случае нелинейные функции от [math]x[/math] .
Решением дифференциального уравнения n-го порядка называется функция [math]y(x)[/math] , непрерывная на некотором интервале [math](a,b)[/math] вместе со своими производными до [math](n-1)[/math] порядка включительно, имеющая производную [math]y^(x)[/math] и такая, что подстановка [math]y(x)[/math] в уравнение обращает его в тождество.
График решения дифференциального уравнения называется интегральной кривой.
Одной из важнейших задач в теории и приложениях дифференциальных уравнений является задана Коши (начальная задана), в которой требуется найти решение дифференциального уравнения, удовлетворяющее заданным начальным условиям. Для уравнения (6.2) она записывается следующим образом:
y_0,y'_0,\ldots,y_0^[/math] — заданные числа.
Теорема 6.1 (о существовании и единственности решения задачи Коши (6.4)). Пусть выполнены следующие условия:
а) функция [math]f(x,y,\ldots,y^)[/math] определена и непрерывна в некоторой замкнутой области [math]\overline[/math] , а также имеет в [math]\overline[/math] ограниченные частные производные по переменным [math]y,y',\ldots,y^[/math] ;
б) точка [math](x_0,y_0,y'_0,\ldots,y_0^)[/math] лежит внутри области [math]\overline[/math] .
Тогда решение задачи Коши (6.4) существует и единственно.
Общим решением дифференциального уравнения л-го порядка в области [math]G\subset \overline[/math] ( [math]\overline[/math] — область, в которой выполнены условия теоремы 6.1) называется функция [math]y=y(x,C_1,\ldots,C_n)[/math] , зависящая от [math]n[/math] произвольных постоянных, и такая, что при подстановке в уравнение она обращает его в тождество при любых значениях [math]C_1,\ldots,C_n[/math] . Геометрически общее решение в области [math]G[/math] представляет собой семейство непересекающихся интегральных кривых, полностью покрывающих всю область.
Общим интегралом дифференциального уравнения называется соотношение вида [math]\varphi(x,y,C_1,\ldots,C_n)=0[/math] , неявно определяющее общее решение.
При конкретных значениях [math]C_1,\ldots,C_n[/math] , включая [math]\pm\infty[/math] , из общего решения выделяется частное решение, а общий интефал становится частным интегралом. В каждой точке [math](x,y)[/math] частного решения или частного интефала выполняются условия теоремы 6.1.
Наряду с проблемой решения дифференциальных уравнений л-го порядка на практике возникает проблема решения систем обыкновенных дифференциальных уравнений первого порядка, связывающих независимую переменную [math]x[/math] , неизвестные функции [math]y_1(x),\ldots, y_n(x)[/math] и их производные [math]y'_1(x),\ldots, y'_n(x)[/math] .
В случае, если уравнения разрешимы относительно производных, систему можно записать в нормальной форме Коши (где [math]f_i(x,y_1,\ldots,y_n),
i=\overline[/math] — известные функции):
Решением системы (6.5) называется совокупность [math]n[/math] функций [math]y_1(x),\ldots, y_n(x)[/math] , непрерывных на некотором интервале (д,6), такая, что подстановка этих функций в (6.5) обращает все уравнения в тождества.
Задача Коши для системы (6.5) состоит в нахождении решения системы, удовлетворяющего начальным условиям (где [math]y_,y_,\ldots,y_[/math] — известные числа):
В векторной форме задача Коши (6.5),(6.6) имеет вид
Теорема 6.2 (о существовании и единственности решения задачи Коши (6.5),(6.6)). Пусть выполнены следующие условия:
а) функции [math]f_i(x,y_1,\ldots,y_n),
i=\overline[/math] , определены и непрерывны в некоторой замкнутой области [math]\overline[/math] , а также имеют в [math]\overline[/math] ограниченные частные производные по переменным [math]y_1,\ldots,y_n[/math] ;
б) точка [math](x_0,y_,y_,\ldots,y_)[/math] лежит внутри области [math]\overline[/math] .
Тогда решение задачи Коши (6.5),(6.6) существует и единственно.
1. Во многих практических приложениях независимая переменная обозначается через [math]t[/math] и имеет смысл времени, поэтому задача Коши называется начальной задачей.
2. Понятия общего и частного решений, общего и частного интегралов для уравнения первого порядка и систем совпадают по форме, если заменить функцию [math]y(x)[/math] на вектор-функцию [math]Y(x),
f(x,y)[/math] на [math]F(x,Y)[/math] , а [math]y_0[/math] — на [math]Y_0[/math] .
Численные методы, рассматриваемые в данном разделе, пригодны для решения задач Коши, записанных в форме (6.5),(6.6). Чтобы решить задачу Коши (6.4) этими методами, ее необходимо привести к системе [math]n[/math] уравнений первого порядка, т.е. к виду (6.5),(6.6).
Получим точное решение модельного примера, используемого далее для демонстрации применения различных численных методов.
Пример 6.1. Найти аналитическое решение задачи Коши [math]Ty'+y=1,
y(0)=0[/math] , где [math]T>0[/math] — известное число, называемое постоянной времени.
Решение задачи Коши найдем с помощью известной методики.
1. Определим общее решение однородного уравнения [math]Ty'+y=0[/math] . Поскольку корень [math]\lambda=-1\!\!\not\,T[/math] соответствующего характеристического уравнения [math]T\lambda+1=0[/math] действительный, то [math]y_0(x)=Ce^= Ce^[/math] — общее решение однородного уравнения.
2. Частное решение неоднородного уравнения ищется в форме [math]y_=A[/math] , где [math]A=\text[/math] . После подстановки в решаемое уравнение получаем [math]y_=1[/math] .
3. Общее решение неоднородного уравнения получается как сумма общего решения однородного уравнения и частного решения неоднородного уравнения:
Ему соответствует семейство интегральных кривых, характеризующееся параметром [math]C[/math] , который может принимать произвольные значения.
4. Частное решение неоднородного уравнения находим из начального условия: [math]y(0)=C+1=0[/math] . Отсюда [math]C=-1[/math] и [math]y(x)=1-e^[/math] .
Пример 6.2. Записать дифференциальное уравнение второго порядка [math]2y''+y'+4y=6\sin,
y'(0)=2[/math] , в виде системы дифференциальных уравнений в нормальной форме Коши.
y_2(x)=y'_1(x)[/math] и запишем уравнение в форме [math]y''=-\fracy'-2y+3\sin[/math] , разрешенной относительно старшей производной. Тогда получим
Численные и приближенно-аналитические методы решения задачи Коши в отличие от аналитических методов позволяют найти искомую функцию [math]y(x)[/math] лишь приближенно. Но при этом численные методы являются более универсальными, так как с их помощью можно приближенно решать многие из задач, точные решения которых аналитическими методами не могут быть найдены. Аналитическими методами в основном решаются только линейные задачи и некоторые типы нелинейных задач, в то время как для численных методов эти ограничения отсутствуют.
Чтобы упростить изложение и в силу того, что численные методы легко обобщаются на системы уравнений, в дальнейшем будем рассматривать решение задачи Коши для уравнения первого порядка
Численное решение задачи ищется в узлах сетки [math]\Omega_n= \[/math] , где [math]h_=x_-x_,
i=\overline[/math] — расстояние между соседними узлами, называемое шагом интегрирования (параметром сетки). Если [math]h_=h=\text[/math] , сетка называется равномерной (регулярной), а если [math]h_=\text[/math] — неравномерной (нерегулярной). В случае равномерной сетки узлы находятся по формуле [math]x_= x_0+ih,
i=\overline[/math] , а в случае неравномерной (где [math]\delta_= \frac[/math] — параметр нерегулярности):
Решение находится в виде последовательности значений [math]\widehat_0,\widehat_1, \widehat_2,\ldots, \widehat_n[/math] , являющихся приближением значений [math]y_0,y(x_1),y(x_2),\ldots,y(x_n)[/math] точного решения [math]y(x)[/math] в узлах сетки [math]\Omega_[/math] (рис. 6.1).
Сеточное представление [math]y(x_i),
i=\overline[/math] , известной функции [math]y(x)[/math] (точного решения задачи Коши) называется проекцией [math]y(x)[/math] на сетку [math]\Omega_n[/math] .
Дискретные и непрерывно-дискретные методы
Численные методы решения обыкновенных дифференциальных уравнений делятся на две группы:
– дискретные методы , позволяющие найти решение только в узлах сетки. Эти методы называются еще разностными методами или методами сеток;
– непрерывно-дискретные методы , основанные на использовании дискретных методов и сплайн-функций для восполнения численных результатов. Они позволяют найти непрерывные решения дифференциальных уравнений.
Дискретные методы (методы сеток) подразделяются на явные и неявные. Значение [math]\widehat_[/math] , на (i+1)-м шаге может определяться явно:
где [math]\Phi(.)[/math] — некоторая функция, зависящая от конкретного метода (кроме последней рассчитанной точки [math](x_,\widehat_)[/math] могут использоваться еще [math](k-1)[/math] предыдущих точек), или неявно:
где искомая величина [math]\widehat_[/math] входит одновременно и в левую, и в правую часть.
Явные и неявные методы делятся также на одношаговые и многошаговые (k-шаговые). В одношаговых методах для расчета очередной точки [math](x_,\widehat_)[/math] требуется информация только о последней рассчитанной точке [math](x_,\widehat_)[/math] . В k-шаговых методах для нахождения точки [math]x_,\widehat_[/math] требуется информация о [math]k[/math] предыдущих точках (рис. 6.2).
Формулы (6.10),(6.11) в общем случае представляют собой нелинейные уравнения относительно [math]\widehat_[/math] и называются разностными схемами.
Численный алгоритм (метод) называется устойчивым, если численные результаты непрерывно зависят от входных данных и если погрешность остается ограниченной при заданных пределах изменения параметров численного алгоритма (шагов сетки, числа итераций и др.).
Сходимость приближенных методов является основной проблемой, от успешного преодоления которой зависит точность решения всей задачи. Численный алгоритм называется сходящимся, если при стремлении его параметров к определенным предельным значениям, например, при [math]h\to0[/math] (или при [math]s\to\infty[/math] , где [math]s[/math] — число итераций), результаты стремятся к точному решению.
При известном точном решении некоторой модельной задачи сходимость может быть проверена следующим образом. Фиксируется некоторая точка [math]x>x_0[/math] и строится последовательность сеток [math]\Omega_[/math] , таких, что [math]h\to0,[/math] [math]x=x_= x_+nh[/math] . Здесь для простоты считаем, что все сетки, образующие указанную последовательность, являются равномерными. Тогда, если [math]|\widehat_-y(x_)|\to0[/math] при [math]h\to0
(n\to\infty)[/math] , то метод является сходящимся в точке [math]x[/math] . Если метод сходится в каждой точке [math]x\in[c,d]\subset (a,b)[/math] , то он сходящийся на [math][c,d][/math] .
Локальная и глобальная ошибки
Локальной ошибкой численного метода на (i+1) -м шаге называется величина
где [math]y(x_)[/math] — значение точного решения при [math]x=x_[/math] , а [math]\widehat_[/math] — приближенное решение, получаемое по формуле (6.10) или (6.11) при условии, что вместо приближенных значений [math]\widehat_, \widehat_,\ldots, \widehat_[/math] используются значения, соответствующие точному решению, т.е. [math]y(x_), y(x_),\ldots, y(x_)[/math] .
Глобальной ошибкой называется величина [math]e_(h)= \widehat_-y(x_)[/math] , где [math]\widehat_[/math] — значение, получаемое по формулам (6.10) или (6.11) при [math]i=n-1[/math] .
Глобальная ошибка определяется:
а) ошибками округления и ошибками арифметических действий, обусловленными числом разрядов компьютера и характером выполняемых операций для расчета значения искомой функции в очередной точке [math]x_[/math] ;
б) методическими ошибками, определяемыми выбранным алгоритмом;
в) переходными ошибками, обусловленными тем, что при расчете значения р/+1 вместо точных значений [math]y(x_), y(x_),\ldots, y(x_)[/math] берутся приближенные значения [math]\widehat_, \widehat_,\ldots, \widehat_[/math] , полученные на предыдущих шагах.
Локальные ошибки "переносятся" в точку [math]x_n[/math] и формируют глобальную ошибку.
Число [math]p[/math] называется порядком (точностью) численного метода, если его глобальная ошибка есть [math]O[/math] большое от [math]h^p[/math] , то есть [math]e_n(h)= O(h^p)[/math] .
На практике в качестве характеристики точности метода часто используется величина [math]\varepsilon(h)= \max_\bigl|\widehat_-y(x_)\bigr|[/math] .
Рассмотрим введенные понятия более подробно на примере явных одношаговых методов, построенных для задачи (6.9). При этом формулу (6.10) представим в виде
где [math]\Psi(\widehat_,x_,h)[/math] — некоторая функция, определяемая конструкцией того или иного метода.
Обозначим [math]y(x,x_,\widehat_)[/math] — решение задачи Коши [math]u'=f(x,u),
u(x_)= \widehat_[/math] . Тогда локальная ошибка определяется выражением
Геометрическая интерпретация возникновения локальных и глобальной ошибок изображена на рис. 6.3.
Можно показать, что если локальная ошибка имеет порядок [math](p+1)[/math] , то есть [math]\varepsilon_(h)= O(h^)[/math] , то глобальная погрешность имеет на единицу меньший порядок, т.е. [math]e_(h)= O(h^p)[/math] .
Перейдем теперь к рассмотрению устойчивости численных методов. Она проверяется на "тестовом примере"
где [math]\mu[/math] — в общем случае комплексная константа. Дифференциальное уравнение в (6.12) является простейшим линейным уравнением, и для него можно получить значимые критерии устойчивости в явной форме.
Устойчивость методов решения задачи Коши
Метод называется устойчивым (ограниченно устойчивым), если существует такое число [math]h_>0[/math] , что при использовании метода для решения задачи (6.12), где [math]\operatorname\mu<0[/math] , с шагом [math]0<h<h_[/math] при [math]i\to\infty[/math] глобальная ошибка ограничена. Величина [math]h_[/math] называется критическим шагом. Если [math]h>h_[/math] , глобальная ошибка может неограниченно возрастать.
В ограниченно устойчивых методах при задании величины шага [math]h[/math] необходимо учитывать значение критического шага [math]h_[/math] . Для сложных дифференциальных уравнений и систем нахождение [math]h_[/math] является самостоятельной задачей, а свойство ограниченной устойчивости предупреждает вычислителя о возможных проблемах. Поэтому на практике становится актуальной задача конструирования таких методов, которые были бы устойчивы при любом значении шага, а его величина выбиралась бы только исходя из желаемой точности расчетов (при этом класс решаемых задач может быть ограничен).
Метод называется A-устойчивым , если при его применении с любым фиксированным положительным шагом [math]h[/math] все численные решения задачи (6.12) с комплексной константой [math]\mu
(\operatorname\mu<0)[/math] стремятся к нулю при [math]i\to\infty[/math] .
Область A-устойчивости — совокупность значений [math]h[/math] и [math]\mu[/math] , удовлетворяющих условию [math]\operatorname(h\mu)<0[/math] . Она изображена на рис. 6.4,а. Выполнение свойства A-устойчивости является желательным, поскольку если решение задачи (6.12) асимптотически устойчиво (в силу условия [math]\operatorname\mu<0[/math] корень характеристического уравнения находится в левой полуплоскости), то погрешность численного решения стремится к нулю при любой величине шага [math]h>0[/math] .
При исследовании устойчивости численного метода необходимо использовать соответствующую ему разностную схему для решения задачи (6.12) и привести ее к линейному разностному уравнению с постоянными коэффициентами:
Известно, что критерием устойчивости решения линейного разностного уравнения является требование расположения корней [math]\lambda_[/math] соответствующего характеристического уравнения [math]a_n \lambda^n+ a_\lambda^+ \ldots+a_0=0[/math] внутри круга единичного радиуса с центром в начале координат, т.е.
1. Численные методы, которые можно представить в виде (где | [math]|\alpha_0|+|\beta_0|\ne0,
f_= f(x_, \widehat_)[/math] называются линейными k-шаговыми методами)
\sigma(\xi)= \sum\limits_^ \beta_\xi^>[/math] . Линейный многошаговый метод является устойчивым, если для фиксированного значения [math]h\mu[/math] корни уравнения
лежат внутри круга единичного радиуса с центром в начале координат.
2. Для ограниченно устойчивых методов важной задачей является нахождение величины критического шага [math]h_[/math] . Если константа [math]\mu[/math] в уравнении (6.12) действительная, то можно найти интервал устойчивости .
3. Существуют определения, смягчающие свойство A-устойчивости. Приведем одно из них. Метод называется A(α)-устойчивым, [math]\alpha\in(0,\pi\!\!\not\,2)[/math] , если при его применении все численные решения уравнения (6.12) с фиксированным положительным шагом [math]h[/math] стремятся к нулю при [math]i\to\infty[/math] для всех [math]\mu[/math] , удовлетворяющих условию [math]|\arg(-\mu)|<\alpha,
|\mu|\ne0[/math] , где [math]\arg(-\mu)[/math] — аргумент комплексного числа [math](-\mu)[/math] . Область A(α)-устойчивости показана на рис. 6.4,5. Это условие применимо и для линейных систем с постоянными коэффициентами [math]y'=Ay[/math] , где [math]A[/math] — матрица коэффициентов, имеющая собственные значения [math]\lambda_,
i=\overline[/math] . Геометрическая интерпретация изображена на рис. 6.4,в.
4. Можно показать, что явные линейные многошаговые методы не могут быть A-устойчивыми.