Методы конечных разностей для численного решения обыкновенных дифференциальных уравнений и их систем. 


Мы поможем в написании ваших работ!



ЗНАЕТЕ ЛИ ВЫ?

Методы конечных разностей для численного решения обыкновенных дифференциальных уравнений и их систем.



 

Метод Эйлера имеет 1-й порядок точности с ошибкой на каждом шаге интегрирования порядка h2. Тогда, пренебрегая в уравнении (13) слагаемыми при h2 и выше, при замене исходного непрерывного отрезка [a,b] на совокупность узлов сетки {xi}, , с учетом определения yi¢ из заданного дифференциального уравнения получим:

 

. (14)

 

Модифицированный метод Эйлера является методом 2-го порядка точности с ошибкой порядка h3 и сохраняет в ряде Тейлора слагаемые при h и h2. При этом вторая производная аппроксимируется разностью y¢¢»Dy¢/Dx=(y¢i+1i)/h, тогда

 

,

 

где – уточненное значение функции в i+1 -м узле; – предварительное (прогнозное), которое определяется по методу Эйлера. Графическая интерпретация метода Эйлера представлена на рис. 2а, а его модификации на рис. 2б, где ai – угол наклона касательной к искомой функции y в точке xi (y¢i=tgai), aср =(ai+ai+1)/2, – точное значение функции в в i+1 -м узле, e М – погрешность метода:

 

y y

 

e М e М

Dy Dy

a ia ср

Dx=hDx=h

xixi+1 xxixi+1 x

а) б)

Рис. 2. Графическая интерпретация метода Эйлера и его модификации

 

Семейство методов РунгеКутты относится к методам 4-го порядка точности с ошибкой . Повышение точности обеспечивается за счет сохранения еще большего числа членов ряда Тейлора. Тогда в обобщенном виде для функции одной переменной можно записать [2, 15]:

 

(15)

где , , ,

,

aj, pj, bjn – фиксированные числа ().

При m=s=4 (s – порядок погрешности метода) наиболее употребительной для данного семейства является следующая совокупность расчетных формул:

 

, (16)

 

где ;

.

 

Для дальнейшего повышения эффективности метода может быть использовано автоматическое изменение шага. При этом после вычисления с шагом h все вычисления повторяются с шагом h/2. При выполнении условия | | <e переходят к определению с прежним шагом, в противном случае шаг уменьшают, если условие слишком сильное (<<e) – шаг увеличивают.

Многошаговые методы базируются на использовании значений функции и аргумента на нескольких предыдущих шагах (узлах). В сравнении с одношаговыми методами многошаговые алгоритмы могут быть более эффективными при той же точности или более точными при том же объеме вычислений. Однако при их реализации возникают и серьезные затруднения:

1) данные методы не обладают свойством самостартования, поэтому перед началом вычислений необходимо найти определенное количество значений искомой функции в начале отрезка интегрирования каким-либо одношаговым методом;

2) многошаговые методы (как они были первоначально описаны) могут применяться только с постоянным шагом h;

3) никакой многошаговый метод не может быть абсолютно устойчивым, если его порядок выше второго.

Ряд разработанных алгоритмов позволяет создавать численные модели, обходящие трудности (1) – (3) в том числе предложенный Нордсиком подход, получивший название многозначного метода, а также семейство жестко-устойчивых методов (формул) Гира [9].

В обобщенном виде, например для трехшагового метода, можно записать:

 

. (17)

 

Константы ai, bi подбираются априори. Если b0 =0, то метод будет явным, ненулевое b0 приводит к неявному методу. Аналогично описываются k -шаговые методы общего вида.

Метод Адамса–Бэшфорта является явным четырехшаговым методом 4-го порядка точности с ошибкой порядка h5 и определяется выражением

. (18)

 

Метод Адамса–Бэшфорта–Моултона является развитием предыдущего, реализующий так называемую схему прогноза и коррекции, в которой значение функции в очередной точке вначале определяется при помощи явной прогнозной формулы, а затем уточняется по неявной формуле коррекции:

 

(формула прогноза); (19)

 

(формула коррекции). (20)

 

Трехшаговый метод Гира эквивалентен одному из его четырехзначных методов и реализуется при помощи следующей неявной формулы:

 

. (21)

 

Точность и устойчивость численных методов связаны с обоснованным выбором шага h, при котором достигается снижение ошибки до величины, не превосходящей допустимой погрешности e. Параметр, который необходимо при этом контролировать называется глобальной ошибкой. Он определяется разностью между точным значением и значением, полученным при численном моделировании или . Эта разность в свою очередь имеет два источника.

1) Локальная ошибка (ошибка на шаге), или ошибка усечения, в точности совпадает с суммой остаточных членов в формуле Тейлора, влияние которых для конкретного метода было принято пренебрежимо малым. Например, для метода Эйлера это слагаемые при h2 и выше.

2) Распространяемая ошибка определяется множителем перехода (коэффициент усиления), который, например, для метода Эйлера может быть оценен как (1+hJ), где Jякобиан (матрица Якоби) [9].

Следовательно, глобальная ошибка метода для i+1 -го узла включает локальную ошибку для этого i+1 -го узла и переносимую с множителем перехода глобальную ошибку в i -м узле. Поэтому, если |1+hJ|<1, то ошибки в методе не возрастают и метод считается устойчивым. В противном случае метод неустойчив.

Неравенство |1+hJ|<1 эквивалентно неравенствам –2 < hJ<0, которые задают так называемый интервал устойчивости метода Эйлера. Если ОДУ неустойчиво, то эти неравенства никогда не выполняются и неустойчивость метода проявляется всегда. Если же ОДУ устойчиво (J<0), то эти два неравенства приводят к условию h<|2/J|.

Таким образом, ОДУ является устойчивым или неустойчивым вне зависимости от того, каким методом решается. Численный же метод интегрирования может быть устойчивым или неустойчивым для конкретной задачи (вида дифференциального уравнения, области интегрирования) в зависимости от размера шага.

На практике мы можем работать только с локальной ошибкой, причем желательно контролировать ее относительную величину, т.е. получать решение с определенным числом верных знаков, обычно от трех до пяти, вне зависимости от величины самого решения. Поэтому заменяем e на e|yi |. Тогда стратегия выбора шага сводится к следующему выражению:

 

, где . (22)

 

В целом погрешность рассмотренных методов e M может быть оценена с помощью правила Рунге, в соответствии с которым для заданной области интегрирования проводят вычисления с шагом h (i=1,N), а затем с шагом h/2 (i=1,2N1). Тогда

 

, (i=1,N). (23)

 

Краткая оценка методов. Одношаговые методы Рунге–Кутты просты, легко программируемы и для нежестких задач при обычных требованиях к точности (»10-5) не уступают многошаговым и многозначным методам. Однако эффективно решать жесткие задачи, которые чаще всего встречаются на практике, эти методы не способны. Еще одна проблема связана с выводом решения в случае автоматического выбора шага слишком большой длины. Программы многошаговых, или многозначных, методов не связаны подобными ограничениями и могут быть ориентированы как для решения жестких, так и нежестких задач [9].

 

Обыкновенным дифференциальным уравнением n -го порядка называется уравнение вида , где F - известная функция n +2 переменных, определенная в области ОÌR3; x – неизвестная переменная, заданная на отрезке [a,b]; y=y(x) – неизвестная функция; y¢,y¢¢,…,y(n) – ее производные, n – порядок уравнения.

Функция y(x) называется решением ДУ, если она при всех xÎ[a,b] удовлетворяет уравнению . Задача , или для xÎ[a,b], при y(xo)=yo, y¢(xo)=yo,1 , , y(n-1)(xo)=yo,n-1) называется задачей Коши или задачей с начальными условиями и имеет единственное решение.

Задача Коши для дифференциального уравнения n -го порядка при помощи вспомогательных функций (y1=y, y2=y1¢=y¢, …, yn=yn-1¢=y(n-1)) может быть сведена к задаче Коши для системы из n ДУ 1-го порядка. Данная система в векторной форме может быть записана в виде Y¢=F(x,Y), Y(x0)=Y0, где Y ¢= (y1¢(x), y2¢(x), …, yn¢(x)), F(x,Y) = (y2,y3,…,yn, f(x,y1,…,yn)).

Численное решение задачи Коши для этой системы состоит в построении таблицы приближенных значений координат yi1, yi2,…,yin вектора решения Yi(x) в точках x1, x2, …, xN (N – число узлов сетки). Эти координаты для каждой вспомогательной функции y1, y2, …,yn могут быть найдены при помощи методов, рассмотренных в лаб. раб. № 1.

Пусть необходимо решить ДУ второго порядка для xÎ[a,b] с начальными условиями y(a)=yo, y¢(a)=yo,1. Тогда, используя вспомогательные функции y1=y и y2=y1¢=y¢, с учетом y2¢=y¢¢, получим следующую систему из двух ДУ первого порядка:

 

(24)

 

При решении полученной системы методом Рунге–Кутты (см. лаб. раб. № 1) набор используемых формул примет следующий вид:

 

; (25)

;

где ;

;

;

;

;

;

;

.

 

Метод Рунге–Кутты–Фельберга (с автоматическим изменением шага) имеет погрешность порядка h5 и реализуется с учетом Y={y1,y2,… yj,…,ym}, K0={k01,…k0m}, K1={k11,…k1m}, K2={k21,…k2m}, K3={k31,…k3m}, K4={k41,…k4m}, K5={k51,…k5m} (m – число ДУ в системе, j – номер ДУ) при помощи следующих выражений:

 

;

;

;

;

;

;

 

; (26)

. (27)

 

Тогда погрешность

 

.

 

Если , то вычисления продолжаются с шагом h, если – шаг h уменьшается вдвое, если - шаг увеличивается в два раза.

 

Устойчивость уравнения и систем

Устойчивость уравнения определяется сближением кривых семейства решений сe-x уравнения y¢=f(x,y) с ростом х (якобиан J<0), если кривые расходятся, то уравнение называется неустойчивым (J>0). Уравнение общего вида может демонстрировать на различных участках отрезка интегрирования оба типа поведения. Например, уравнение y¢=2a(x1)y обладает семейством решений y=c×exp(a(x1)2). Кривые из этого семейства на отрезке [0,2] вначале будут удаляться друг от друга (якобиан J=¶f/¶y=2a(x1) при x<1 положителен), а затем сближаться (J=2a(x1) при x>1 отрицателен).

Жестким называют сверхустойчивое уравнение (J<<0).

Устойчивость систем ДУ связана с собственными значениями матрицы J. Собственные значения – это числа l1,…,ln, для которых матрица JliI будет вырожденной, т.е. определитель матрицы равен нулю (для одного уравнения l1=J).

 

В общем случае li – это комплексные числа, которые не постоянны и зависят от x и решения y. Доказано, что наибольшее влияние на устойчивость системы оказывает собственное значение с максимальной вещественной частью [9]. Для большинства практических задач наличие у всех собственных значений отрицательных вещественных частей обусловливает устойчивость решений. Положительные вещественные части обычно соответствуют областям с расходящимися неустойчивыми решениями. Системы, собственные значения которых имеют как положительные, так и отрицательные вещественные части, или все li нулевые также в большинстве случаев неустойчивы. Жесткость системы определяется наличием большого (по модулю) отрицательного собственного значения.

В ряде математических пакетов, таких как MathCAD, MATLAB и др. для решения обыкновенных дифференциальных уравнений (ОДУ или ODE), а также их систем эти методы реализованы в виде решателейОДУ [10].



Поделиться:


Последнее изменение этой страницы: 2017-02-07; просмотров: 370; Нарушение авторского права страницы; Мы поможем в написании вашей работы!

infopedia.su Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. Обратная связь - 3.90.255.22 (0.048 с.)