авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |   ...   | 9 |

«Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ ...»

-- [ Страница 5 ] --

Дополнительное требование L -устойчивости приводит к равенству p1 = a. (7.26) Это условие в случае (7.23) не выполняется. Таким образом, с исполь зованием (7.22) можно определить коэффициенты схемы (7.8), при ко торых она будет иметь порядок точности p, p = m + 1. Такой подход представляет интерес, когда требуется разработать алгоритм на основе одной схемы заданного порядка точности. Если говорить об алгоритме переменного порядка и шага, то построенные схемы могут оказаться малоэффективными. При создании алгоритма с автоматическим выбо ром численной схемы требуется наличие набора методов, в которых используется одна матрица Dn. В противном случае переход с одной схемы на другую будет приводить к существенному увеличению вы числительных затрат. Если параметры (7.8) определены из (7.22), то коэффициент a является функцией числа стадий m, что приводит при изменении m к перевычислению матрицы Dn. От этого ограничения можно избавиться, если не требовать максимального порядка точности.

A-устойчивые методы переменного порядка и шага Построим набор численных схем (7.8) с различным числом стадий m, m 1, порядка p = m и таких, чтобы параметр a в формуле вы числения Dn не зависел от числа стадий. Это будет выполнено, если 7.5. Методы решения линейных задач коэффициенты pi, 1 i m, определять из первого соотношения (7.22) при некотором значении параметра a. Тогда для всего набора формул может быть использована одна и та же матрица Dn. Параметр a естественно выбрать таким образом, чтобы весь набор численных формул обладал свойством A -устойчивости.

Возьмем неравенство для контроля точности такого набора чис ленных формул. При каждом m локальная ошибка n,m имеет вид (7.20) при p = m. Тогда для контроля точности можно использовать оценку ошибки вида (2.30). Для получения этой ошибки требуется оценка функции (t ) вида (t ) = f ( m1) f. Поступим следующим об, xm )T так, чтобы вы разом. Выберем компоненты вектора X m = ( x1, полнялось m xi ki = m hm f n(m1) f n + O(hm+1 ), (7.27) i = где m – постоянная. С использованием (7.13) нетрудно видеть, что вектор X m должен удовлетворять следующей системе линейных урав нений:

Bm1,m X m = 0, (7.28) где элементы Bm1,m определены в (7.9). Система (7.28) состоит из (m 1) -го уравнения относительно m неизвестных компонент вектора X m. Задаваясь произвольными значениями одной компоненты, ос тальные однозначно получим из (7.28). В дальнейшем при всех m бу дем полагать xm = 1. Тогда для определения xi, 1 i m 1, получим линейную систему Bm1,m1 X m 1 = U m с невырожденной матрицей Bm1,m1, где U m1 = (b1,m,..., bm1,m )T.

После определения X m значение m в (7.27) можно вычислить p p = a p 1 b pj x j, (7.29) j = Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ где p = m, а bmj определены в (7.9). Учитывая (7.27), оценку вида (2.30) можно вычислить по следующей формуле:

p n, p = p xi ki / p, (7.30) i = где p и p определяются формулами (7.21) и (7.29) соответственно.

Теперь для контроля точности p -стадийной схемы p -го порядка p xi ki | m / m |.

можно использовать, например, неравенство i = Обозначим через M максимальное число стадий в построенном набо ре численных формул. Тогда, учитывая, что при всех m, 1 m M, имеет место xm = 1, коэффициенты pi, 1 i m, схемы (7.8) и коэф фициенты xi, 1 i m 1, для оценки ошибок n,m, 1 m M, можно хранить в виде матрицы размерности M, в которой в m -й строке и столбце находятся компоненты векторов Pm и X m1 соответственно.

Постоянные | m / m | будем хранить в виде вектора размерности M.

L-устойчивые методы m-го порядка Теперь перейдем к построению методов вида (7.8), для которых выполняется условие Q ( x ) 0 при x. Из (7.24) следует, что это требование будет выполнено, если к условиям аппроксимации доба вить соотношение (7.26). Нетрудно видеть, что построить L устойчивую m -стадийную схему (7.8) выше m -го порядка точности, вообще говоря, нельзя, потому что система (7.22), (7.26) может быть несовместной.

Из сравнения (7.17) и (7.15) следует, что схема (7.8) будет иметь m -й порядок точности и для нее будет выполняться условие Q ( x ) при x, если ее параметры удовлетворяют системе (7.22) и урав нению (7.26). Выражая из линейной системы (7.22) вектор Pm и под ставляя его во второе соотношение (7.22), получим уравнение m -й степени относительно параметра a. Определяя один из корней полу ченного многочлена, например, методом деления отрезка пополам, из (7.22) получим коэффициенты pi, 1 i m, при которых схема (7.8) 7.5. Методы решения линейных задач будет иметь порядок точности p, p = m. Ниже для m, 1 i 10, при ведено по одному значению коэффициентов a (m), т. е.

a(1) = 1;

a(2) = 0, 2928932188;

a(3) = 0,1589839000;

a(4) = 0, 2204284103;

a(5) = 0, 2780538411;

a(6) = 0,3341423671;

a(7) = 0, 2040669280;

a(8) = 0, 2343731596;

a(9) = 0, 2643073554;

a(10) = 0, 2939936770, при которых (7.8) будет L -устойчивой схемой m -го порядка точности.

Оценку ошибки n,m для контроля точности вычислений снова можно вычислить по формуле (7.30) при p = m, где при определении m по формуле (7.21) используются вычисленные значения парамет ров a и pi, 1 i m. При выборе величины шага интегрирования можно использовать неравенство m D1 jn xi ki | m / m |, 1 jn 2.

n i = Данные методы, как и A -устойчивые схемы максимального поряд ка, представляют интерес, когда алгоритм интегрирования строится на основе одной численной формулы заданного порядка точности. Но при построении алгоритма переменного порядка и шага требуется набор схем, в которых используется одна и та же матрица Dn. Этого можно добиться, если при каждом m рассматривать методы ( m 1) -го поряд ка точности.

L-устойчивые методы переменного порядка и шага Из (7.24) следует, что схема (7.8) будет иметь (m 1) -й порядок точности и для нее будет выполнено Q( x) 0 при x, если ее коэффициенты удовлетворяют линейной системе алгебраических уравнений вида Bm1,m Pm = Vm1 (a), p1 = a, (7.31) где a – свободный параметр, a 0, а элементы матрицы Bm1,m и вектора Vm 1 определяются из (7.9) и (7.18) соответственно. Задаваясь Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ при всех m одним значением параметра a, из (7.31) получим кон кретные наборы коэффициентов pi, 1 i m, m 2, причем при вы числении стадий ki в (7.8) используем одну и ту же матрицу Dn при каждом значении параметра m. Параметр a, естественно, должен удовлетворять условию A -устойчивости.

Оценим величину ошибки для контроля точности вычислений. Ло кальная ошибка m -стадийной схемы имеет вид (7.20) при p = m 1, где при вычислении m1 по формуле (7.21) применяется решение (7.31). Тогда для контроля точности можно использовать ошибку n,m 1 вида (2.30). Для получения оценки n,m 1 поступим следующим образом. Определим компоненты вектора X m1 таким образом, чтобы выполнялось равенство m xi ki = m1hm1 f n(m2) f n + O(hm ), i = где m1 определяется по формуле (7.29) при p = m 1. С использова нием (7.13) видим, что вектор X m1 должен удовлетворять следующей системе уравнений:

Bm 2,m1 X m1 = 0. (7.32) Задаваясь произвольным значением одной компоненты вектора X m1, остальные однозначно определим из (7.32). Ниже будем полагать xm1 = 1 при всех m. Тогда для определения xi, 1 i m 2, получим линейную систему Bm 2,m2 X m2 = U m 2, где m 2, U m 2 = (b1 m1,..., bm2 m1 )T. Для контроля точности m -стадийной формулы (7.8), (7.31) можно использовать неравенство m1 D1 jn xi ki | m1 / m 1 |, 1 jn 2, n i =1 где Dn и ki, 1 i m, определены в (7.8), параметры m1 и m1 вы числяются при p = m 1 по формулам (7.21) и (7.29) соответственно.

Коэффициенты Pm, X m1 и | m1 / m1 | можно хранить по ана логии с параметрами A -устойчивых схем. Отметим, что все построен ные выше методы, порядок точности которых больше 2, имеют второй 7.6. Методы решения нелинейных задач. Схемы с одним вычислением правой части порядок точности, если их применять для решения нелинейной систе мы (7.2). Этот вывод следует из сравнения (3.31) и (7.15). Поэтому при разумном выборе величины шага интегрирования их можно применять для решения задачи Коши (7.2).

7.6. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ ЗАДАЧ СХЕМЫ С ОДНИМ ВЫЧИСЛЕНИЕМ ПРАВОЙ ЧАСТИ Как уже отмечалось выше, в случае решения задачи Коши для не линейной системы обыкновенных дифференциальных уравнений (7.2) нельзя построить метод вида (7.8) выше второго порядка точности.

Поэтому приведем методы первого и второго порядка.

Пусть m = 1. В этом случае (7.8) совпадает с методом типа Розен брока yn+1 = yn + p1k1, Dn k1 = hf ( yn ). (7.33) Используя разложение k1 в ряд Тейлора (7.13), получим yn+1 = yn + p1hf n + ap1h 2 f n f n + O(h3 ).

(7.34) Положим yn = y (tn ). Тогда из сравнения (3.31) и (7.34) получим, что при p1 = 1 и a = 1 / 2 формула (7.33) имеет максимальный поря док, равный 2. Применяя (7.33) для решения задачи y = y, y (0) = y0, Re() 0, будем иметь yn+1 = Q( x) yn, x = h, где Q ( x) = [1 + ( p1 a ) x ] / (1 ax). Из вида функции устойчивости Q ( x) имеем, что формула максимального порядка L -устойчивой не являет ся. Если отказаться от требования второго порядка, то схема (7.33) при p1 = a = 1 является L -устойчивой формулой первого порядка точности.

Пусть m = 2, т. е. рассмотрим схему вида yn+1 = yn + p1k1 + p2 k2, Dn k1 = hf ( yn ), Dn k2 = k1. (7.35) Подставляя представления k1 и k2 из (7.13) в первую формулу (7.35), получим yn+1 = yn + ( p1 + p2 )hf n + a ( p1 + 2 p2 ) h 2 f n f n + + a 2 ( p1 + 3 p2 )h3 f n 2 f n + O(h 4 ).

(7.36) Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ Полагая снова yn = y (tn ) и сравнивая (3.31) и (7.36) до членов с h включительно, получим условия второго порядка точности p1 + p2 = 1, a ( p1 + 2 p2 ) = 1 / 2. (7.37) Исследуем устойчивость (7.35). Применяя ее для решения задачи y = y, имеем yn+1 = Q( x) yn, x = h, где функция устойчивости Q ( x) записывается в виде Q( x) = 1 + ( p1 + p2 2a ) x + a( a p1 ) x 2 (1 ax) 2.

Тогда схема (7.35) будет L -устойчивой, если p1 = a. Подставляя это соотношение в (7.37), получим коэффициенты метода второго порядка точности, т. е.

p1 = a, p2 = 1 a, (7.38) где a определяется из условия L -устойчивости:

a 2 2a + 0,5 = 0. (7.39) Сравнивая (3.31) и (7.36) до членов с h3 включительно, получим, что локальная ошибка n,2 формулы (7.35) с коэффициентами (7.38) имеет вид n,2 = h3 (a 1 / 3) f 2 f + f 2 / 6 + O(h 4 ).

f (7.40) Уравнение (7.39) имеет два корня a1 = 1 2 / 2, a2 = 1 + 2 / 2. Выбе рем a = a1 – в этом случае коэффициент меньше в главном члене ло кальной ошибки (7.40).

Рассмотрим одновременно формулу типа Розенброка yn+1 = yn + ak1 + (1 a )k2, Dn k1 = hf ( yn ), Dn k2 = hf ( yn + ak1 ), где a = 1 0,5 2. В ней используются два вычисления правой части исходной задачи на шаге, она имеет второй порядок точности и L -устойчива. Построенная здесь формула (7.35) с коэффициентами (7.38) тоже обладает вторым порядком точности и L -устойчивостью, а их локальные ошибки (6.29) и (7.40) отличаются незначительно.

В то же время (7.35) требует на каждом шаге на одно вычисление функции f меньше, чем формула (6.12), при прочих равных затратах, что в некоторых случаях делает ее предпочтительнее.

7.7. Методы решения нелинейных задач. Схемы с двумя вычислениями правой части 7.7. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ ЗАДАЧ СХЕМЫ С ДВУМЯ ВЫЧИСЛЕНИЯМИ ПРАВОЙ ЧАСТИ Прежде чем перейти к исследованию (m, k ) -методов с двумя вы числениями правой части задачи, докажем следующее утверждение.

Теорема 7.3. Пусть в формуле (7.3) имеет место k = 2. Тогда при любом выборе множеств (7.1) и при любом значении параметра m нельзя построить метод (7.3) выше четвертого порядка точности.

Для простоты доказательство проведем для скалярной задачи = f ( y ), y (t0 ) = y0, t0 t tk, точное решение y (tn +1 ) которой мож y но записать в виде 1 f y (tn+1 ) = y (tn ) + hf + h 2 f + h3 f 2 f + f 2 6 + f f + h 4 f 3 f + 4 f 2 + f 3 24 + h5 f 4 f + 4 f 3 + ff ff + 5 f 2 f 2 + f 2 f 3 + f IV f 4 120 + O(h6 ), f (7.41) где элементарные дифференциалы вычислены на точном решении y (tn ).

С использованием (7.13) имеем, что при любом выборе множества M k второе вычисление функции f будет осуществляться в точке yn,c = yn + ci hi f n(i 1) f n + O( h5 ), i = где величина ci, 1 i 4, определяется через коэффициенты схемы (7.3). Поэтому для доказательства теоремы достаточно показать, что в представлении функции hf ( yn,c ) в виде ряда Тейлора не содержится слагаемого h5 f n 2 f n. В этом случае в разложении (7.41) имеется сла гаемое h5 f 2 f 3, а в соответствующих представлениях ki, 1 i m, оно отсутствует. Разлагая hf ( yn,c ), имеем hf ( yn,c ) = hf n + c1h 2 f n f n + h3 c2 f n 2 f n + c1 f n f n2 + + h 4 c3 f n3 f n + c1c2 f n f n f n2 + c1 f n f n / 6 + 3 Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ + h5 c4 f n 4 f n + c1c3 f n2 f n f n2 + c1 c2 f n f n f n / 2 + 2 + c1 f nIV f n4 / 24 + O(h6 ), что завершает доказательство. Теорема доказана.

Перейдем к исследованию (m, 2) -методов (7.3). Выберем множест ва (7.1) M 2 = {1, 2}, J i = {1}, 2 i m, т. е. рассмотрим схемы вида m yn +1 = yn + pi ki, Dn k1 = hf ( yn ), (7.42) i = Dn k2 = hf ( yn + 21k1 ) + 21k1, Dn ki = ki 1 + i1k1, 3 i m, где Dn определяется по второй формуле (7.3).

A-устойчивый метод третьего порядка Пусть m = 2. Подставим разложения k1 и k2 из (7.7) в первую формулу (7.42). Полагая yn = y (tn ) и сравнивая (3.31) и полученное представление приближенного решения до членов с h3 включительно, получим условия третьего порядка точности схемы (7.42), т. е.

p1 + (1 + 21 ) p2 = 1, ap1 + (a + 21 + 2a 21 ) p2 = 1 / 2, (7.43) ( ) a 2 p1 + a 2 + 2a21 + 3a 2 21 p2 = 1 / 6, 2 p2 = 1 / 3.

Исследуем совместность (7.43). Умножим первое уравнение (7.43) на 2a, а второе – на 1 и сложим. Затем умножим первое уравнение на 3a 2, а третье – на 1 и сложим. Будем иметь ap1 + ap2 21 p2 = 2a 1 / 2, 2a 2 p1 + 2a 2 p2 2a21 p2 = 3a 2 1 / 6.

Умножим на 2a первое из полученных соотношений и сложим со вторым:

a 2 a + 1 / 6 = 0. (7.44) Пусть 21 – свободный параметр. Тогда из последнего соотноше ния (7.43) выразим p2, а из первого и второго уравнений, которые 7.7. Методы решения нелинейных задач. Схемы с двумя вычислениями правой части линейны относительно p1 и 21, вычислим оставшиеся коэффи циенты. В результате определим p1 = 0,5[21 (12a21 321 + 2) 2a ], p2 =, (3 6a21 2) 21 = 21 21, (7.45) 2a где a вычисляется из равенства (7.44).

Исследуем устойчивость схемы (7.42) при m = 2. Применяя (7.42) для решения задачи y = y, получим yn+1 = Q( x) yn, x = h, где Q ( x) = 1 + (1 2a ) x + (a 2 2a + 0,5) x 2 / (1 ax)2. (7.46) При получении (7.46) использовались соотношения (7.43). Из (7.46) следует, что условие L -устойчивости схемы (7.42) имеет вид a 2 2a + 1 / 2 = 0, которое противоречит (7.44). Таким образом, при m = 2 схема (7.42) третьего порядка точности L -устойчивой не является.

Далее, уравнение (7.44) имеет два корня: a1 = (3 + 3) / 6 и a2 = ( 3 3 ) / 6. Так как при x имеем Q ( x) (a 2 2a + 0,5) / a 2, то нетрудно видеть, что при a = a1 схема (7.42) с коэффициентами (7.45) является A -устойчивой.

Свободный коэффициент 21 можно использовать для уничтожения одного слагаемого в главном члене локальной ошибки. Учитывая (7.7), видим, что при условии 3 p2 = 1 / 4 слагаемое h 4 f 3 в локальной f ошибке отсутствует. Используя это соотношение, из (7.43) будем иметь p1 = a 1 (76a 3) / 54, p2 = 16 / 27, 21 = 3 / 4, 21 = a 1 (3 54a ) / 32, где a определяется из (7.44).

L-устойчивый метод третьего порядка Исследуем схему (7.42) при m = 3. Подставим разложения ki, 1 i 3, из (7.7) в первую формулу (7.42). Полагая yn = y (tn ) и срав нивая (3.31) и полученное представление приближенного решения до Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ членов с h3 включительно, получим условия третьего порядка точно сти схемы (7.42):

p1 + (1 + 21 ) p2 + (1 + 21 + 31 ) p3 = 1, ap1 + (a + 21 + 2a 21 ) p2 + +(2a + 21 + 3a 21 + 2a31 ) p3 = 1 / 2, (7.47) ( ) a 2 p1 + a 2 + 2a21 + 3a 2 21 p2 + + ( 3a 2 + 3a21 + 6a 2 21 + 3a 2 31 ) p3 = 1 / 6, 2 ( p2 + p3 ) = 1 / 3.

y = y, получим условие Применяя (7.42) для решения задачи L -устойчивости a( p1 + p2 ) a 2 p221 = 0. (7.48) Функция устойчивости не приводится в силу ее громоздкости.

Исследуем совместность нелинейной системы (7.47), (7.48). Умно жим первое уравнение (7.47) на 6a 2, а второе – на 3a и сложим. За тем умножим первое уравнение (7.47) на 3a 2, а третье – на 1 и сло жим. Получим 3a 2 p1 + 3a 2 p2 3a21 p2 3a 2 21 p3 3a21 p3 = 6a 2 3a / 2, 2a 2 p1 + 2a 2 p2 2a21 p2 3a 2 21 p3 3a21 p3 = 3a 2 1 / 6.

Умножим второе из полученных равенств на 1 и сложим с пер вым. Учитывая (7.48), будем иметь a3 3a 2 + 3a / 2 1 / 6 = 0. (7.49) Пусть 21 и 21 – свободные параметры. Умножим первое урав нение (7.47) на 2a и сложим со вторым, учитывая (7.48). Из получен ного соотношения выразим p3. Затем последовательно определим из четвертого соотношения (7.47) параметр p2, из равенства (7.48) – 7.7. Методы решения нелинейных задач. Схемы с двумя вычислениями правой части коэффициент p1, а из первого условия (7.47) – значение 31. В ре зультате будем иметь p1 = (ac3 + c4 c5 ) / c3, ( ) ( 621c2 ), 2 p3 = 0,5c1 / c2, p2 = 2c2 321c 31 = {2c2 [ (1 a)c3 c4 + c5 2a (1 + 21 )c2 ]} (c2 c3 ), где параметр a определяется из условия L -устойчивости (7.49), c1 = 1 4a + 2a 2, c2 = 21 + a 21, c3 = 6a21c2, c4 = 2(21 a )c2, c5 = 321 (21 a)c1, а 21 и 21 – свободные коэффициенты.

Таким образом, при m = 3 можно построить L -устойчивую фор мулу (7.42) третьего порядка. Легко убедиться, что при любом значе нии m нельзя построить L -устойчивую формулу вида (7.42) выше третьего порядка точности. Для этого достаточно записать для произ вольного значения m условия аппроксимации четвертого порядка и потребовать L -устойчивости. Простейшие исследования полученной системы показывают, что она несовместна. Однако если не требовать L -устойчивости, то можно построить формулу вида (7.42) четвертого порядка.

A-устойчивый метод четвертого порядка Итак, пусть m = 4, т. е. рассмотрим схему вида yn+1 = yn + pi ki, Dn kn1 = hf ( yn ), i = Dn k2 = hf ( yn + 21k1 ) + 21k1, (7.50) Dn k3 = k2 + 31k1, Dn k4 = k3 + 41k1, где Dn определяется второй формулой (7.3). Подставим разложения ki, 1 i 4, из (7.7) в первую формулу (7.50). Полагая yn = y (tn ) и сравнивая (3.31) и полученное представление приближенного решения Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ до членов с h 4 включительно, получим условия третьего порядка точ ности схемы (7.50):

p1 + (1 + 21 ) p2 + (1 + 21 + 31 ) p3 + +(1 + 21 + 31 + 41 ) p4 = 1, ap1 + (a + 21 + 2a 21 ) p2 + (2a + 21 + 3a 21 + 2a31 ) p3 + + (3a + 21 + 4a 21 + 3a31 + 2a 41 ) p4 = 1 / 2, ( ) a 2 p1 + a 2 + 2a31 + 3a 2 21 p2 + ( ) + 3a 2 + 3a21 + 6a 2 21 + 3a 2 31 p3 + ( ) + 6a 2 + 4a21 + 10a 2 21 + 6a 2 31 + 3a 2 41 p4 = 1 / 6, (7.51) 2 ( p2 + p3 + p4 ) = 1 / 3, a2 ( p2 + p3 + p4 ) = 1 / 8, 21 a2 (0,5 p2 + p3 + 1,5 p4 ) = 1 / 24, 3 ( p2 + p3 + p4 ) = 1 / 4, 21 3 2 a p1 + a (321 + a + 4a 21 ) p2 + a (621 + 4a + 10a 21 + +4a31 ) p3 + a 2 (1021 + 10a + 2a 21 + 10a31 + 4a 41 ) p4 = 1 / 24.

Система (7.51) является нелинейной относительно pi и ij, ее иссле дование представляется достаточно трудоемкой задачей. Поэтому пе репишем (7.50) в таком виде, чтобы условия аппроксимации имели ли нейный вид относительно коэффициентов pi и ij. Рассмотрим вспомогательную схему yn +1 = yn + ci di, Dn d1 = hf ( yn ), (7.52) i = Dn d5 = hf ( yn + 21d1 ), Dn di = di 1, 2 i 4, 6 i 7.

Подставляя di, 1 i 7, в формулу (7.50), получим yn +1 = yn + p1d1 + ( 21 p2 + 31 p3 + 41 p4 )d 2 + + ( 21 p3 + 31 p4 )d3 + 21 p4 d 4 + p2 d5 + p3d6 + p4 d7.

7.7. Методы решения нелинейных задач. Схемы с двумя вычислениями правой части Отсюда видим, что коэффициенты (7.50) и (7.52) связаны соотноше p1 = c1, p2 = c5, p3 = c6, p4 = c7, 21 p4 = c4, ниями 21 p3 + 31 p4 = c3, 21 p2 + 31 p3 + 41 p4 = c2, последовательным исключением разрешая которые относительно pi и ij, будем иметь p1 = c1, p2 = c5, p3 = c6, p4 = c7, 21 = c4 / c7, 31 = (c3c7 c4 c6 ) / c7, (7.53) ( ) 2 2 41 = c2c7 c4 c5c7 c3c6 c7 + c4 c6 c7.

Теперь перейдем к исследованию численной схемы (7.52). Для это го требуются разложения стадий di, 1 i 7, в ряды Тейлора в окре стности точки yn до членов с h 4 включительно, которые легко полу чить из (7.7). Подставим их в первую формулу (7.52). Полагая yn = y (tn ) и сравнивая (3.31) и полученное представление приближен ного решения, получим условия четвертого порядка точности схемы (7.52), т. е.

ci = 1, a(c1 + 2c2 + 3c3 + 4c4 ) + (21 + a )c5 + i = +(2a + 21 )c6 + (3a + 21 )c7 = 1 / 2, ( ) a 2 (c1 + 3c2 + 6c3 + 10c4 ) + a 2 + a21 c5 + ( ) ( ) + 3a 2 + 3a21 c6 + 6a 2 + 4a21 c7 = 1 / 6, 2 (c5 + c6 + c7 ) = 1 / 3, ( ) a3 (c1 + 4c2 + 10c3 + 20c4 ) + a3 + 3a 221 c5 + (7.54) ( ) ( ) + 4a3 + 6a 221 c6 + 10a3 + 10a 221 c7 = 1 / 24, a2 (c5 + c6 + c7 ) = 1 / 6, a21 (0,5c5 + c6 + 1,5c7 ) = 1 / 24, 3 (c5 + c6 + c7 ) = 1 / 4.

Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ Исследуем совместность (7.54). Из четвертого и шестого уравнений получаем a = 3 / 8, из четвертого и восьмого имеем 21 = 3 / 4. Под ставляя полученные значения a и 21 в шестое и седьмое уравнения (7.54), имеем c5 + c6 + c7 = 16 / 27, c5 / 2 + c6 + 3c7 / 2 = 16 / 81. (7.55) Пусть c7 есть свободный параметр. Тогда подставляя соотношения (7.55) в первое, второе, третье и пятое уравнения (7.54), будем иметь 1 c1 1 11 / 1 c 20 / 1 2 3 4 2 = (7.56) 1 3 6 10 c3 80 / 81 c 4 10 20 c4 128 / 81 5c 1 В левой части системы (7.56) стоит матрица, составленная из элемен тов (7.9), т. е. она не вырожденная. Тогда из линейной системы (7.56) определим ci, 1 i 4, а из (7.55) – коэффициенты c5 и c6. В резуль тате имеем 21 = 3 / 4, c1 = (60 + 81c7 ) / 81, c2 = (18 324c7 ) / 81, c3 = (405c7 64) / 81, (7.57) c4 = (19 162c7 ) / 81, c5 = (64 + 81c7 ) / 81, c6 = (16 + 162c7 ) / 81, где c7 – свободный параметр. Подставляя (7.57) в (7.53), получим ко эффициенты схемы (7.50) четвертого порядка точности.

Отметим, что описанный способ исследования численной схемы (7.50) может быть применен для изучения других (m, k ) -методов, в частности описанных выше. В этом случае по части параметров ij условия аппроксимации линеаризуются, что значительно упрощает исследование.

L-устойчивый метод четвертого порядка Из приведенных выше рассуждений следует, что в классе (m, k) методов с двумя вычислениями правой части исходной системы диф ференциальных уравнений можно построить метод четвертого порядка 7.7. Методы решения нелинейных задач. Схемы с двумя вычислениями правой части точности. Однако он не обладает требуемыми свойствами устойчиво сти, что приводит к очевидным трудностям при решении жестких сис тем. Если иначе организовать вычисления, то можно построить L устойчивую численную схему четвертого порядка точности при k = 2.

Выберем множества M 2 = {1, 3}, J i = {2}, 3 i m, т. е. рассмотрим численные схемы вида m yn+1 = yn + pi ki, Dn k1 = hf ( yn ), Dn k2 = k1, i = Dn k3 = hf ( yn + 31k1 + 32 k2 ) + 32 k2, (7.58) Dn ki = ki 1 + i 2 k2, 4 i m, где матрица Dn определяется по второй формуле (7.3).

Пусть m = 4. Подставим разложения ki, 1 i 4, из (7.6) в первую формулу (7.58). Полагая yn = y (tn ) и сравнивая (3.31) и полученное представление приближенного решения до членов с h 4 включительно, получим условия четвертого порядка точности схемы (7.58), которые имеют вид p1 + p2 + (1 + 32 ) p3 + (1 + 32 + 42 ) p4 = 1, ap1 + 2ap2 + (a + 31 + 32 + 3a32 ) p3 + +(2a + 31 + 32 + 4a32 + 3a 42 ) p4 = 1 / 2, ( ) a 2 p1 + 3a 2 p2 + a 2 + 2a31 + 3a32 + 6a 2 42 p3 + ( ) + 3a 2 + 3a31 + 4a32 + 10a 2 32 + 6a 2 42 p4 = 1 / 6, ( ) a3 p1 + 4a3 p2 + a3 + 3a 231 + 6a 232 + 10a332 p3 + (7.59) ( ) + 4a3 + 6a 231 + 10a 232 + 2a332 + 10a3 42 p4 = 1 / 24, (31 + 32 ) 2 ( p3 + p4 ) = 1 / 3, a(31 + 32 )(31 + 232 )( p3 + p4 ) = 1 / 8, Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ a(31 + 32 )2 (0,5 p1 + p4 ) = 1 / 24, (31 + 32 )3 ( p3 + p4 ) = 1 / 4.

Система (7.59) является нелинейной относительно коэффициентов pi и ij, что вызывает определенные трудности с ее исследованием. По этому при изучении методов (7.58) поступим по аналогии с изучением численной формулы (7.50). Рассмотрим вспомогательную численную схему yn+1 = yn + ci di, Dn d1 = hf ( yn ), Dn d 2 = d1, i = Dn d3 = hf ( yn + 31d1 + 32 d 2 ), Dn d 4 = d 2, (7.60) Dn d5 = d3, Dn d6 = d 4.

Подставляя di, 1 i 6, из (7.60) в формулу (7.58), получим yn+1 = yn + p1d1 + p2 d 2 + p3d3 + +(32 p3 + 42 p4 )d 4 + p4 d5 + 32 p4 d6.

Сравнивая полученное соотношение с (7.60), видим, что коэффициен ты численных формул (7.58) и (7.60) связаны соотношениями p1 = c1, p2 = c2, p3 = c3, p4 = c5, (7.61) 32 = c6 / c5, 42 = (c4 c5 c3c6 ) / c5.

Теперь перейдем к изучению (7.60). Для этого требуются разложе ния стадий di, 1 i 6, в ряды Тейлора в окрестности точки yn до членов с h4 включительно, которые легко получить из (7.6). Подста вим их в первую формулу (7.60). Полагая yn = y (tn ) и сравнивая полу ченное представление приближенного решения и (3.31), получим ус ловия четвертого порядка точности схемы (7.60), т. е.

ci = 1, ac1 + 2ac2 + (a + 31 + 32 )c3 + 3ac4 + i = +(2a + 31 + 32 )c5 + 4ac6 = 1 / 2, 7.7. Методы решения нелинейных задач. Схемы с двумя вычислениями правой части ( ) a 2c1 + 3a 2 c2 + a 2 + 2a31 + 3a32 c3 + 6a 2 c4 + ( ) + 3a 2 + 3a31 + 4a32 c5 + 10a 2 c6 = 1 / 6, ( ) a3c1 + 4a3c2 + a3 + 3a31 + 6a32 c3 + 10a3c4 + ( ) + 4a3 + 6a 231 + 10a 232 c5 + 20a3c6 = 1 / 24, (7.62) a(31 + 32 )(31 + 232 )(c3 + c5 ) = 1 / 8, a(31 + 32 ) 2 (0,5c3 + c5 ) = 1 / 24, (31 + 32)3 (c3 + c5 ) = 1 / 4, (31 + 32 )2 (c3 + c5 ) = 1 / 3.

Исследуем совместность (7.62). Из четвертого и седьмого уравнений имеем 31 + 32 = 3 / 4, c3 + c5 = 16 / 27. (7.63) Тогда из пятого и шестого соотношений (7.62) получим 32 и c3 соот ветственно. Из (7.63) выразим 31 и c5. Учитывая первое уравнение (7.62), подставим полученные выражения для 31, 32, c3 и c5 во второе, в третье и восьмое равенства (7.62). Получим линейную систе му относительно параметров c2, c4 и c6, которая имеет вид ac2 + 2ac4 + 3ac6 = (22a + 5) / 54, 2a 2c2 + 7a 2 c4 + 16a 2c6 = (64a 2 + 29a 30) / 54, a 2 c2 + 3a 2 c4 + 6a 2 c6 = (32a 2 11a 6) / 54.

Разрешив данную систему, из первого уравнения (7.62) получим c1.

В конечном результате будем иметь c1 = a 2 (76a 2 29a + 3) / 27, c2 = a 2 (89a 146a 2 12) / 27, c3 = a 1 (32a 4) / 27, c4 = a 2 (72a 2 59a + 10) / 18, Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ c5 = a 1 (4 16a) / 27, c6 = a 2 (19a 18a 2 4) / 18, (7.64) 31 = a 1 (48a 9) / 32, 32 = a 1 (9 24a ) / 32.

y = y, получим условие L Применяя (7.60) для решения устойчивости a(c1 + c3 ) a 2 c331 = 0. (7.65) Выражение функции устойчивости не приводится в силу ее громоздко сти. Подставляя в (7.65) значение коэффициентов (7.64), будем иметь a 4 4a3 + 3a 2 2a / 3 + 1 / 24 = 0. (7.66) Подставляя теперь (7.64) в (7.61), получим коэффициенты схемы (7.58), т. е.

76a 2 29a + 3 146a 2 + 89a 12 32a p1 =, p2 =, p3 =, 2 2 27 a 27 a 27 a 4 16a 48a 9 9 24a p4 =, 31 =, 32 =, (7.67) 27 a 32a 32a 54a 2 + 57 a 12 864a3 + 828a 2 288a +, 42 =, 32 = (7.68) a (4 16a ) 8a 32a при которых она имеет четвертый порядок точности. Параметр a определяется из условия L -устойчивости (7.66). Заметим, что уравне ние (7.66) имеет два вещественных корня a1 = 0,5728160624821 и a2 = 3,100316735116.

Из приведенных выше рассуждений можно сделать следующие вы воды. Во-первых, свойство устойчивости (m, k ) -методов зависит от выбора множеств (7.1) или, что то же самое, от выбора способа реали зации численной схемы. Это следует из сравнения численных схем (7.50), (7.53), (7.57) и (7.58), (7.67), (7.68), которые фактически при одинаковых вычислительных затратах обладают различными свойст вами устойчивости.

Во-вторых, при k = 2 в классе (m, k ) -методов можно построить численную схему, которая по свойствам точности не уступает неявно му методу типа Рунге–Кутты с двумя вычислениями правой части 7.8. Замораживание матрицы Якоби в (3, 2)-методе решения жестких задач дифференциальной задачи, а при линейном анализе устойчивости она также не хуже. В то же время при записи (7.58) сразу заложен способ ее реализации, т. е. до начала расчетов можно оценить вычислитель ные затраты на шаг интегрирования.

Что касается неявных методов типа Рунге–Кутты, то для них вы числительные затраты в сильной степени зависят от способа реализа ции и, в частности, использование двухстадийной схемы вовсе не оз начает, что на каждом шаге будет дважды вычисляться правая часть исходной задачи. Поэтому на некоторых задачах (m, k ) -методы будут эффективнее неявных формул типа Рунге–Кутты.

7.8. ЗАМОРАЖИВАНИЕ МАТРИЦЫ ЯКОБИ В (3, 2)-МЕТОДЕ РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ Легко доказать, что нельзя построить (2, 2) -метод выше второго порядка с замораживанием матрицы Якоби. Нетрудно убедиться так же, что максимальный порядок точности (m, 2) -методов с заморажи ванием аналитической или численной матрицы Якоби не может быть выше третьего порядка. Ниже построена (3, 2)-схема, которую можно использовать с замораживанием как аналитической, так и численной матрицы Якоби.

Исследование (3, 2)-метода Для решения (7.2) рассмотрим численную формулу yn +1 = yn + p1k1 + p2 k2 + p3k3, Dn = E ahAn, (7.69) Dn k1 = hf ( yn ), Dn k2 = k1, Dn k3 = hf ( y n +1) + 32 k2.

Матрица An представима в виде An = f n + hBn + O(h 2 ), а схему y n +1 = yn + 31k1 + 32 k2 (7.70) называют внутренней или промежуточной схемой метода (7.69).

Для построения численной формулы третьего порядка разложим ki, 1 i 3, в ряды Тейлора в окрестности точки yn до членов с h включительно:

Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ k1 = hf n + ah 2 f n f n + a 2 h3 f n2 f n + ah3 Bn f n + O(h 4 ), k2 = hf n + 2ah 2 f n f n + 3a 2 h3 f n2 f n + 2ah3 Bn f n + O(h 4 ), k3 = (1 + 32 )hf n + (a + 31 + 32 + (7.71) ( ) +3a32 )h 2 f n f n + a 2 + 2a31 + 3a32 + 6a 2 32 h3 f n 2 f n + +(31 + 32 ) 2 h3 f n f n2 / 2 + (1 + 332 )ah3 Bn f n + O(h 4 ).

Подставляя (7.71) в первую формулу (7.69), запишем yn+1 = yn + [ p1 + p2 + (1 + 32 ) p3 ] hf n + + [ ap1 + 2ap2 + (a + 31 + 32 + +3a32 ) p3 ] h 2 f n f n + ( + a 2 p1 + 3a 2 p2 + a 2 + 2a31 + 3a32 + (7.72) ) + 6a 232 p3 h3 f n 2 f n + (31 + 32 )2 p3h3 f n f n2 / 2 + + [ ap1 + 2ap2 + (a + 3a32 ) p3 ] h3 Bn f n + O (h 4 ), где элементарные дифференциалы вычислены на приближенном ре шении yn.

Разложение точного решения y (tn +1 ) в ряд Тейлора в окрестности точки tn до членов с h3 включительно имеет вид ( ) y (tn+1 ) = y (tn ) + hf + h 2 f / 2 + h3 f 2 f + f 2 / 6 + O(h 4 ), (7.73) f f где элементарные дифференциалы вычислены на точном решении y (t n ).

Сравнивая (7.72) и (7.73) до членов с h3 включительно при усло вии yn = y (tn ), получим условия третьего порядка точности схемы (7.69), т. е.

p1 + p2 + (1 + 32 ) p3 = 1, ap1 + 2ap2 + (a + 31 + 32 + 3a32 ) p3 = 1 / 2, 7.8. Замораживание матрицы Якоби в (3, 2)-методе решения жестких задач ( ) a 2 p1 + 3a 2 p2 + a 2 + 2a31 + 3a32 + 6a 2 32 p3 = 1 / 6, (7.74) (31 + 32 ) 2 p3 = 1 / 3, ap1 + 2ap2 + (a + 3a32 ) p3 = 0.

Последнее соотношение (7.74) возникает за счет применения «испор ченной» матрицы Якоби. Исследуя совместность нелинейной алгеб раической системы (7.74), получим следующий набор коэффициентов:

p1 = 5 / 4 + 332 / 4, p2 = 1 332 / 2, p3 = 3 / 4, 32 = 2 / 3 31, (7.75) a 2 + 3a / 2 + 3a 2 32 / 4 3a31 / 4 = 1 / 6, где 31 и 32 – свободные коэффициенты.

Исследуем устойчивость схемы (7.69) на линейном скалярном уравнении y = y, y (0) = y0, (7.76) где – комплексное число, Re( ) 0. Применяя (7.69) для решения (7.76), получим yn+1 = Q3 ( x) yn, x = h, где функция устойчивости Q3 ( x) имеет вид Q3 ( x) = a3 + a 2 p1 + a 2 p3 a31 p3 x3 / (1 ax)3 + + 3a 2 2ap1 ap2 + (31 + 32 2a ) p3 x 2 / (1 ax)3 + + {[ 3a + p1 + p2 + (1 + 32 ) p3 ] x + 1} (1 ax)3.

Из вида Q3 ( x) следует, что для L -устойчивости (7.69) необходимо a 2 ap1 ap3 + 31 p3 = 0. (7.77) Запишем условие L -устойчивости промежуточной схемы (7.70).

Применяя (7.70) для решения (7.76), получим y n +1 = Q2 ( x) yn, где Q2 ( x) имеет вид { )} ( Q2 ( x) = 1 + (31 + 32 2a) x + a 2 a31 x 2 (1 ax) 2.

Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ Из вида Q2 ( x) следует, что промежуточная схема (7.70) будет L устойчивой, если 31 = a. С учетом соотношения (7.77) из (7.75) име ем, что коэффициенты метода (7.69), обладающего L -устойчивостью основной и промежуточной численных схем, записываются следую щим образом:

p1 = a, p2 = 2a + 3 / 2, p3 = 3 / 4, 32 = 2 / 3 a, 31 = a, (7.78) 32 = 4a / 3 5 / 3, a3 3a 2 + 3a / 2 1 / 6 = 0.

Схема (7.69) будет A -устойчива при 1 / 3 a 1,0685790. Послед нее уравнение (7.78) имеет три вещественных корня. Для практических расчетов рекомендуется корень a = a3 = 0, 43586652150846.

Контроль точности вычислений Для контроля точности вычислений и автоматического выбора ве личины шага интегрирования используем идею вложенных методов.

Для этого введем дополнительную стадию Dn k4 = k3. Теперь на осно ве уже вычисленных стадий k1, k2, k3 и k4 построим численную схе му второго порядка точности yn+1,2 = yn + b1k1 + b2 k2 + b3k3 + b4 k4. (7.79) Будем контролировать точность, основываясь на разности приближе ний, полученных по формулам (7.69) и (7.79), т. е.

|| yn+1 yn +1,2 ||, (7.80) где – требуемая точность расчетов, || || – некоторая норма в R N.

Для нахождения коэффициентов b1, b2, b3 и b4 запишем разло жение k4 в ряд Тейлора в окрестности точки yn до членов с h3 вклю чительно, т. е.

k4 = (1 + 32 )hf n + (2a + 31 + 32 + 4a32 ) h 2 f n f n + ( ) + 3a 2 + 3a31 + 4a32 + 10a 232 h3 f n2 f n + +(31 + 32 ) 2 h3 f n f n 2 / 2 + (2 + 432 )ah3 Bn f n + O(h 4 ).

7.8. Замораживание матрицы Якоби в (3, 2)-методе решения жестких задач Теперь, учитывая (7.71), условия второго порядка точности схемы (7.79) записываются следующим образом:

b1 + b2 + (1 + 32 )(b3 + b4 ) = 1, ab1 + 2ab2 + (a + 31 + 32 + 3a32 )b3 + (7.81) +(2a + 31 + 32 + 4a32 )b4 = 1 / 2.

Напомним, что схема (7.69) построена с учетом возможности замо раживания матрицы Якоби и ее численной аппроксимации. Ошибки, вносимые «испорченной» матрицей Якоби, устраняются за счет по следнего соотношения (7.74). Аналогичное требование к формуле (7.79) приводит к дополнительному уравнению на коэффициенты bi, 1 i 4, т. е.

ab1 + 2ab2 + (1 + 332 )ab3 + (2 + 432 ) ab4 = 0. (7.82) Решая совместно (7.81) и (7.82), получим b1 = 2a 1 / (4a / 3 2 / 3)b3, b2 = 2 3a + (4a / 3 2 / 3)b3, b4 = 3 / 4 b3, где b3 – свободный параметр. Положив b3 = 0, получим следующий набор ко эффициентов: b1 = 2a 1 / 2, b2 = 2 3a и b3 = 3 / 4. Применение нера венства для контроля точности (7.80) не приводит к дополнительным вычислениям правой части или обращениям матрицы Якоби. Введение дополнительной стадии k4, используемой только для контроля точно сти, увеличивает вычислительные затраты на один обратный ход мето да Гаусса на шаг интегрирования. Отметим одну важную особенность построенной оценки ошибки (7.80). Схема (7.69) L -устойчивая. Сле довательно, Q3 ( x) 0 при x. Так как для точного решения y (tn+1 ) = exp( x) y (tn ) задачи (7.76) выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки ошибки.

Однако для разности [ yn +1 yn +1,2 ] это свойство не выполняется, по тому что вспомогательная схема (7.79) L -устойчивой не является.

С целью исправления асимптотического поведения оценки ошибки вместо (7.80) будем контролировать следующее неравенство:

|| D1 jn ( yn +1 yn+1,2 ) ||, 1 jn 2. (7.83) Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ В этом случае поведение оценки ошибки при jn = 2 будет согласовано с поведением точного решения задачи (7.76) при x. Подчерк нем, что в смысле главного члена оценки ошибок, применяемые в не равенствах (7.80) и (7.83), совпадают при любом значении jn.

Использование неравенства (7.83) вместо (7.80) фактически не приводит к увеличению вычислительных затрат. Это связано с тем, что (7.83) при jn = 2 проверяется только в том случае, если оно нарушено при jn = 1. Такая ситуация встречается достаточно редко, в основном при быстром росте величины шага интегрирования. Однако это позво ляет выбирать шаг более правильно, тем самым уменьшается количе ство неоправданных повторных вычислений решения (возвратов).

Глава ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ Н иже исследованы (m, 3)-методы решения жестких задач. До казана теорема о максимальном порядке точности. Получены коэффициенты A -устойчивых и L -устойчивых численных схем пято го (максимального) порядка. При k = 3 рассматриваются (m, k ) методы следующего вида:

m yn+1 = yn + pi ki, Dn = E ahf n, i = Dn k1 = hf ( yn ), Dn ki = ki 1, 2 i r 1, r m, r Dn kr = hf yn + rj k j + r,r 1kr 1, j = Dn ki = ki 1 + i,r 1kr 1, r + 1 i s 1, (8.1) s Dn k s = hf yn + s, j k j + s,r 1kr 1 + s,s 1ks 1, j = Dn ki = ki 1 + i,r 1kr 1 + i,s 1ks 1, s + 1 i m.

В случае произвольных r, s и m описываются всевозможные пред ставления (m, 3) -методов.

8.1. ОБОЗНАЧЕНИЯ При исследовании (m, 3) -методов будут использоваться обозначе ния j r r 1 r rj i, rj, b2 = jrj, b3 = b1 = j =1 i = j =1 j = Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ r 1 s 1 s j rj 0,5i (i + 1), d1 = sj, d 2 = ( j r + 1)sj, b4 = j =1 i =1 j =r j =r j r +1 j r + s 1 s 1 sj i, d 4 = sj 0,5i (i + 1), d3 = j =r j =r i =1 i = r 1 s 1 s 1 sj + d1 + k,r 1 sj, c1 = j =1 k =r j =k s r 1 s k,r 1 j sj, c2 = jsj + d 2 + k =r j =k j = s 1 j r 1 j s sj i + d3 + k,r 1 sj i, c3 = k =r i = j =1 i =1 j =k r 1 j sj 0,5i (i + 1) + d 4 + c4 = j =1 i =1 s 1 s 1 j + k, r 1 sj 0,5i (i + 1), (8.2) j =k k =r i =1 r 1 s 1 s 1 m m p j + k,r 1 p j + k,r 1 p j + R1 = k =s j =1 k =r j =k j =k m s 1 m + l,r 1 k,s 1 p j, k =s l =r j =k r 1 s 1 s 1 k,r 1 ( j k + r ) p j + R2 = jp j + k =r j =1 j =k 8.1. Обозначения m m + k,r 1 ( j k + r ) p j + k =s j =k m s 1 m + l,r 1 k, s 1 ( j k + s ) p j, k =s l =r j =k s 1 j k + r r 1 s j p j i + k,r 1 p j R3 = i + k =r i = j =1 j =k i = m j k + r m + k,r 1 p j i + k =s i = j =k m j k + s m s + l,r 1 k,s 1 p j i, k =s j =k i =1 l =r s 1 j k + r r 1 s j p j 0,5i (i + 1) + k,r 1 p j 0,5i (i + 1) + R4 = j =1 i =1 k =r j =k i = m j k + r m + k,r 1 p j 0,5i (i + 1) + j =k k =s i = m j k + s m s + l,r 1 k,s 1 p j 0,5i (i + 1), k =s j =k l =r i = r 1 j p j i(i + 1)(i + 2) / 6 + R5 = j =1 i = s 1 j k + r s + k,r 1 p j i (i + 1)(i + 2) / 6 + j =k k =r i = Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ s 1 j k + r s + k,r 1 p j i (i + 1)(i + 2) / 6 + j =k k =r i = m j k + r m + k,r 1 p j i (i + 1)(i + 2) / 6 + j =k k =s i = m j k + s m s + l,r 1 k,s 1 p j i (i + 1)(i + 2) / 6, k =s j =k l =r i = s 1 m m k,s 1 p j, R6 = pj + j =r k =s j =k s 1 m m ( j r + 1) p j + k,s 1 ( j k + 2) p j, R7 = k =s j =r j =k m j k + j r + s 1 m p j i + k, s 1 p j i, R8 = i =1 k = s i = j =r j =k j r + s 1 pj R9 = 0,5i (i + 1) + j =r i = m j k + m + k, s 1 p j 0,5i (i + 1), j =k k =s i = j r + s 1 pj R10 = i (i + 1)(i + 2) / 6 + j =r i = m j k + m m + k, s 1 p j i (i + 1)(i + 2) / 6, R11 = p j, j =k k =s i =1 j =s j s + m m pj ( j s + 1) p j, R12 = R13 = i, j =s i = j =s 8.2. A-устойчивый (m, 3)-метод пятого порядка точности j s +1 m p j 0,5i(i + 1), R14 = j =s i = j s + m pj R15 = i (i + 1)(i + 2) / 6.

j =s i = 8.2. A-УСТОЙЧИВЫЙ (m, 3)-МЕТОД ПЯТОГО ПОРЯДКА ТОЧНОСТИ Построим схему пятого порядка точности. Заметим, при m = 4 и любом выборе множеств M 3 нельзя построить такую (m,3) -схему из за несовместности алгебраических систем, обеспечивающих методам пятый порядок точности. Поэтому пусть m = 5. Положим M 3 = {1, 4, 5}, J 4 = {3} и J 5 = {4}, т. е. рассмотрим схему вида yn+1 = yn + pi ki, Dn = E ahf n, i = Dn k1 = hf ( yn ), Dn k2 = k1, Dn k3 = k2, (8.3) Dn k4 = hf ( yn + 41k1 + 42 k2 + 43k3 ) + 43k3, Dn k5 = hf ( yn + 51k1 + 52 k2 + 53k3 + 54 k4 ) + 54 k4.

Запишем разложения в ряд Тейлора для стадий ki, 1 i 5, в окрест ности точки yn и подставим в первую формулу (8.3). Сравнивая при yn = y (tn ) с представлением точного решения y (tn +1 ) до членов с h включительно, получим условия пятого порядка точности схемы (8.3), а именно 1) R1 + R6 + R11 = 1;

2) a( R2 + R7 + R11 ) + b1R6 + c1R11 = 1 / 2;

3) a 2 ( R3 + R8 + R11 ) + a [b1R7 + b2 R6 + (c1 + c2 ) R11 ] + d1b1R11 = 1 / 6;

2 4) b1 R6 + c1 R11 = ;

Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ 5) a3 ( R4 + R9 + R11 ) + a 2 [b1R8 + b2 R7 + b3 R6 + (c1 + c2 + c3 ) R11 ] + + ad1 (2b1 + b2 ) R11 = 1 / 24;

( ) 2 2 6) a b1 R7 + c1 R11 + d1b1 R11 = 1 / 12;

7) a (b1b2 R6 + c1c2 R11 ) + d1b1c1R11 = 1 / 8;

3 8) b1 R6 + c1 R11 = ;

9) a 4 ( R5 + R10 + R11 ) + a3 [b1R9 + b2 R8 + b3 R7 + b4 R6 + + (c1 + c2 + c3 + c4 ) R11 ] + a 2 d1 (3b1 + 2b2 + b3 ) R11 = 1 / 120;

(8.4) 10) a 2 b1 R8 + c1 R11 + 2ad1b1 R11 = 1 / 60;

2 2 11) c1b1 d1R11 = 1 / 15;

( ) 3 3 12) a b1 R7 + c1 R11 + d1b1 R11 = 1 / 20;

13) a 2 (b1b2 R7 + c1c2 R11 ) + ad1b1 (c1 + b2 ) R11 = 1 / 40;

4 14) b1 R6 + c1 R11 = 1 / 5;

( ) 2 2 15) a b1 b2 R6 + c1 c2 R11 + d1b1c1 R11 = 1 / 10;

16) a 2b2 R6 + (ac2 + d1b1 ) 2 R11 = 1 / 20;

17) a 2b1b3 R6 + ac1 [ ac3 + d1 (b1 + b2 ) ] R11 = 1 / 30.

Исследуем совместность этой системы. Умножим 4-е уравнение (8.4) на b1 и сложим его с 8-м;

8-е уравнение (8.4) умножим на b1 и сложим с 14-м. Получаем 2 12c1 (c1 b1 ) R11 = (3 4b1 ), 20c1 (c1 b1 ) R11 = (4 5b1 ). (8.5) Умножая первое соотношение (8.5) на c1 и складывая со вторым, имеем c1 = 0,6(4 5b1 ) / (3 4b1 ). (8.6) 8.2. A-устойчивый (m, 3)-метод пятого порядка точности Из первого соотношения (8.5) получаем выражение для параметра R11, т. е.

R11 = (3 4b1 ) 12c1 (c1 b1 ).

(8.7) 2 R6 = (1 3c1 R11 ) / (3b1 ).

Из 4-го уравнения (8.4) определяем Из 11-го уравнения (8.4) находим ( ) d1 = 1 15c1b1 R11. (8.8) Умножая 6-е уравнение (8.4) на b1 и складывая с 12-м, запишем R11 = (3 5b1 ) 60ac1 (c1 b1 ).

(8.9) Сравнивая (8.7) и (8.9), имеем b1 = (0,6 3a) / (1 4a). (8.10) ( ) ( ) Из 6-го уравнения (8.4) определяем R7 = 1 12 d1b1 + ac1 R11 12ab1.

2 2 Умножая 7-е уравнение (8.4) на c1 и складывая с 15-м, получаем вы ражение для параметра b2 = (5c1 4) [ 40ab1 (c1 b1 ) R6 ]. Из 7-го уравне ния (8.4) определяем параметр c2 = [1 8(ab1b2 R6 + d1b1c1R11 ) ] / (8ac1R11 ).

Из 1-го уравнения (8.4) имеем R1 = 1 R6 R11. Из 2-го уравнения (8.4) R2 = [ 0,5 b1R6 c1R11 a ( R7 + R11 ) ] a. Из 3-го уравнения находим (8.4) с учетом выражения для R8 в (8.2), определяем параметр R3 :

{ R3 = 1 / 6 a 2 (2 R7 R6 + R11 ) a [b1R7 + b2 R6 + (c1 + c2 ) R11 ] d1b1R11} a 2. (8.11) Умножим 6-е уравнение (8.4) на 2a и сложим с 10-м. С учетом 4-го равенства (8.4) получим 20a 2 10a + 1 = 0. (8.12) Заметим, что R4 = 3R3 3R2 + R1 + 43 (3R7 2 R6 ). Тогда с учетом 1-, 2- и 3-го равенств (8.4), а также выражений для R8 и R9 в (8.2), 5-е уравнение (8.4) можно привести к виду Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ a3 43 (3R7 2 R6 ) + a 2 (b3 R6 + c3 R11 ) = R16, (8.13) где R16 определяется по формуле R16 = 1 / 24 a / 2 + 3a 2 / 2 a 3 a 2 [b1 (2 R6 R7 ) + +b2 ( R7 3R6 ) + (c1 2c2 ) R11 ] ad1 (b2 b1 ) R11.

Заметим, что имеет место соотношение R5 = 4 R4 6 R3 + 4 R2 R1 + + 43 (R7 R6 ). Тогда, с учетом 1-, 2-, 3- и 5-го уравнений (8.4), а так же выражений для b4, c4, R8, R9 и R10 в (8.2), 9-е равенство (8.4) приводится к виду a 2 (a 43 + b3 ) [ a( R7 R6 ) + d1R11 ] = R17, (8.14) где R17 = 1 / 120 a / 6 + a 2 2a3 + a a 2 (b1 2b2 ) [ a ( R7 R6 ) + d1R11 ].

Перепишем 17-е уравнение (8.4) следующим образом:

a 2 (b1b3 R6 + c1c3 R11 ) = R18, (8.15) где R18 = 1 / 30 ac1 (b1 + b2 )d1R11. Итак, получили линейную систему трех уравнений (8.13)–(8.15) относительно неизвестных 43, b3 и c3.

Решение этой системы имеет вид 43 = {c1R16 R18 R17 R6 (c1 b1 ) / [ a ( R7 R6 ) + {a3 [3c1(R7 R6 ) + b1R6 ]}, + d1R11 ]} b3 = R17 {a 2 [ a ( R7 R6 ) + d1R11 ]} a 43, ( ) ( a 2c1R11 ).

c3 = R18 a 2b1b3 R Подставляя найденные значения параметров в 13-е и 16-е уравнения (8.4), убеждаемся в том, что эти равенства также справедливы.

Вернемся к начальным обозначениям в (8.4). Зная R7 и R6, реша ем линейную систему уравнений относительно неизвестных p4 и 54.

8.2. A-устойчивый (m, 3)-метод пятого порядка точности Находим pi, 1 i 3, через известные значения Ri, 1 i 3. Решая линейную систему с известными параметрами bi, 1 i 3, находим коэффициенты 4 j, 1 j 3. Зная ci, 1 i 3, и значение 43, реша ем линейную систему уравнений относительно неизвестных 5 j, 1 j 3. В итоге получаем коэффициенты:

0,6(4 5b1 ) 3 4b 0,6(1 5a), c1 =, R 11= b1 =, 1 4a 3 4b1 12c1 (c1 b1 ) 2 2 1 / 3 c1 R11 1 / 12 (d1b1 + ac1 ) R, d1 =, R7 = R6 =, 2 2 b1 15c1b1 R11 ab 5c1, c2 = [1 / 8 (ab1b2 R6 + d1b1c1R11 ) ] (ac1R11 ), b2 = 40ab1 (c1 b1 ) R 0,5 b1R6 c1R11 a( R7 + R11 ) R1 = 1 R6 R11, R2 =, a { a 2 (2 R7 R6 + R11 ) a[b1R7 + b2 R6 + R3 = }a, +(c1 + c2 ) R11 ] d1b1R R16 = 1 / 24 0,5a + 1,5a 2 a3 a 2 [b1 (2 R6 R7 ) + +b2 ( R7 3R6 ) + (c1 2c2 ) R11 ] ad1 (b2 b1 ) R11, R17 = 1 / 120 a / 6 + a 2 2a3 + a 4 a 2 (b1 2b2 ) [ a ( R7 R6 ) + d1R11 ], R18 = 1 / 30 ac1 (b1 + b2 )d1R11, 43 = {c1R16 R18 R17 R6 (c1 b1 ) / [ a( R7 R6 ) + {a3 [3c1( R7 R6 ) + b1R6 ]}, + d1R11 ]} {a2 [ a( R7 R6 ) + d1R11 ]} a43, b3 = R Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ ( ) ( a 2c1R11 ), c3 = R18 a 2b1b3 R R R p4 = 2 R6 R7, 54 = 7, p5 = R11, R p3 = R3 2 R2 + R1 3 43 R7, 54 = d1, p2 = R2 R1 43 (2 R6 + R7 ) 2 p3, p1 = R1 43 R6 p3 p2, 43 = b3 2b2 + b1, 42 = 2b3 + 5b2 3b1, 41 = b3 3b2 + 3b1, 53 = c3 2c2 + c1 354 43, 52 = c2 c1 354 43 253, 51 = c1 54 43 54 53 52, где параметр a удовлетворяет уравнению (8.12). Применяя (8.3) к ли нейному скалярному уравнению y = y с комплексным, Re() 0, получаем yn+1 = Q ( z ) yn, z = h, Q( z ) = {1 z (5a 1) + z 2 (20a 2 10a + 1) / 2 z 3 (60a 60a 2 + 15a 1) / 6 + z 4 (120a 4 240a3 + 120a 2 20a + 1) / z 5 a5 a 4 ( p1 + p4 + p5 ) + a3 ( 41 p4 (51 + 54 ) p5 ) } a 24154 p5 (1 az )5.

Отсюда следует необходимое условие L -устойчивости схемы (8.3), т. е.

a5 a 4 ( p1 + p4 + p5 ) + + a3 [41 p4 + (51 + 54 ) p5 ] a 24154 p5 = 0.

Это соотношение можно записать в виде a5 5a 4 + 5a3 5a 2 / 3 + 5a / 24 1 / 120 = 0. (8.16) 8.3. L-устойчивый (m, 3)-метод пятого порядка точности Данное условие не совпадает с (8.12), т. е. схема (8.3) не является L-устойчивой. Уравнение (8.12) имеет два вещественных корня a1 = 1 / 4 + 5 / 20 и a2 = 1 / 4 5 / 20. Ниже приведены коэффициенты при a = a2 :

a1 = 0,36180339887499;

p1 = 0,67495669696944;

p2 = 0,19157093400385;

p3 = 0, 46186417743554;

p4 = 0, 40323911084859;

p5 = 0,61380586637395;

41 = 0,67335567072963;

42 = 0, 28140414756465;

43 = 0,13065037833070;

43 = 1, 26110782817440;

51 = 1,11273790471360;

52 = 0,91147229526022;

53 = 0, 47464957019196;

54 = 0,14445594674080;

54 = 0,54170980452313.

Итак, схема (8.3) с коэффициентами (8.16) при a = 0, 25 5 / является A -устойчивой схемой пятого порядка. Построить L устойчивый (5, 3)-метод пятого порядка точности нельзя. Если запи сать при m = 5 и любом выборе множества M 3 условия аппроксима ции пятого порядка и потребовать L-устойчивость, то полученная сис тема окажется несовместной.

8.3. L-УСТОЙЧИВЫЙ (m, 3)-МЕТОД ПЯТОГО ПОРЯДКА ТОЧНОСТИ Перейдем к построению L -устойчивой схемы пятого порядка точ ности. Пусть m = 6. Положим M 3 = {1, 4, 6}, J 4 = {3}, т. е. рассмотрим схему вида yn+1 = yn + pi ki, Dn = E ahf n, i = Dn k1 = hf ( yn ), Dn k2 = k1, Dn k3 = k2, Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ Dn k4 = hf ( yn + 41k1 + 42 k2 + 43k3 ) + 43k3, Dn k5 = k4, (8.17) Dn k6 = hf ( yn + 61k1 + 62 k2 + 63k3 + 64 k4 + 65 k5 ).

Разложения в ряд Tейлора для ki, 1 i 4, в окрестности точки yn совпадают с их разложениями в (8.3). С учетом обозначений (8.2) вы пишем разложения в ряд Tейлора в окрестности точки yn для стадий k5 и k6. Подставляя полученное разложения в первую формулу (8.17) и сравнивая при yn = y (tn ) с представлением точного решения y (tn +1 ) до членов с h5 включительно, получаем условия пятого порядка точ ности схемы (8.17), т. е.

1) R1 + R6 + R11 = 1;

2) a( R2 + R7 + R11 ) + b1R6 + c1R11 = 1 / 2;

3) a 2 ( R3 + R8 + R11 ) + a [b1R7 + b2 R6 + (c1 + c2 ) R11 + d1b1 ] R11 = 1 / 6;

2 4) b1 R6 + c1 R11 = ;

5) a3 ( R4 + R9 + R11 ) + a 2 [b1R8 + b2 R7 + b3 R6 + +(c1 + c2 + c3 ) R11 ] + a [ d1 (b1 + b2 ) + d 2b1 ] R11 = 1 / 24;

( ) 2 2 6) a b1 R7 + c1 R11 + d1b1 R11 = 1 / 12;

7) a (b1b2 R6 + c1c2 R11 ) + d1b1c1R11 = 1 / 8;

3 8) b1 R6 + c1 R11 = ;

(8.18) 9) a 4 ( R5 + R10 + R11 ) + a3 [b1R9 + b2 R8 + b3 R7 + b4 R6 + + (c1 + c2 + c3 + c4 ) R11 ] + a 2 [ ( d1 + d 2 + d3 )b1 + +(d1 + d 2 )b2 + d1b3 R11 ] = 1 / 120;

( ) 10) a 2 b1 R8 + c1 R11 + a (d1 + d 2 )b1 R11 = 1 / 60;

2 2 11) c1b1 d1R11 = 1 / 15;

( ) 3 3 12) a b1 R7 + c1 R11 + d1b1 R11 = 1 / 20;

8.3. L-устойчивый (m, 3)-метод пятого порядка точности 13) a 2 (b1b2 R7 + c1c2 R11 ) + ad1b1 (c1 + b2 ) R11 = ;

4 14) b1 R6 + c1 R11 = 1 / 5, ( ) 2 2 15) a b1 b2 R6 + c1 c2 R11 + d1b1c1 R11 = 1 / 10;

16) a 2b2 R6 + (ac2 + d1b1 ) 2 R11 = 1 / 20;

17) a 2b1b3 R6 + ac1 (ac3 + d1b2 + d 2b1 ) R11 = 1 / 30.

Исследуем совместность этой системы. Поступая по аналогии с тем, как это было проделано для системы (8.4), получаем выражения для параметров c1, R11, R6, d1, b1, R7, b2, c2, R1, R2 и R3, кото рые совпадают с их выражениями в (8.6)–(8.8), (8.10)–(8.11) соответст венно. Умножим 6-е уравнение (8.18) на a и сложим с 10-м. С уче том уравнений для параметра R8 в (8.2) получаем выражение для ( 60ab12 R11 ).


d 2 = 1 5a 60a 2b1 ( R7 R6 ) Как и ранее, подставляя найденные значения параметров в 13-е и 16-е уравнения (8.18), убеж даемся в том, что эти равенства также справедливы. Принимая во вни мание выражения для параметров R8 и R9 в (8.2), 5-е уравнение (8.18) приведем к виду a3 43 (3R7 2 R6 ) + a 2 (b3 R6 + c3 R11 ) = R16, (8.19) где R16 = 1 / 24 a / 2 + 3a 2 / 2 a3 a 2 [b1 (2 R6 R7 ) + +b2 ( R7 3R6 ) + (c1 2c2 ) R11 ] a (d1b2 2d1b1 + d 2b1 ) R11.

Девятое уравнение (8.18), с учетом 1, 2, 3 и 5-го соотношений (8.18), а также выражений для b4, c4, d3, R10, R9 и R8 в (8.2), будет выглядеть следующим образом:

a3 43 [ a ( R7 R6 ) + (3d 2 2d1 ) R11 + (c1 2c2 ) R11 ] + (8.19a) + a 2b3 [ a ( R7 R6 ) + d1R11 ] = R17, Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ где R17 = 1 / 120 a / 6 + a 2 2a3 + a 4 a3 (b1 2b2 )( R7 R6 ) a 2 [b1 (2d1 d 2 ) + b2 (d 2 3d1 ) ] R11.

Перепишем 17-е уравнение (8.18) так:

a 2 (b1b3 R6 + c1c3 R11 ) = R18, (8.20) где R18 = 1 / 30 ac1 (d1b2 + d 2b1 ) R11. Получили линейную систему (8.19), (8.19a) и (8.20) относительно неизвестных 43, b3 и c3. Решая эту систему, получим 43 = {(c1R16 R18 ) [ a( R7 R6 ) + d1R11 ] { R6 R17 (c1 b1 )} a3 [ c1 (3R7 2 R6 ) ( a ( R7 R6 ) + d1R11 ) } R6 (c1 b1 )(a ( R7 R6 ) + ( 3d 2 2d1 ) R11 ), { b3 = R17 a3 43 [ a ( R7 R6 ) + {a2 [ a( R7 R6 ) + d1R11 ]}, + (3d 2 2d1 ) R11 ]} ( ) ( a 2c1R11 ).

c3 = R18 a 2b1b3 R Вернемся к начальным обозначениям в (8.2). Зная R7 и R6, реша ем линейную систему двух уравнений относительно неизвестных p4 и p5. Находим pi, 1 i 3, через известные параметры Ri, 1 i 3, и 43. Решая линейную систему с известными параметрами d1 и d 2, находим параметры 64 и 65. Зная ci, 1 i 3, решаем линейную систему уравнений относительно неизвестных 6 j, 1 j 3. В итоге получаем следующие коэффициенты:

3(4 5b1 ) 3 4b 0,6(1 5a), c1 =, R11 = b1 =, 1 4a 5(3 4b1 ) 12c1 (c1 b1 ) 8.3. L-устойчивый (m, 3)-метод пятого порядка точности ( ) 2 2 1 / 12 d1b1 + ac1 R 1 / 3 c1 R11 R6 =, d1 =, R7 =, 2 2 b1 15c1b1 R11 ab 5c1, c2 = [1 / 8 (ab1b2 R6 + d1b1c1R11 ) ] (ac1R11 ), b2 = 40ab1 (c1 b1 ) R R1 = 1 R6 R11, R2 = [1 / 2 b1R6 c1R11 a ( R7 + R11 )] a, { R3 = 1 / 6 a 2 (2 R7 R6 + R11 ) a [b1R7 + b2 R6 + (c1 + c2 ) R11 ] d1b1R11} a 2, ( 60ab12 R11 ), d 2 = 1 5a 60a 2b1 ( R7 R6 ) R16 = 1 / 24 a / 2 + 3a 2 / 2 a3 a 2 [b1 (2 R6 R7 ) + b2 ( R 3R6 ) + (c1 2c2 ) R11 ] a(d1b2 2d1b1 + d 2b1 ) R11, R17 = 1 / 120 a / 6 + a 2 2a3 + a 4 a3 (b 2b2 )( R7 R6 ) a 2 [b1 (2d1 d 2 ) + b2 (d 2 3d1 )] R11, R18 = 1 / 30 ac1 (d1b2 + d 2b1 ) R11, (8.21) 43 = {(c1R16 R18 ) [ a ( R7 R6 ) + d1R11 ] R6 R17 (c {a3 c1(3R7 2R6 ) ( a(R7 R6 ) + d1R11 ) b1 )} } R6 (c1 b1 ) ( a( R7 R6 ) + (3d 2 2d1 ) R11 ), { b3 = R17 a3 43 [ a ( R7 R6 ) + {a2 [a( R7 R6 ) + d1R11 ]}, +(3d 2 2d1 ) R11 ]} ( a 2c1R11 ), c3 = R18 a 2b1b3 R6 p5 = R7 R6, Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ p4 = 2 R6 R7, p6 = R11, p3 = R3 2 R2 + R1 3 43 R7, p2 = R2 R1 43 (2 R6 + R7 ) 2 p3, p1 = R1 43 R6 p3 p2, 65 = d 2 d1, 64 = 2d1 d 2, 63 = c3 2c2 + c1 3d 2 43, 62 = c2 c1 43 (2d1 + d 2 ) d 2 + d1 263, 61 = c1 43d1 d1 63 62, 43 = b3 2b2 + b1, 42 = 2b3 + 5b2 3b1, 41 = b3 3b2 + 3b1.

Применяя (8.17) к задаче y = y с комплексным, Re() 0, по лучаем yn+1 = Q( z ) yn, z = h, Q( z ) = {1 z (6a 1) + z 2 (30a 2 12a + 1) / z 3 (120a3 90a 2 + 18a 1) + + z 4 (360a 4 480a3 180a 2 24a + 1) / 24 z 5 (720a5 1800a 4 + +180a 2 24a + 1) / 24 z 5 (720a5 1800a 4 + +1200a3 300a 2 + 30a 1) / 120 + z 6 [a a5 ( p1 + p4 + p6 ) + a 4 ( 41 p4 + (61 + 64 ) p6 ) } a34164 p6 (1 az )6.

Необходимое условие L -устойчивости схемы (8.17) имеет вид a 6 a5 ( p1 + p4 + p6 ) + + a 4 [41 p4 + (61 + 64 ) p6 ] a34164 p6 = 0.

8.4. Теорема о максимальном порядке точности (m, 3)-методов С учетом проделанных выше выкладок это уравнение перепишется в виде 15 4 10 3 5 2 1 a 6 6a5 + a a + a a+ = 0. (8.22) 2 3 8 20 Данное уравнение имеет четыре вещественных корня. Схема (8.17), (8.21) есть L -устойчивая схема пятого порядка точности, если a ко рень (8.22).

8.4. ТЕОРЕМА О МАКСИМАЛЬНОМ ПОРЯДКЕ ТОЧНОСТИ (m, 3)-МЕТОДОВ Теорема 8.1. Пусть в методах (7.3) имеет место k = 3. Тогда при любом выборе множества M 3 и любом значении параметра m нельзя построить (m, k ) -метод выше пятого порядка точности.

Для доказательства рассмотрим (m, k ) -методы вида (8.1). Разлагая ki, 1 i m, в ряды Тейлора до членов с h6 включительно и подстав ляя их в первую формулу (8.1), получаем следующий вид приближен ного решения:

( ) yn+1 = yn + c11hf n + c21h 2 f n f n + h3 c31 f n 2 f n + c32 f n f n2 + ( ) + h 4 c41 f n3 f n + c42 f n f n f n2 + c43 f n f n f n2 + c44 f n f n + ( + h5 c51 f n4 f n + c52 f n2 f n f n2 + c53 f n 2 f n + c54 f n f n f n + 3 +c55 f n f n f n f n2 + c56 f nIV f n4 + c57 f n f n f n + c58 f n f n f n f n f n + ) ( +c59 f n f n 2 f n2 + h6 c61 f n5 f n + c62 f n3 f n f n2 + c63 f n f n 2 f n + + c64 f n2 f n f n + c65 f n 2 f n f n f n2 + c66 f n f n f n4 + c67 f n f nIV f n4 + 3 +c68 f n f n f n f n + c69 f n f n + c610 f n 2 f n f n + c611 f n 2 f n2 f n f n + 3 V5 3 +c612 f n f n f n 2 f n2 + c613 f n f n f n f n f n f n + c614 f nIV f n f n4 + Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ +c615 f n f n f n4 + c616 f n f n 2 f n + c617 f n f n f n f n f n2 + c618 f n f n3 f n2 + ) +c619 f n f n 2 f n f n f n + c620 f n f n f n f n + O(h7 ), где cij зависят от коэффициентов схемы (8.1). Сравнивая представле ние точного решения y (tn +1 ) в виде ряда Тейлора с рядом для при ближенного решения yn +1, видим, что выполнение уравнений 1 1 1 1 1 c32 =, c44 =, c53 =, c56 =, c66 =, c615 = (8.23) 6 24 30 120 144 является необходимым условием шестого порядка точности схемы (8.1). Нетрудно видеть, что с учетом (8.2) имеет место ( ) 12 1) c32 = b1 R6 + c1 R11 ;

( ) 13 2) c44 = b1 R6 + c1 R11 ;

( ) 14 3) c56 = b1 R6 + c1 R11 ;

(8.24) 4) c53 = b1 c1d1R11;

5) c66 = b1 c1d1R11, 1 6) c615 = b1 c1 d1R11.

Доказательство теоремы следует из несовместности уравнений (8.24) при условии (8.23). Теорема доказана.

8.5. ТЕОРЕМА О МАКСИМАЛЬНОМ ПОРЯДКЕ ТОЧНОСТИ (m, 3)-МЕТОДОВ С ЗАМОРАЖИВАНИЕМ МАТРИЦЫ ЯКОБИ В рамках (m, k ) -методов значительно проще решается проблема замораживания матрицы Якоби, а также некоторых других видов ап проксимаций. На основе исследования (m, k ) -методов с одним, двумя и тремя вычислениями правой части сделано предположение, что мак 8.5. Теорема о максимальном порядке точности (m, 3)-методов симальный порядок точности (m, k ) -методов с замораживанием со ставляет (k + 1), однако строгие доказательства отсутствуют. Ниже показано, что максимальный порядок точности (m, 3) -методов с замо раживанием и/или численным вычислением матрицы Якоби равен 4.

Ниже будут рассматриваться численные методы решения задачи Коши для автономной системы обыкновенных дифференциальных уравнений вида y = f ( y ), y (t0 ) = y0, t0 t tk. (8.25) Произвольный (m, 3) -метод можно записать в виде m yn+1 = yn + pi ki, Dn k1 = hf ( yn ), Dn ki = ki 1, 1 i s1, i = ( ) s Dn ks1 = hf yn + i =1 s1,i ki + s1,s1 1ks1 1, Dn ki = ki 1 + i,s1 1ks1 1, s1 i s2, (8.26) s2 Dn k s2 = hf yn + s2,i ki + s2, s1 1k s1 1 + s2,s2 1ks2 1, i = Dn ki = ki 1 + i,s1 1ks11 + i,s2 1ks2 1, i s2.

В случае замораживания или численной аппроксимации матрицы Якоби предполагается, что при реализации методов используется мат рица Dn, представимая в виде Dn = E ahAn, An = f n + hBn,1 + h 2 Bn,2 + h3 Bn,3 + O(h 4 ), (8.27) где Bn,i – некоторые матрицы, не зависящие от размера шага интегри рования.

Разложение матрицы Dn 1 в ряд Тейлора в окрестности точки yn до членов с h 4 включительно записывается следующим образом:

Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ ( ) Dn 1 = E + ahf n + ah 2 Bn,1 + af n 2 + ah3 ( Bn,2 + af n Bn,1 + aBn,1 f n + ( ) + a 2 f n3 + ah 4 Bn,3 + af n Bn,2 + aBn,2 f n + aBn,1 + a 2 f n Bn,1 f n + ) + a 2 f n 2 Bn,1 + a 2 Bn,1 f n 2 + a3 f n4 + O(h5 ).

Прежде чем сформулировать теорему о максимальном порядке точно сти, докажем следующие леммы.

Лемма 8.1. Стадии (m, 3) -метода (8.26) представимы в виде ли нейной комбинации стадий следующих (m, k ) -схем:

m yn+1 = yn + piI kiI, Dn k1 = hf ( yn ), Dn kiI = kiI1, i 1, I (8.28) i = m yn+1 = yn + piII kiII, kiII = kiI, 1 i s1, i = (8.29) s1 Dn ks = hf yn + s,i kiII II II, Dn kiII = kiII1, i s1, 1 i = m yn+1 = yn + piIII kiIII, kiIII = kiII, 1 i s2, i = s2 1 s2 Dn ks = hf yn + s,,iI kiI + s,,iII kiII III III III, (8.30) 2 2 i =1 i = s Dn kiIII = kiIII, i s2.

При доказательстве легко заметить, что первые ( s1 1) стадии ме тодов (8.26) и (8.28) совпадают, т. е.

ki = kiI, 1 i s1. (8.31) Математической индукцией покажем, что ki, s1 i s2, представимы в виде i j,s11ksI1+i j.

ki = kiII + (8.32) j = s 8.5. Теорема о максимальном порядке точности (m, 3)-методов II II Положив s, j = s1, j, получим, что правая часть в стадиях ks и ks 1 вычисляется в одной и той же точке. Отсюда, c учетом (8.31), запишем II I ks1 = ks + s1,s1 1ks. Легко видеть, что полученное выражение пред 1 ставимо в виде (8.32). Предположив, что соотношение (8.32) выполне но для некоторой стадии ki, s1 i s2, найдем выражение для стадии I ki +1. Из (8.26) следует, что Dn ki +1 = ki + i +1,s1 1ks 1. Подставляя сю да выражение (8.32), получим I i ki +1 = Dn 1 kiII + j,s1 1ks +i j + i +1,s1 1ks 1.

I j = s После несложных преобразований запишем i + j,s11Dn 1ksI1+i j.

ki +1 = Dn 1kiII + j = s Полученное выражение с учетом (8.28) и (8.29) представимо в виде (8.32).

Аналогично покажем, что стадии ki, i s2, представимы в виде s2 II i n, s2 1 ks +i n + j, s1 1kiI+ s + s n j + ki = kiIII + 2 n = s j = s I + n,s1 1ks +i n. (8.33) Для этого подберем s,,iI и iIII, II таким образом, чтобы правая часть III III в стадиях ks2 и ks вычислялась в одной и той же точке. Приравнивая III соответствующие выражения для стадий ks2 и ks и воспользовав шись формулой (8.32), имеем Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ i s,,iI = s2,i, 1 i s1, s,,iI = s2,i j,s11, III III 2 j = s s1 i s2, s,,iII = s2,i, s1 i s2.


III В результате можно записать s2 1 s2 k s2 = Dn 1hf yn + s,,iI kiI + s,,iII kiII III III + 2 i =1 i = s + Dn 1 s2,s1 1k s1 1 + Dn 1 s2,s2 1k s2 1.

Отсюда, с помощью (8.28)–(8.30) и (8.32) получим s2 j,s11ksI1+ s2 j.

III II I k s2 = k s + s2, s2 1k s + s2,s1 1ks + s2,s2 2 2 j = s Нетрудно видеть, что данное выражение представимо в виде (8.33).

Предположим теперь, что формула (8.33) верна для некоторого i, i s2. Подставляя в соответствующую строку (8.26) выражения (8.32), (8.33) и воспользовавшись формулами (8.28)–(8.30), запишем s2 II i n, s2 1 ks +(i +1)n + j,s11k s + s + (i +1)n j + ki +1 = kiIII + I + 2 n = s j = s II s2 1 + n,s1 1k s +(i +1)n + i +1, s1 1k s + i +1,s2 1 k s + j, s1 1k s + s j.

I I I 2 j =s 1 1 Полученное выражение представимо в виде (8.33), т. е. формула (8.33) верна. Из формул (8.31)–(8.33) следует, что стадии (8.26) представимы в виде линейной комбинации стадий (8.28)–(8.30). Лемма доказана.

Лемма 8.2. Разложения стадий (m, k ) -методов вида (8.28) с од ним вычислением правой части при условии (8.27) представимы в виде kiI = hf n + i,1ah 2 f n f n + i,2 a 2 h3 f n 2 f n + i,1ah3 Bn,1 f n + + i,1ah 4 Bn,2 f n + i,2 a 2 h 4 f n Bn,1 f n + i,2 a 2 h 4 Bn,1 f n f n + 8.5. Теорема о максимальном порядке точности (m, 3)-методов + i,3a3h 4 f n3 f n + i,1ah5 Bn,3 f n + i,2 ah5 f n Bn,2 f n + (8.34) + i,2 ah5 Bn,2 f n f n + i,2 a 2 h5 Bn f n + i,3a3h5 f n Bn,1 f n f n + + i,3a3h5 f n 2 Bn,1 f n + i,3a3h5 Bn,1 f n 2 f n + i,4 a 4 h5 f n 4 f n + O(h6 ), где i, j – некоторые постоянные, а элементарные дифференциалы вычислены на приближенном решении yn.

Доказательство леммы проведем с помощью метода математиче I ской индукции. Разлагая стадию k1 в ряд Тейлора до членов с h5, запишем k1 = hf n + ah 2 f n f n + a 2 h3 f n 2 f n + ah3 Bn,1 f n + ah 4 Bn,2 f n + I + a 2 h 4 f n Bn,1 f n + a 2 h 4 Bn,1 f n f n + a3h 4 f n3 f n + ah5 Bn,3 f n + + a 2 h5 f n Bn,2 f n + a 2 h5 Bn,2 f n f n + a 2 h5 Bn f n + a3h5 f n Bn,1 f n f n + + a3h5 f n 2 Bn,1 f n + a3h5 Bn,1 f n2 + a 4 h5 f n 4 f n + O( h6 ), I т. е. разложение стадии k1 представимо в виде (8.34). Пусть теперь разложение i -й стадии имеет вид (8.34). Применяя слева Dn 1 к (8.34), запишем kiI+1 = hf n + ( i,1 + 1)ah 2 f n f n + ( i,2 + i,1 + 1)a 2 h3 f n2 f n + +( i,1 + 1)ah3 Bn,1 f n + ( i,1 + 1)ah 4 Bn,2 f n + +( i,2 + i,1 + 1)a 2 h 4 f n Bn,1 f n + ( i,2 + i,1 + 1)a 2 h 4 Bn,1 f n f n + +( i,3 + i,2 + i,1 + 1)a3h 4 f n3 f n + ( i,1 + 1)ah5 Bn,3 f n + +( i,2 + i,1 + 1)a 2 h5 f n Bn,2 f n + ( i,2 + i,1 + 1)a 2 h5 Bn,2 f n f n + +( i,2 + i,1 + 1)a 2 h5 Bn f n + ( i,3 + i,2 + i,1 + 1)a3h5 f n Bn,1 f n f n + 2 Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ +( i,3 + i,2 + i,1 + 1)a3h5 f n2 Bn,1 f n + ( i,3 + i,2 + i,1 + +1)a3h5 Bn,1 f n 2 f n + ( i,4 + i,3 + i,2 + i,1 + 1)a 4 h5 f n 4 f n + O(h6 ).

Полученное выражение представимо в виде (8.34), при этом коэффи циенты i, j задаются рекуррентными формулами:

i +11 = i,1 + 1, i +1,2 = i,2 + i,1 + 1,, i +1,3 = i,3 + i,2 + i,1 + 1, i +1,4 = i,4 + i,3 + i,2 + i,1 + 1, (8.35) 11 = 1,2 = 1,3 = i,4 = 1.

, Формулы (8.35) можно переписать в следующем виде:

i +1, j = i, j + i, j 1, 1 j 4, 1, j = 1, 1 j 4, i,0 = 1 i.

Здесь i,0 имеет смысл коэффициента при элементарном дифферен циале hf n в i -й стадии. Лемма доказана.

Лемма 8.3. В (m, 2) -методах (8.29) при условии (8.27) разложения стадий kiII, i s1, имеют вид kiII = hf n + (i,1 + b0 ) h 2 f n f n + ( ) + b0 h3 f n f n2 + i,2 + (i + 1 s1 )ab0 h3 f n 2 f n + +i,1h3 Bn,1 f n + i,3 + a 2b0 (i + 2 s1 )(i + 1 s1 ) h 4 f n3 f n + + ab0b1h 4 f n f n f n2 + ab0 (i + 1 s1 )h 4 f n f n f n2 + + b0 h 4 f n f n + i,2 h 4 f n Bn,1 f n + i,5 h 4 Bn,1 f n f n + 3 +i,1h 4 Bn,2 f n + i,4 h5 f n 4 f n + a 2b1 h5 f n f n f n f n f n + 8.5. Теорема о максимальном порядке точности (m, 3)-методов + a 2b0b2 h5 f n f n2 f n2 + ab0 b1h5 f n f n f n + 1 + a (i + 1 s1 )b0 h5 f n f n f n + b0 h5 f nIV f n4 + 3 3 (8.36) 6 + a 2 (i + 1 s1 )b0b1h5 f n f n f n f n2 + 1 + a 2 i + 1 s1 + (i + 1 s1 )(i s1 ) b0 h5 f n2 f n f n2 + 2 +i,2 h5 f n Bn,2 f n + i,3h5 f n 2 Bn,1 f n + i,6 h5 f n Bn,1 f n f n + 1 + ab0b1h5 f n Bn,1 f n2 + ab0b1h5 f n f n Bn,1 f n + 2 +i,7 h5 Bn,1 f n2 f n + (i + 1 s1 )ab0 h5 Bn,1 f n f n2 + +i,5 h5 Bn,2 f n f n + i,1h5 Bn,3 f n + i,8 h5 Bn,1 f n + O (h6 ), где коэффициенты i, j заданы в лемме 8.2, i, j – некоторые посто янные, а bi определяются следующим образом:

s1 1 s1 1 s1 1 s1 s1, j, s1, j j,1, s1, j j,2, b3 = s1, j j,3.

II II II II b0 = b1 = b2 = j =1 j =1 j =1 j = Доказательство проведем с помощью математической индукции.

Пользуясь тем, что разложения стадий kiI, 1 i s1, известны из лем мы 8.2 и задаются формулами (8.34) с рекуррентными коэффициента II ми (8.35), найдем разложение стадии k s в ряд Тейлора в окрестности точки yn до членов с h5 включительно. Разложение функции s1 hf yn + s, j k I II j j = Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ в ряд Тейлора имеет вид s1 1 s1 1 II hf yn + s, j k j = hf n + hf n s, j k j + II j =1 1 j = 2 s1 1 II s1 1 II 1 + hf n s, j k j + h f n s, j k j + j =1 1 j =1 1 2 s1 1 II s, j k j + O(h6 ).

+ hf nIV (8.37) j =1 1 Раскрыв скобки в (8.37) с учетом леммы 8.2 и выписав из получивше гося выражения члены с h5 включительно, получим s hf yn + s +1, j k I II 2 3 = hf n + b0 h f n f n + ab1h f n f n + j j = + b0 h3 f n f n2 + a 2b2 h 4 f n3 f n + ab1h 4 f n Bn,1 f n + ab0b1h 4 f n f n f n2 + + b0 h 4 f n f n + ab1h5 f n Bn,2 f n + a 2b2 h5 f n 2 Bn,1 f n + 3 + a 2b2 h5 f n Bn,1 f n f n + ab0b1ah5 f n Bn,1 f n2 + 1 + ab0b1ah5 f n f n Bn,1 f n + a3b3h5 f n 4 f n + a 2b1 h5 f n f n f n f n f n + 2 + a 2b0b2 h5 f n f n2 f n2 + ab0 b1h5 f n f n f n + + b0 h5 f nIV f n4 + O(h6 ). (8.38) Применяя слева Dn 1 к (8.38), запишем k s = hf n + (a + b0 )h 2 f n f n + b0 h3 f n f n2 + II ( ) + ab1 + ab0 + a 2 h3 f n 2 f n + ah3 Bn,1 f n + 8.5. Теорема о максимальном порядке точности (m, 3)-методов ( ) + a 2b2 + a 2b1 + a 2b0 + a3 h 4 f n3 f n + ab0b1h 4 f n f n f n2 + ( ) 12 + ab0 h 4 f n f n f n2 + b0 h 4 f n f n + ab1 + a 2 h 4 f n Bn,1 f n + 2 ( ) ( + ab0 + a 2 h 4 Bn,1 f n f n + ah 4 Bn,2 f n + a3b3 + ) + a3b2 + a3b1 + a 3b0 + a 4 h5 f n4 f n + a 2b1 h5 f n f n f n f n f n + 12 + a 2b0b2 h5 f n f n2 f n2 + ab0 b1h5 f n f n f n + ab0 h5 f n f n f n + 3 2 1 4 5 IV b0 h f n f n + a 2b0b1h5 f n f n f n f n2 + + ( ) + a 2b0 h5 f n2 f n f n2 + ab1 + a 2 h5 f n Bn,2 f n + 2 ( ) ( ) + a 2b2 + a 2b1 + a3 h5 f n 2 Bn,1 f n + a 2b2 + a 2b0 + a3 h5 f n Bn,1 f n f n + 1 + ab0b1h5 f n Bn,1 f n2 + ab0b1h5 f n f n Bn,1 f n + 2 ( ) + a 2b1 + a 2b0 h5 Bn,1 f n 2 f n + ab0 h5 Bn,1 f n f n2 + ( ) + ab0 + a 2 h5 Bn,2 f n f n + ah5 Bn,3 f n + a 2 h5 Bn,1 f n + O(h6 ).

II Легко заметить, что стадия ks представима в виде (8.36), причем s1,1 = a, s1,2 = ab1 + a 2, s1,3 = a 2b2 + a 2b1 + a 3, s1,4 = a3b3 + a3b2 + a3b1 + a3b0 + a 4, s1,5 = ab0 + a 2, s1,6 = a 2b2 + a 2b0 + a3, s1,7 = a 2b1 + a 2b0 + a 3, s1,8 = a 2.

Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ Предположим теперь, что стадия kiII, i s1, представима в виде (8.36).

Умножая (8.36) слева на Dn 1, получим kiII1 = hf n + (i,1 + b0 + a )h 2 f n f n + b0 h3 f n f n2 / 2 + + + i,2 + (i + 1 s1 )ab0 + a(i,1 + b0 ) + a 2 h3 f n 2 f n + +(i,1 + a )h3 Bn,1 f n + i,3 + (i + 2 s1 )(i + 1 s1 )a 2b0 + ai,2 + +(i + 1 s1 )a 2b0 + a 2 (i,1 + b0 ) + a3 h 4 f n3 f n + ab0b1h 4 f n f n f n2 + 1 + a(i + 1 s1 )b0 + ab0 h 4 f n f n f n2 + b0 h 4 f n f n + 2 2 ( ) + i,2 + ai,1 + a 2 h 4 f n Bn,1 f n + i,5 + a(i,1 + b0 ) + a 2 h 4 Bn,1 f n f n + + (i,1 + a )h 4 Bn,2 f n + i,4 + ai,3 + (i + 1 s1 )(i + 2 s1 )a3b0 + + a 2 i,2 + (i + 1 s1 ) a3b0 + a3 (i,1 + b0 ) + a 4 h5 f n 4 f n + 1 + a 2b1 h5 f n f n f n f n f n + a 2b0b2 h5 f n f n f n2 + ab0 b1h5 f n f n f n + 2 2 1 + (i + 1 s1 )ab0 + ab0 h5 f n f n f n + b0 h5 f nIV f n4 + 3 3 6 + a 2 (i + 1 s1 )b0b1 + a 2b0b1 h5 f n f n f n f n2 + 2 1 + a 2b0 i + 1 s1 + (i + 1 s1 )(i s1 ) + a 2 (i + 1 s1 )b0 + a 2b 2 ( ) ( ) h5 f n2 f n f n2 + i,2 + ai,1 + a 2 h5 f n Bn,2 f n + i,3 + ai,2 + a 2 i,1 + a h5 f n2 Bn,1 f n + i,6 + ai,5 + a 2 (i,1 + b0 ) + a3 h5 f n Bn,1 f n f n + 8.5. Теорема о максимальном порядке точности (m, 3)-методов 1 + ab0b1h5 f n Bn,1 f n2 + ab0b1h5 f n f n Bn,1 f n + 2 + i,7 + ai,2 + (i + 1 s1 )a 2b0 + a 2 (i,1 + b0 ) + a h5 Bn,1 f n 2 f n + a(i + 1 s1 )b0 + ab0 h5 Bn,1 f n f n2 + 2 + i,5 + a(i,1 + b0 ) + a 2 h5 Bn,2 f n f n + (i,1 + a) h5 Bn,3 f n + ( ) + i,8 + ai,1 + a 2 h5 Bn,1 f n + O (h6 ).

Полученное выражение представимо в виде (8.36), при этом коэффи циенты i, j определяются по следующим рекуррентным формулам:

i +11 = a + i,1, i +1,2 = i,2 + ai,1 + a 2,, i +1,3 = i,3 + ai,2 + a 2i,1 + a3, i +1,4 = i,4 + ai,3 + a 2 i,2 + a3i,1 + a 4 + a3 (i + 3 s1 )(i + 2 s1 )b0 / 2, i +1,5 = i,5 + ai,1 + ab0 + a 2, i +1,6 = i,6 + ai,5 + a 2 i,1 + a3 + a 2b0, i +1,7 = i,7 + ai,2 + a 2 i,1 + a3 + (i + 2 s1 )a 2b0, i +1,8 = i,8 + ai,1 + a 2.

Лемма доказана.

Лемма 8.4. В методах (8.30) при условии (8.27) разложения стадий kiIII, i s2, имеют вид kiIII = hf n + ( i,1 + c0 )h 2 f n f n + c0 h3 f n f n2 + ( ) + i,2 + (i + 1 s2 ) ac0 + b0 c0 + c1 h3 f n 2 f n + i,1h3 Bn,1 f n + II +i,3h 4 f n3 f n + ac0 (i + 1 s2 ) / 2 + b0 c0 h 4 f n f n f n2 + 2 II Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ ( ) II +(i,2 + c1 ) h 4 f n Bn,1 f n + c0 c1 + b0c0 h 4 f n f n f n2 + +c0 h 4 f n f n / 6 + i,2 + (i + 1 s2 )ac0 h 4 Bn,1 f n f n + i,1h 4 Bn,2 f n + 3 3 II +i,4 h5 f n4 f n + [ac0 (i + 1 s2 )(c1 + b0 c0 ) + II + ab0b1c0 ]h5 f n f n f n f n2 + i,5 h5 f n 2 f n f n2 + 3 II + a(i + 1 s2 )c0 + b0 c0 h5 f n f n f n / 6 + i,6 h5 f n2 Bn,1 f n + 3 3 +i,7 h5 f n Bn,1 f n f n + (i,2 + c1 )h5 f n Bn,2 f n + ( ) II 2 II h5 f n f n f n f n f n / 2 + c0b0 c0 h5 f n 2 f n / 2 + + c1 + b0c ( ) +c0 c1 + b0 c0 h5 f n f n f n / 2 + c0 h5 f nIV f n4 / 24 + II 2 3 +i,8 h5 Bn,1 f n 2 f n + (i + 1 s2 )ac0 h5 Bn,1 f n f n2 / 2 + + i,2 + (i + 1 s2 )ac0 h5 Bn,2 f n f n + i,1h5 Bn,3 f n + +i,2 h5 Bn,1 f n + c0 (c2 + ab0 d1 )h5 f n f n2 f n2 + c0 c1h5 f n Bn,1 f n2 / 2 + 2 +c0c1h5 f n f n Bn,1 f n / 2 + O (h6 ), (8.39) где i, j – некоторые постоянные, bi и i, j заданы в предыдущих леммах, а ci и di определяются следующим образом:

s2 1 s2 ( j + 1 s1 ) s,,jII, d 2 = ( j + 1 s1 )( j s1 )s2,,jII, III III d1 = j = s1 j = s s2 1 s2 1 s2 s,,jI, c0 = s,,jII, ciI = s2,,jI j,i, I III II III III c0 = 0 i 3, 2 j =1 j = s1 j = s2 s2,,jII j,i, 0 i 3, ci = ai ciI + ciII, 0 i 3.

ciII = III j = s 8.5. Теорема о максимальном порядке точности (m, 3)-методов Доказательство проведем методом математической индукции. Для III этого получим выражение для стадии k s. Разлагая (8.34) и (8.36), за пишем s2 1 s2 s2,,iI kiI + s2,,iII kiII III III = i =1 i = s ( ) II 2 II = c0 hf n + c1 + b0 c0 h 2 f n f n + b0 c0 h3 f n f n2 / 2 + ( ) +(c2 + ab0 d1 )h3 f n 2 f n + c1h3 Bn,1 f n + c3 + a 2b0 d 2 / 2 + a 2b0 d II h 4 f n3 f n + ab0b1c0 h 4 f n f n f n2 + ab0 d1h 4 f n f n f n2 / 2 + 3 II +b0 c0 h 4 f n f n / 6 + c2 h 4 f n Bn,1 f n + 3 ( ) II I + c5 + a 2c2 h 4 Bn,1 f n f n + c1h 4 Bn,2 f n + O(h5 ).

s2 1 s2 Разлагая функцию hf yn + s,,iI kiI + s,,iII kiII III III в ряд Тейлора в 2 i =1 i = s окрестности yn и выписывая члены до h5, получим s2 1 s2 hf yn + s,,iI kiI + s,,iII kiII III III = hf n + c0 h 2 f n f n + 2 i =1 i = s ( ) 12 1 2 II + c1 + b0c0 h3 f n 2 f n + c0 h3 f n f n2 + b0 c0 h 4 f n f n f n2 + II 2 +(c2 + ab0 d1 )h 4 f n3 f n + c1h 4 f n Bn,1 f n + ( ) ( ) 13 + c0 c1 + b0c0 h 4 f n f n f n2 + c0 h 4 f n f n + c3 + a 2b0 d 2 + a 2b0 d II 6 h5 f n4 f n + ab0b1c0 h5 f n f n f n f n2 + ab0 d1h5 f n 2 f n f n2 + II ( ) 1 3 II + b0 c0 h5 f n f n f n + c2 h5 f n 2 Bn,1 f n + c5 + a 2 c 3 II I Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ ( ) II h5 f n Bn,1 f n f n + c1h5 f n Bn,2 f n + c1 + b0 c0 h5 f n f n f n f n f n + 1 2 II + b0 c0 c0 h5 f n 2 f n + c0 (c2 + ab0 d1 )h5 f n f n2 f n + 3 ( ) 1 1 + c0 c1h5 f n Bn,1 f n2 + c0 c1h5 f n f n Bn,1 f n + c0 c1 + b0 c II 2 2 1 4 5 IV h 5 f n f n f n + 3 c0 h f n f n + O( h6 ).

Умножая данное выражение слева на Dn 1, запишем k s = hf n + ( a + c0 )h 2 f n f n + c0 h3 f n f n2 + III ( ) II + ac0 + ac1 + b0 c0 + a 2 h3 f n 2 f n + ah3 Bn,1 f n + ( ) II + c2 + ac1 + a 2c0 + ab0 d1 + ab0 c0 + a3 h 4 f n3 f n + ( ) ( ) ( ) ac0 + b0 c0 h 4 f n f n f n2 + ac1 + a 2 h 4 f n Bn,1 f n + c0 c1 + b0c 2 2 II II + ( ) h 4 f n f n f n2 + c0 h 4 f n f n + a 2 + ac0 h 4 Bn,1 f n f n + ah 4 Bn,2 f n + + c3 + ac2 + a 2c1 + a3c0 + a 2b0 d 2 + 2a 2b0 d1 + b0 c0 + a 4 h5 f n4 f n + II ( ) ( ) + ab0b1c0 + ac0 c1 + b0c0 h5 f n f n f n f n2 + ab0 d1 + ab0 c0 + a 2 c 2 2 II II II ( ) ( ) 1 3 II h5 f n2 f n f n2 + b0 c0 + ac0 h5 f n f n f n + c2 + ac1 + a3 h5 f n 2 Bn,1 f n + ( ) ( ) + c5 + a 2c2 + a 2 c0 + a3 h5 f n Bn,1 f n f n + c1 + a 2 h5 f n Bn,2 f n + II I 8.5. Теорема о максимальном порядке точности (m, 3)-методов ( ) II 1 1 2 II c1 + b0 c0 h5 f n f n f n f n f n + b0 c0 c0 h5 f n 2 f n + + 2 ( ) 12 + c0 c1 + b0 c0 h5 f n f n f n + c0 h5 f nIV f n4 + II 2 ( ) + ac1 + a 2 c0 + ab0 c0 + a3 h5 Bn,1 f n 2 f n + ac0 h5 Bn,1 f n f n2 + II ( ) + ac0 + a 2 h5 Bn,2 f n f n + ah5 Bn,3 f n + a 2 h5 Bn,1 f n + +c0 (c2 + ab0 d1 )h5 f n f n2 f n2 + c0c1h5 f n Bn,1 f n2 + + c0 c1h5 f n f n Bn,1 f n + O (h6 ).

III Получили, что k s представима в виде (8.39). Домножим (8.39) слева на Dn 1 :

kiIII = hf n + (i,1 + a + c0 )h 2 f n f n + c0 h3 f n f n2 + + ( ) + i,2 + ai,1 + c1 + ac0 + (i + 1 s2 ) ac0 + b0 c II h3 f n 2 f n + (i,1 + a)h3 Bn,1 f n + i,3 + ai,2 + a 2 i,1 + a3 + ac1 + ( ) + a 2 c0 + a (i + 1 s2 ) ac0 + b0 c0 h 4 f n3 f n + ac0 (i + 1 s2 ) + ac0 + 2 II ( ) +b0 c0 h 4 f n f n f n2 + i,2 + ai,1 + a 2 + c1 h 4 f n Bn,1 f n + 2 II ( ) +c0 c1 + b0 c0 h 4 f n f n f n2 + c0 h 4 f n f n + II Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ + i,2 + ai,1 + a 2 + (i + 1 s2 )ac0 + ac0 h 4 Bn,1 f n f n + + (i,1 + a )h 4 Bn,2 f n + i,4 + ai,3 + a 2 i,2 + ( i,1 + a) h 4 Bn,2 f n + + i,4 + ai,3 + a 2 i,2 + a3i,1 + a 4 + a 2 c1 + ( ) II + a3c0 + a 2 (i + 1 s2 ) ac0 + b0c0 h5 f n 4 f n + ( ) ( ) + ac0 (i + 1 s2 ) c1 + b0c0 + ab0b1c0 + ac0 c1 + b0 c II II II 2 II h5 f n f n f n f n2 + i,5 + a 2 c0 (i + 1 s2 ) + a 2 c0 + b0 c0 h5 f n 2 f n f n2 + + a(i + 1 s2 )c0 + b0 c0 + ac0 h5 f n f n f n + 3 3 II 3 ( ) + i,6 + ai,2 a 2 i,1 + a3 + ac1 h5 f n 2 Bn,1 f n + + i,7 + ai,2 + a 2i,1 + a3 + (i + 2 s2 )a 2 c0 h5 f n Bn,1 f n f n + ( ) ( ) + i,2 + ai,1 + a 2 + с1 h5 f n Bn,2 f n + c1 + ac0 h5 f n f n f n + 3 ( ) + i,6 + a i,2 a 2 i,1 + a 3 + ac1 h5 f n2 Bn,1 f n + ac0 h5 f n f n f n + 3 + i,6 + ai,2 a 2 i,1 + a3 + ac1 h5 f n2 Bn,1 f n + ( ) 12 + c0 c1 + b0 c0 h5 f n f n f n + c0 h5 f nIV f n4 + II 2 ( ) + i,8 + ai,2 + a 2 i,1 + a3 + ac1 + a 2c0 h5 Bn,1 f n 2 f n + 8.5. Теорема о максимальном порядке точности (m, 3)-методов + (i + 2 s2 )ac0 h5 Bn,1 f n f n2 + i,2 + ai,1 + a 2 + 2 +(i + 2 s2 )ac0 h5 Bn,2 f n f n + (i,1 + a )h5 Bn,3 f n + ( ) + i,2 + ai,1 + a 2 h5 Bn,1 f n + c0 (c2 + ab0 d1 )h5 f n f n2 f n2 + 2 1 + c0c1h5 f n Bn,1 f n2 + c0c1h5 f n f n Bn,1 f n + O( h6 ).

2 Итак, стадии метода (8.29) представимы в виде (8.39), где коэффици енты i, j, i s2, определяются следующим образом:

i +11 = i,1 + a, s2,1 = a, i +1,2 = i,2 + ai,1 + a 2, s2,2 = a 2,, ( ) i +1,3 = i,3 + ai,2 + a 2 i,1 + a3 + ac1 + a 2 c0 + a (i + 1 s2 ) ac0 + b0c0, II s2,3 = c2 + ac1 + a 2 c0 + a3 + ab0 d1 + ab0 c0, II i +1,4 = i,4 + ai,3 + a 2 i,2 + a3i,1 + a 4 + ( ) + a 2 c1a3c0 + a 2 (i + 1 s2 ) ac0 + b0 c0, II s2,4 = c3 + ac2 + a 2 c1 + a3c0 + a 4 + a 2b0 d 2 + 2a 2b0 d1 + a 2b0 c0, II 1 2 1 2 II i +1,5 = i,5 + (i + 1 s2 )a 2 c0 + b0 c0, 2 ( ) ab0 d1 + a 2c0 + ab0 c0, 2 2 2 II s2,5 = i +1,6 = i,6 + ai,2 + a 2i,1 + a3 + ac1, s2,6 = c2 + ac1 + a3, i +1,7 = i,7 + ai,2 + a 2 i,1 + a3 + (i + 2 s2 )a 2 c0, Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ s2,7 = a 2c2 + c5 + a 2 c0 + a3, I II ( ) i +1,8 = i,8 + ai,2 + a 2 i,1 + a3 + ac1 + a 2 c0 + a (i + 1 s2 ) ac0 + b0c0, II s2,8 = ac1 + a 2 c0 + a3 + ab0 c0.

II Лемма доказана.

Теорема 8.2. Максимальный порядок точности произвольного (m, 3) -метода (8.26) при условии использования замороженной или численной матрицы Якоби равен 4.

Для доказательства разложим точное решение задачи (8.25) в ряд Тейлора в окрестности точки t0, который имеет вид ( ) 1 y (tn+1 ) = y (tn ) + hf + h 2 f + h3 f 2 f + f 2 + f f 2 ( ) h f 3 f + f 2 + 3 f 2 + f 3 + + ff ff f h ( f f + f 2 f 2 + 4 f 2 f 3 + f 3 + + f ff +3 f f 2 + f 4 + 6 f 3 + 3 f + 4 f f 2 f 2 ) + O( h6 ). (8.40) ff f ff f ff f Применяя леммы 1, 3 и 4, запишем требования пятого порядка для диф ференциалов h3 f n f n2, h 4 f n f n f n2, h 4 f n f n, h5 f n 2 f n, h5 f nIV f n4, т. е.

3 s2 1 2 m i pi c0 + b0 j,s2 1 =, h3 f n f n2 : b 2 pi + (8.41) i = s1 i = s2 j = s2 s2 h f n f n f n2 : ab 4 pi (i + 1 s1 ) + i = s 2 m i pi ac0 (i + 1 s2 ) + b0 c0 + ab0 j,s2 1 ( j + 1 s1 ) =, (8.42) 2 II + i = s2 j = s 8.5. Теорема о максимальном порядке точности (m, 3)-методов s2 1 3 3 i m pi c0 + b0 j, s2 1 =, h 4 f n f n : b 3 pi + (8.43) i = s2 i = s1 j = s m h5 f n 2 f n : c0b0 c 3 2 II pi =, (8.44) i = s s2 1 4 m i pi c0 + b0 j,s2 1 =.

h5 f nIV f n4 : b 4 pi + (8.45) i = s1 i = s2 j = s2 Поскольку в разложении точного решения (8.40) отсутствуют эле ментарные дифференциалы с участием матриц Bi,n, имеем несколько дополнительных условий на пятый порядок точности метода. Одно из них, соответствующее элементарному дифференциалу h5 Bn,1 f n f n2, имеет вид s2 1 m ( ) 2 pi (i + 1 s1 ) + pi ac0 (i + 1 s2 ) + ab i = s1 i = s (8.46) i + ab0 j,s2 1 ( j + 1 s1 ) = 0.

j = s Вычитая (8.46) из (8.42), запишем m 2 II pi = b0 c0. (8.47) i = s Подставляя (8.47) в (8.44), имеем c0 = 4 / 5. Умножим (8.41) на b и вычтем из (8.43). С учетом c0 = 4 / 5 получим m 16 4 b0 pi = b0. (8.48) 25 5 i=s Умножим (8.43) на b0 и вычтем из (8.45). В результате запишем m 4 16 4 1b b0 pi = 0. (8.49) 5 25 5 i=s Г л а в а 8. ИССЛЕДОВАНИЕ (m, 3)-МЕТОДОВ Поделив (8.49) на (8.48), приходим к соотношению b0 = 0, что про тиворечит (8.47). Деление правомерно при выполнении соотношений m 4 pi 0;

b0 ;

b0.

5 i = s Первое условие выполняется в силу (8.47), два других проверяются непосредственной подстановкой в (8.48) и (8.49). Из равенства b0 = следует несовместность уравнений (8.41)–(8.46). Поэтому порядок точности (m, 3) -метода (8.26) при условии (8.27) не может быть выше четвертого. Пример L -устойчивого (m, 3) -метода четвертого порядка точности при требовании (8.27) известен. Теорема доказана.

Глава ГИБРИДНЫЕ СИСТЕМЫ Т ерминология гибридных систем в современной литературе рассматривается в основном на содержательном уровне и достаточно противоречива в разных источниках. Ниже вводятся стро гие математические определения ГС и системообразующих категорий.

На большом количестве простых и доступных для понимания приме ров типичных гибридных систем введенные определения иллюстри руются детально.

9.1. ДИСКРЕТНО-НЕПРЕРЫВНАЯ МОДЕЛЬ Гибридной называется система, которая отличается как непрерыв ным, так и дискретным поведением. Непрерывное поведение системы на временном интервале [t0, tk ] характеризуется вектором состояния x(t ) R n и задается отображением C : R m R m, где R – множество вещественных чисел. В дальнейшем для краткости вектор-функцию с начальным условием x0 = x(t0 ) будем обозначать x(t, x0 ). Для нее в общем случае выполняются условия:

• функция x(t, x0 ) непрерывна по совокупности переменных;

• x(t, x0 ) = x0 ;

• x ( t2 ;

x(t1, x0 ) ) = x(t0 + t1 + t2, x0 ).



Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |   ...   | 9 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.