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

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

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


Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 14 |

«Светлой памяти моих ро- дителей Марии Ивановны и Сергея Дмитриевича по- свящается В.С. ...»

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

( (3.2.5) r = Если последовательность значений векторов x(0), x(1), x(2), …, x(k), … является сходящейся, то есть у нее существует предел = lim x (k ), то этот предел является решением системы k (3.2.4). Действительно, имеем, что ( ) = lim x (k +1) = lim + x (k ) = +. (3.2.6) k k В связи с тем, что вектор является решением системы (3.2.4), которая получена из системы (3.2.1), то будет являться решением и этой системы.

Требования (3.2.6) сходимости итерационного процесса на хождения значений x 1, x 2, …, xn накладывает жесткие условия на коэффициенты системы (3.2.1). Однако, если detA 0, то с помощью элементарных преобразований системы (3.2.1) ее можно заменить эквивалентной системой (3.2.3), такой, что ус ловия сходимости будут выполнены.

В работе [104] приводится следующее условие сходимости итерационного процесса (3.2.5):

n jr 1.

max 1 j n r = Итерационный процесс заканчивается, если для двух при ближений в расчетной схеме (3.2.5) выполнено условие:

n x (jk +1) x (jk ) (k ) =, (3.2.7) j = где (k) – точность вычисления корней на k-м шаге;

– задан ная точность решения системы (3.2.1).

Пусть это условие выполнилось на шаге N = k + 1. Тогда в качестве корней системы (3.2.1) принимаются значения.

3.2.2. Метод Зейделя Данный метод представляет собой некоторую модифика цию метода простых итераций. Основная его идея состоит в том, что при вычислении (k +1)-го приближения неизвестной x i учитываются уже вычисленные ранее (k +1)-е приближения не известных x 0, x 1, …, x i–1. Иначе говоря, найденное (k +1)-е при ближение сразу же используется для получения (k +1)-го при ближения последующих координат (Рис. 3.2).

(k +1) (k +1) x x … i (k +1) xi (k +1) (k +1) … xi +1 xn Рис. 3. Предположим, что k-е приближения корней системы (3.2.3) известны. Тогда (k +1)-е приближения корней будут на ходиться по следующим итерационным формулам метода Зей деля:

n x1k +1) = 1 + ( 1 j x(jk ) ;

k = 0, 1, 2, j = n x2k +1) = 2 + 21x1k +1) + ( ( 2 j x(jk ) ;

j =..................................................... (3.2.8) i 1 n (k +1) = + x (k +1) + x (k ) ;

x i ij j ij j i j =1 j =i........................................................

n (k +1) = + x (k +1) + x (k ) +.

n nj j xn nn n n j = Расчетная схема (3.2.8) этого метода в общей форме имеет вид:

j + jr xr k ), j = (1, n ), k = 0, 1, 2,...(3.2.9) n x (jk +1) = j + jr xr k +1) ( ( r =1 r =i Достаточные условия сходимости метода Зейделя опреде ляются выполнением для всех и параметра q следующего неравенства [104]:

aij q aii, (3.2.10) j i где a ij – элементы матрицы А.

Условие завершения итерационного процесса (3.2.9) опи сывается неравенством (3.2.7). Таким же образом определяются искомые значения корней системы уравнений (3.2.1).

Пример 3. Построить расчетную схему метода простых итераций и метода Зейделя для численного решения следующей системы уравнений:

8 x 5 y + z = 1;

x + 6 y 2z = 7;

(3.2.11) x y + 4 z = 9.

Заметим, что эта система имеет точное решение x = 1;

y = 2;

z = 3.

Из матрицы a11 a12 a13 8 5 A = a21 a22 a23 = 1 6 a31 a32 a33 1 1 видно, что модули диагональных коэффициентов a ii в каждом уравнении отличны от нуля и больше суммы модулей всех ее остальных коэффициентов. Разделим каждое уравнение систе мы (3.2.11) на соответствующий диагональный коэффициент, сформируем столбец (x, y, z) в левой части и перенесем осталь ные слагаемые в правую часть. Проведя эти действия, получим расчетную схему вида (3.2.5) метода простых итераций для ре шения системы (3.2.11):

x (k +1) = + y (k ) z (k );

15 88 y (k +1) = x (k ) + z (k );

71 (3.2.12) 66 z (k +1) = + x (k ) + y (k ), k = 0,1,2,.

91 44 Начальное приближение в этой итерационной схеме обыч но выбирают равным столбцу свободных членов преобразован ( ) 1 7 ной системы (3.2.12), то есть x (0), y (0), z (0) =,,.

8 6 В связи с тем, что 5 1 1 1 1 1 max jr = max +, +, + = 1, 8 8 6 3 4 4 1 j 3 r = приведенное в Разд. 3.2.1 условие сходимости итерационного процесса (3.2.12) выполняется.

Расчетная схема метода Зейделя (3.2.9) для решения систе мы (3.2.11) запишется как:

x (k +1) = + y (k ) z (k );

15 88 y (k +1) = x (k +1) + z (k );

71 (3.2.13) 66 z (k +1) = + x (k +1) + y (k +1), k = 0,1,2,.

91 44 Здесь, как и выше, начальное приближение будет иметь вид:

(x(0), y (0), z (0) ) = 1, 7, 9 (0,125;

1,667;

2,250 ).

8 6 Для проверки сходимости метода Зейделя при решении системы (3.2.11) запишем неравенства, входящие в условия (3.2.10). С учетом состава элементов матрицы А эти неравенст ва примут вид:

5 + 1 q 8;

1 + 2 q 6;

1 + 1 q 4.

Эти неравенства выполняются при q = 0,9 1. Последнее означает, что итерационный процесс (3.2.13) будет сходиться к точному решению системы (3.2.11).

Результаты вычисления с точностью = 0,001 по расчетной схеме (3.2.12) приведены в табл. 3.1.

Из таблицы видно, что условия (3.2.7) выполняются на 7-м шаге. Приближенные значения корней равны x = 0,999;

y = 2,000;

z = 2,999.

Таблица 3.1.

(k) (k) (k) (k) k x y z 0 1,125 1,667 2,250 ––– 1 0,573 1,896 2,573 1, 2 0,988 1,929 2,867 0, 3 0,972 1,958 2,979 0, 4 0,976 1,998 2,982 0, 5 1,001 1,998 2,993 0, 6 1,000 1,998 3,000 0, 7 0,999 2,000 2,999 0, 3.3. Методы решения нелинейных алгебраических и трансцендентных уравнений Пусть дано нелинейное уравнение, которое в общем виде записывается как:

f ( x) = 0, (3.3.1) где f(x) – нелинейная алгебраическая или трансцендентная функция, либо функция, включающая в себя различные комби нации таких функций.

Число называется решением уравнения (3.3.1), если при его подстановке в левую часть (3.3.1) получается тождество.

Предполагается, что функция f(x) определена и непрерывна на некотором конечном или бесконечном интервале (a, b).

Будем считать, что уравнение (3.3.1) имеет только изолиро ванные корни, то есть для каждого k-го корня уравнения (3.3.1) существует окрестность:

k xk k, k = 0,1, 2,, которая не содержит других корней этого уравнения.

Процесс нахождения изолированных корней уравнения (3.3.1) можно разделить на два этапа:

1) отделение корней, то есть нахождение интервалов ( k ;

k ), в которых находится один и только один корень;

2) уточнение корней, которое состоит в определении при ближенных значений корней с определенной степенью точности.

Метод простых итераций. Пусть дано уравнение (3.3.1) и известно, что f(x) имеет на отрезке [a, b] один корень.

Уравнение (3.3.1) заменим эквивалентным уравнением вида:

x = (x). (3.3.2) Выберем некоторое начальное приближение ис помощью итерационной формулы:

xi +1 = ( xi ), i = 0, 1, 2,... (3.3.3) будем вычислять значения x 1, x 2, ….

Метод итераций, заданный расчетной схемой (3.3.3), при условии, что ( x) 1, (3.3.4) сходится к единственному решению уравнения (3.3.2), следова тельно, и (3.3.1) независимо от выбора начального приближе ния x 0 [a, b].

При невозможности непосредственного получения из урав нения (3.3.1) выражения вида (3.3.2) применяется следующий прием.

Умножим обе части выражения (3.3.1) на параметр C и прибавим к ним переменную x. В результате получим равенст во вида:

x + C f ( x) = x. (3.3.5) Обозначив левую часть полученного соотношения через (x), получим формулу (3.3.2).

Сформируем, исходя из формулы (3.3.5), ограничение, на кладываемое на значение постоянной C, такое, что получаемый итерационный процесс будет сходящимся.

Как было указано выше, сходимость итерационного про цесса (3.3.3) будет иметь место при выполнении условия (3.3.4).

Кроме того, заметим, что (x) = x + C f(x), поэтому имеем:

1 + C f ( x) 1.

Проводя несложные преобразования, получаем следующие условия для выбора значения параметра С:

;

0 при f ( x ) 0, x [a, b];

C f ( x ) при f ( x ) 0, x [a, b].

C 0;

f ( x ) При выборе значения С в этих выражениях полагают х = х 0.

После выбора конкретного значения С 0 расчетная схема рассматриваемого метода будет иметь следующий вид:

x i+1 = x i + C 0 f (xi ), i = 0, 1, 2,… Условия окончания итерационного процесса, заданного формулой (3.3.3) на некотором шаге N = i + 1 подразумевают одновременное выполнение следующих неравенств:

xi +1 xi ;

f ( xi +1 ), где и – заданные значения параметров точности решения уравнения (3.3.1).

При этом в качестве корня этого уравнения принимается значение x N.

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

Метод Ньютона (метод касательных). Пусть необходимо решить уравнение (3.3.1). Обозначим через [a, b] корень этого уравнения. Будем считать, что на некотором N-м шаге вычислительного процесса выполняется условие:

xN.

где – заданная точность вычисления корня.

Пусть x k – k-е приближение корня, то есть:

= xk + k, (3.3.6) где k = – x k – ошибка вычислений k-го приближения.

Подставляя (3.3.6) в (3.3.1), получим:

0 = f ( ) = f ( xk + k ) = f ( xk ) + f ( xk ) k, Разрешая это уравнение относительно k, имеем:

f ( xk ) k =. (3.3.7) ( xk ) f Подставляя (3.3.7) в (3.3.6), получаем расчетную схему ме тода Ньютона:

f ( xk ) xk +1 = xk, k = 0, 1, 2,... (3.3.8) f ( xk ) По формуле (3.3.8) при задании начального приближения х получаем последовательность приближений: x 1, x 2, …, xk,… Если эта последовательность сходится, то, как и выше, некото рое N-е приближение считается корнем уравнения (3.3.1).

Метод Ньютона не является самоисправляющимся, то есть он сходится не при любом выборе начального приближения.

Обычно в качестве значения x 0 берется такой конец отрезка [a, b], для которого выполняется условие:

f ( x0 ) f ( x0 ) 0.

К достоинствам данного метода можно отнести его высо кую скорость сходимости, близкую к квадратичной. Кроме то го, расчётная формула (3.3.8) метода достаточно легко реализу ется на ПЭВМ.

Недостатки метода Ньютона состоят в том, что итерацион ный процесс сходится не при любом выборе начального при ближения, а также в том, что в процессе вычислений нужно следить, чтобы на отрезке [a, b] первая и вторая производные функции f(x) в точках x k сохраняли свой знак.

Кроме того, можно заметить, что метод Ньютона является частным случаем метода простых итераций. Действительно, обозначив правую часть формулы (3.3.8) через (x k ), получим расчётную формулу метода итераций (3.3.3).

Модифицированный метод Ньютона. Если на отрезке [a, b] функция f (x) в уравнении (3.3.1) такова, что f'(x) f'(x0), то в формуле (3.3.8) целесообразно принять f'(xk) f'(x0). В этом случае получаем выражение вида:

f ( xk ) xk +1 = xk, k = 0, 1, 2,... (3.3.9) ( x0 ) f Формула (3.3.9) является расчетной схемой модифициро ванного метода Ньютона.

Наглядно различие методов Ньютона и модифицированно го метода Ньютона представлено на Рис. 3.3.

y y y = f(x) y = f(x) f(x0) f(x0) f(x1) f(x1) f(x2) f(x2) f(x3) x x f(x3) x3 x2 0a 0a x3 x2 x x1 b = x0 b = x а б Метод Ньютона Модифицированный метод (метод касательных) Ньютона (метод секущих) Рис. 3. Примеры применения рассмотренных выше методов, а также другие методы решения уравнений вида (3.3.1) приведе ны в работе [26].

3.4. Методы решения систем алгебраических и трансцендентных нелинейных уравнений Классическими системами нелинейных уравнений называ ются выражения вида [26]:

f1 ( x1, x2,, xn ) = 0 ;

f 2 ( x1, x2,, xn ) = 0 ;

(3.4.1)...............................

f n ( x1, x2,, xn ) = 0, где хотя бы одна из функций fi (x 1, …, x n ), является нелинейной функцией переменных x 1, …, x n,.

Такие системы уравнений используются при формировании управлений БЛА в установившихся режимах их полетов (см.

Разд. 5.4).

Решение систем нелинейных уравнений вида (3.4.1) являет ся в общем случае задачей более сложной, чем решение систем линейных алгебраических уравнений и одного нелинейного уравнения вида (3.3.1) [26]. Следует отметить, что в настоящее время отсутствуют методы, которые гарантировали бы опреде ление всех решений любой системы нелинейных уравнений.

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

Для отделения корней системы (3.4.1) общих методов не существует. В прикладных задачах, исследователь обычно име ет представление, где примерно находятся корни системы.

Методы, которым посвящен данный раздел, исходят из то го, что задача отделения корней решена и имеется достаточно малая область n-мерного пространства, содержащая корень, подлежащий определению. Пусть функции fi (x 1, …, xn ), определены в некоторых областях п-мерного про n [8]. Тогда область = i, которая странства i, i = представляет собой пересечение отмеченных выше областей, является областью, где может находиться решение системы (3.4.1).

Наиболее распространенными на практике методами нахо ждения единственного решения систем нелинейных алгебраи ческих и трансцендентных уравнений являются метод простых итераций, метод Ньютона и модифицированный метод Ньюто на [26].

Метод простых итераций. Для применения этого метода от исходной системы (3.4.1) путем соответствующих преобра зований необходимо перейти к эквивалентной системе вида:

x1 = 1 ( x1,, xn );

x2 = 2 ( x1,, xn );

(3.4.2).............................

xn = n ( x1,, xn ).

Итерационный процесс поиска решения, определяемый формулами:

( ) x1 k +1) = 1 x1 k ),, xnk ) ;

( ( ( x2k +1) = 2 (x1 k ),, xnk ) );

( ( ( (3.4.3).....................................

( ) xnk +1) = n x1 k ),, xnk ), k = 0, 1, ( ( ( начинают, задав начальное приближение в некоторой области.

Достаточным условием сходимости итерационного процес са (3.4.3) является выполнение условий [26]:

i x 1, i = (1, n ), n j =1 j i j = (1, n ).

n x 1, i =1 j Распишем первое неравенство:

1 1 при i = 1, ++ + xn x x..............................................................

n n n 1 при i = n.

++ + xn x x Второе условие в развернутой форме имеет вид:

n 1 1 при j = 1, ++ + x x x..............................................................

n 1 1 при j = n.

++ + xn xn xn Рассмотрим один из способов приведения системы (3.4.1) к виду (3.4.2), допускающему сходящийся итерационный процесс.

Пусть задана система уравнений второго порядка:

f1 ( x, y ) = 0 ;

(3.4.4) f 2 ( x, y ) = 0.

Требуется привести ее к виду:

x = 1 ( x, y ), (3.4.5) y = 2 ( x, y ).

Умножим первое уравнение системы (3.4.4) на неизвестную постоянную, второе – на, затем сложим их и прибавим к обеим частям уравнения величину x. Получим первое уравне ние преобразованной системы (3.4.5) в виде:

x + f1 ( x, y ) + f 2 ( x, y ) = x, 1 ( x, y ) = x + f1 ( x, y ) + f 2 ( x, y ).

где Далее, умножим первое уравнение системы (3.4.4) на неиз вестную постоянную, второе – на, затем сложим их и приба вим к обеим частям уравнения величину y. Тогда второе урав нение преобразованной системы (3.4.5) будет иметь вид:

y + f1 ( x, y ) + f 2 ( x, y ) = y, 2 ( x, y ) = y + f1 ( x, y ) + f 2 ( x, y ).

где Неизвестные постоянные,,, определим из приведен ных выше достаточных условий сходимости:

1 2 1 1 и + + 1.

x x y y Запишем эти условия более подробно:

f f f f 1 + 1 + 2 + 1 + 2 1;

x x x x f f f f + 2 + 1 + 1 + 2 1.

y y y y Полагая равными нулю выражения под знаком модуля, по лучим систему, состоящую из четырех линейных алгебраиче ских уравнений с четырьмя неизвестными,,, :

f f 1+ 1 + 2 = 0;

x x f f 1 + 2 = 0;

x x (3.4.6) f1 f + = 0;

y y f f 1+ 1 + 2 = 0.

y y При таком выборе параметров,,, условия сходимости будут выполнены, если частные производные функций f 1 (x, y) и f 2 (x, y) будут изменяться не очень быстро в окрестности точки. Тогда, чтобы решить систему (3.4.4), нуж но задать начальное приближение и вычис f i f i лить значения производных и, x y (0) (0) (0) (0) (x,y ) (x,y ) в этой точке. В противном случае, вычисление,,, осуществляется на каждом k-м шаге итераций с применением численных методов, изложенных в Разд. 3.2, с использованием производных:

f i f i, i = (1,2 ).

, x (x (k ), y(k ) ) x (x (k ), y(k ) ) Метод простых итераций является самоисправляющимся, универсальным и простым для реализации на ПЭВМ. Если сис тема имеет большой порядок, то применение данного метода, имеющего медленную скорость сходимости, не рекомендуется.

В этом случае используют метод Ньютона, который имеет бо лее высокую скорость сходимости.

Пример 3. Построить расчетную схему метода простых итераций для решения следующей системы нелинейных уравнений 2-го порядка:

f1 ( x, y ) = x + y 5 = 0 ;

(3.4.7) f 2 ( x, y ) = xy 6 = при начальном приближении:

( x (0), y (0) ) = (2;

1). (3.4.8) Заметим, что точными решениями системы (3.4.7) являются значения x = 2, y = 3 и x = 3, y = 2.

Для построения рабочих формул метода простых итераций необходимо решить систему линейных алгебраических уравне ний (3.4.6). Для ее решения вычислим частные производные f1 f 2 f1 f при начальном приближении (3.4.8):

,,, x x y y f1 f1 f 2 f = = 1, = y (2;

1) = 1, = x (2;

1) = 2.

x y x y Тогда система (3.4.6) примет следующий вид:

1 + + = 0;

+ = 0;

+ 2 = 0 ;

1 + + 2 = 0.

Решением этой системы являются значения = –2, = = 1, = –1. Тогда рабочие формулы метода простых ите раций для решения системы нелинейных уравнений (3.4.7) примут вид:

x = x 2( x + y 5) + ( xy 6) ;

y = y ( x + y 5) ( xy 6).

Для преобразования в расчетную схему метода вида (3.4.3) эти формулы необходимо переписать в виде:

( )( ) x (k +1) = x (k ) 2 x (k ) + y (k ) 5 + x (k ) y (k ) 6 ;

( )( ) y (k +1) = y (k ) x (k ) + y (k ) 5 x (k ) y (k ) 6 ;

k = 0,1, 2, Метод Ньютона. Пусть требуется решить систему нелиней ных уравнений (3.4.1) с действительными функциями f1, f2, …, fn.

Совокупность их аргументов x 1, x 2, …, x n будем рассматри x x вать как n-мерный вектор x = 2, а совокупность функций x n f f f 1, f 2, …, f n – как вектор-функцию F = 2. Тогда система f n (3.4.1) в векторной форме запишется как F ( x ) = 0. (3.4.9) x (k ) (k ) Предположим, что найдено k-е приближение x (k ) = x (k ) x n одного из изолированных корней = 2 векторного уравне n ния (3.4.9). Тогда точный корень этого уравнения можно пред ставить в виде:

= x (k ) + (k ), (3.4.10) (k ) (k ) где ( k ) = 2 – погрешность определения корня.

(k ) n Подставляя выражение (3.4.10) в формулу (3.4.9), получим:

( ) F x (k ) + (k ) = 0. (3.4.11) Предположим, что функция непрерывно дифференци руема в некоторой области, содержащей и. Разложим ле вую часть уравнения (3.4.11) в ряд Тейлора по степеням, ограничиваясь линейными членами:

( )()() F x (k ) + (k ) = F x (k ) + F x (k ) (k ) = 0, (3.4.12) f1 f1 f x1 x2 xn где F ( x ) = = J ( x ) – матрица Якоби (яко f n f1n f n x xn 1 x2 биан) системы функций f 1, f 2, …, f n относительно переменных x 1, x 2, …, x n [8].

Система (3.4.12) является системой линейных уравнений относительно погрешностей, с матрицей J, поэтому формулу (3.4.12) можно записать в виде:

()() F x (k ) + J x (k ) (k ) = 0. (3.4.13) Если матрица J невырожденная, то существует об ратная матрица J–1 и тогда, умножив обе части соотноше ния (3.4.12) на J–1 слева, получим методом обратной мат рицы следующее решение этой системы линейных уравнений:

( ) ( ).

( k ) = J 1 x ( k ) F x ( k ) Следовательно:

( )( ) x ( k +1) = x ( k ) J 1 x ( k ) F x ( k ), k = 0, 1, 2,... (3.4.14) Соотношение (3.4.14) является векторной записью расчет ной схемы метода Ньютона для решения систем нелинейных уравнений вида (3.4.1).

Если начальное приближение выбрано достаточно близко к решению системы (3.4.9), то итерационная последо вательность векторов сходится к этому ре шению, и сходимость является квадратичной.

Отметим недостатки метода Ньютона:

1) необходимость нахождения обратной матрицы J– на каждом шаге итерационного процесса (3.4.14), 2) возможность выхода приближения за пределы об ласти и связанная с этим расходимость итерационного процесса.

Модифицированный метод Ньютона. Данный метод лик видирует первый недостаток метода Ньютона. Если матрица J–1 непрерывна в окрестности искомого решения, и начальное приближение достаточно близко к, то приближен но можно считать J–1. Тогда формула (3.4.14) J– принимает вид:

( ) ( ).

x ( k +1) = x ( k ) J 1 x (0) F x ( k ) Эта формула является расчетной схемой модифицирован ного метода Ньютона.

Достоинством данного метода является то, что обратная матрица вычисляется только один раз. Но второй недостаток модифицированный метод Ньютона не ликвидирует.

Пример 3. Построить расчетную схему метода Ньютона для численно го решения системы нелинейных уравнений (3.4.7) при началь ном приближении (3.4.8).

Для нахождения обратной матрицы в формуле (3.4.14) не обходимо выполнить следующие действия:

1) найти матрицу частных производных:

f1 f x y 1 J ( x, y ) = = ;

f 2 f 2 y x x y 2) вычислить определитель этой матрицы:

det J ( x, y ) = x y ;

3) сформировать обратную матрицу:

T x y 1 x J ( x, y ) = =.

det J ( x, y ) 1 1 x y y Проведя несложные преобразования с матрицами [17], по лучим для решения задачи конкретизацию формулы (3.4.14) в следующем виде:

(k +1) = x y + 5 x 6 ;

(k ) (k ) (k ) x x(k ) y (k ) (k +1) = x y 5 y + 6, k = 0,1, 2, (k ) (k ) (k ) y x (k ) y (k ) Примеры решения различных систем нелинейных уравне ний рассмотренными методами приведены в работе [26].

3.5. Методы решения нелинейных параметрических уравнений Нелинейные алгебраические и трансцендентные парамет рические уравнения используются в общем случае при получе нии решений вида (2.2.7) задач параметрического нелинейного программирования, рассмотренных в Разд. 2.2. Кроме этого, та кими уравнениями являются системы (2.3.23), (2.3.43), (2.3.47), (2.3.48), которые используются для нахождения векторов u(t) управления динамическими объектами методами теории обрат ных задач управления такими объектами.

Отметим, что в литературе по вычислительной математике отсутствуют специальные методы решения нелинейных пара метрических уравнений.

В общем виде система нелинейных параметрических урав нений задается следующими выражениями:

. (3.5.1) Здесь t – некоторый параметр, при наличии которого иско мые корни этой системы будут функциями вида:

Рассмотрение численных методов решения системы вида (3.5.1) начнем с разработки методов ее решения при п = 1.

Пусть необходимо решить следующее нелинейное парамет рическое уравнение:

F ( x, t ) = 0, t [t1, t k ], (3.5.2) где t – некоторый скалярный параметр, который может прини мать дискретное множество значений t 1, t 2, …, t k или изменять ся непрерывным образом в заданном интервале значений.

В первом случае для решения каждого уравнения:

F ( x, ti ) = 0, i = (1, k ), (3.5.3) можно воспользоваться численными методами решения нели нейных уравнений вида (3.3.1), изложенными в Разд. 3.3. Тру доемкость их использования будет пропорциональна числу k заданных значений параметра t. При непрерывном изменении значения t [t1, tk] существующий практический подход к ре шению уравнения (3.5.2) состоит в сведении ее к первому слу чаю, но со значительно более «густой» сеткой значений пара метра t. В этом случае трудоемкость решения задачи за счет увеличения используемых значений параметра t резко возраста ет.

Преимущество предлагаемого ниже метода состоит в со кращении числа итераций, необходимых для нахождения ре шения в каждом текущем узле сетки значений параметра t, за счет использования решения, полученного в предыдущем узле сетки [109].

Перепишем уравнение (3.5.3) в следующей форме:

Fi ( x) = 0, i = (1, k ), (3.5.4) Fi ( x) = F ( x, ti ).

где Предположим, что каждое уравнение из совокуп ности (3.5.4) имеет единственное решение i [a, b] и F i (x) F j (x), (i j) для всех x [a, b],. Используем результаты теоремы о сходимости метода Ньютона для реше ния этой задачи: «Если Fi (a)Fi (b) 0, а не равны нулю и сохраняют определенные знаки на отрезке [a, b], то, ис ходя из начального приближения x 0, удовлетворяющего нера венству, по методу Ньютона, заданного фор () Fi xni ) ( мулой, можно вычислить единственный ко = xni+ () () xni ) ( Fi xn (i ) рень i уравнения F i (x) = 0 с любой степенью точности» [26]. В качестве начального приближения x 0 выбирается тот конец от резка [a, b], которому отвечают ординаты F i (x 0 ) того же знака, что и.

Будем считать известным значение искомого корня x(i) урав нения Fi(x) = 0 при фиксированном значении. Положим i = i +1 и рассмотрим ситуацию, представленную на Рис. 3.4.

F Fi+1(x) (1) Li+1 ( x ) ( 2) Fi(x) Li+1 ( x ) (i +1) (i +1) x (i ) a b x x2 x Рис. 3. Раскладывая правую часть уравнения:

Fi +1 ( x) = 0 (3.5.5) в ряд Тейлора в окрестности точки [8] с сохранением линейных членов, получаем:

() ( )( ).

Fi +1 ( x) = Fi +1 x (i ) + Fi+1 x (i ) x x (i ) Приравнивая к нулю правую часть этого выражения, получаем уравнение касательной к кривой в точке x = x(i):

() ( )( ) Fi +1 x (i ) + Fi+1 x (i ) x x (i ) = 0.

Предполагая, что корень x(i) является нулевым приближени ем к искомому корню x(i+1), получим первое приближение к не му по формуле:

( ), Fi +1 x0i +1) ( x1i +1) x0i +1) где x0i +1) = x (i ).

= Fi+1 (x0i +1) ) ( ( ( ( Если выполняются неравенства:

( ) Fi +1 x1i +1), (i +1) x1i +1) x1i +1) x0i +1), ( ( ( ( где – требуемая точность решения уравнения (3.5.5), то дан ное приближение принимается за искомое решение уравнения F i+1(x) = 0.

В противном случае в точке строится касательная к кривой Fi+1 (x), которая описывается уравнением вида:

( ) ( )( ) Fi +1 x1i +1) + Fi+1 x1i +1) x x1i +1) = ( ( ( Из этого уравнения получаем формулу для вычисления второго приближения к корню x(i+1):

( ).

Fi +1 x1i +1) ( x2i +1) x1i +1) = Fi+1 (x1i +1) ) ( ( ( Обобщая этот процесс, запишем рекуррентную формулу для вычисления j-го приближения к искомому корню:

( ), + Fi +1 x (ji1 ) x (ji +1) + = j = 1,2,3,.

Fi+1 (x (ji1 ) ) x (ji1 ) (3.5.6) + Начальное условие для нее имеет вид:

x0i +1) = x (i ).

( Общие условия завершения на j-той итерации процесса по иска корня x(i+1) записывается как:

( ) Fi +1 x (ji +1), j = 1,2,3, (3.5.7) x (ji +1) x (ji +1) + (i +1) x (ji1 ).

Для случая, когда кривые F i (x) могут пересекаться, следует применить рекуррентную процедуру, основанную на комбини рованном численном методе, объединяющем метод Ньютона с методом хорд [26].

Аналогично этим выражениям записывается расчётная схе ма для решения системы уравнений (3.5.1).

Для случая, когда по условиям решаемой задачи параметр t должен непрерывно изменяться в заданном интервале [t 1,t k ], для решения уравнения (3.5.2) предлагается использовать ме тод, основанный на разложении левой части этого уравнения в ряд Тейлора [8] с сохранением линейных членов разложения [1]. Проводя это разложение в окрестности точки (t i-1, xi-1 ), ко торая является решением уравнения (3.5.2) при t = t i-1, имеем:

F F ( x xi 1 ) + (t ti 1 ) = 0. (3.5.8) x i 1 t i Отсюда корень уравнения, решаемого при t = ti, вычисляет ся как:

F t i xi = xi 1 (ti +1 ti ). (3.5.9) F x i Отметим, что для применения формулы (3.5.9) шаг сетки по параметру t должен быть достаточно малым.

Общая расчётная схема предлагаемого метода решения сис темы нелинейных параметрических уравнений вида (3.5.1) включает в себя следующие этапы [1]:

1. Строится сетка по параметру t с достаточно большим числом узлов, значения которых определяются как:

t t ti = ti 1 + t, i = (1, k ), t = k 0, N где N – число узлов сетки на отрезке [t 0, t k ].

2. Функции F j раскладываются в ряд Тейлора в окрестности ( ) точки ti 1, x1i 1), x2i 1),, xni 1) с сохранением только линейных ( ( ( членов разложения:

( ) ( ) (i 1), x (i 1),, x (i 1) + F j x (i ) x (i 1) + n x r r F j ti 1, x1 n r =1 r i (3.5.10) F j j = (1, n ), t (ti ti 1 ) = 0, + i где ( ), F j F j ti 1, x1i 1), x2i 1),, xni 1) ( ( ( = x r i 1 xr ( ).

F j = F j ti 1, x1i 1), x2i 1),, xni 1) ( ( ( t i 1 t Заметим, что выражение (3.5.8) является частным случаем соотношения (3.5.10) при п = 1.

3. Система (3.5.10) записывается в виде:

n F F j xri ) = j t, j = (1, n ), x ( (3.5.11) t r =1 r i 1 i где xri ) = xri ) xri 1), t = ti ti 1, r = (1, n ).

( ( ( 4. Используя полученные одним из численных методов, рассмотренных в Разд. 3.2, решения, системы линейных уравнений (3.5.11) находятся искомые решения, исходной нелинейной системы урав нений (3.5.1) по следующим рабочим формулам:

xri ) = xri 1) + xri ), i = (1, k ), r = (1, n ).

( ( ( (3.5.12) Условие завершения на i-той итерации процесса поиска реше () () () записывается как F j xr ), i = 1, k, r = 1, n, j = (1, n ).

ния (i Преимущество данной процедуры состоит в сокращении числа итераций, необходимых для нахождения решения систе мы в каждом узле достаточно «густой» сетки значений пара метра t, за счет использования решения, полученного в преды дущем узле сетки.

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

Пусть требуется получить решение x = x(t), t [t 1, t k ] уравнения (3.5.2). Данный метод предполагает замену этого уравнения сингулярно возмущенным дифференциальным урав нением вида:

dx µ = F ( x, t ), (3.5.13) dt где 0 µ 1 – малый параметр.

Начальное условие при t = для этого уравнения запишется как:

х(, µ) = х0. (3.5.14) Применение метода малого параметра при решении урав нения (3.5.2) обосновывается следующим утверждением, полу ченным на основе результатов работы [110]: «Корень x* = x*(t), полученный путем решения задачи Коши (3.5.13), (3.5.14), бу дет устойчиво стремиться к искомому решению x0 = x0(t) при выполнении следующего условия:

, (3.5.15) где – частная производная функции F по переменной х».

Для выполнения этого условия выражение (3.5.2) может при необходимости быть переписано в форме:

–F(x, t) = 0.

Тогда уравнение (3.5.13) примет вид:

dx µ = F ( x, t ).

dt Начальная точка (, x 0 ) интегральной кривой x = x(t) урав нения (3.5.13) выбирается, исходя из условий:

(3.5.16) Промежуток (, t0), на котором реализуется быстрое измене ние решения задачи Коши (3.5.13), (3.5.14) от начального значе ния х 0 до значений, близких к x*(t), называется пограничным слоем. Множество значений (, х0 ), удовлетворяющих условиям (3.5.16), называется областью притяжения корня x*(t) [110].

На Рис. 3.5 представлено поведение решений задачи Коши (3.5.13), (3.5.14) при различных значениях и х0.

x0 = x0(t) x x* = x*(t) (, x0) (, x0) t1 tk t Рис. 3. Установлено, что при выборе любых начальных условий (3.5.14) из области притяжения можно получить путем интег рирования уравнения (3.5.13) приближенное решение x*(t) уравнения (3.5.2).

Это решение должно удовлетворять условию:

, (3.5.17) где – требуемая точность решения уравнения (3.5.2).

Заметим, что для применения соответствующих численных методов интегрирования (см. Разд. 3.1) уравнение (3.5.13) пред ставляется в следующем виде:

dx = F ( x, t ), t [, t k ]. (3.5.18) dt µ Таким образом, при использовании предлагаемого метода не обходимо выполнить следующую последовательность действий:

1. Выбрать достаточно малое значение параметра µ.

2. Установить область притяжения решения задачи и вы брать начальные значения и х 0, удовлетворяющие условиям (3.5.16).

3. Выбрать и реализовать численный метод интегрирова ния уравнения (3.5.18).

4. В процессе получения решения этого уравнения при t = t 1 проверять выполнение условия (3.5.17). При его выполне нии процесс интегрирования продолжается до достижения па раметром t значения tk.

В противном случае процесс останавливается и произво дится с использованием неравенств (3.5.16) выбор новых зна чений и х 0 с повторным интегрированием уравнения (3.5.18).

Последние действия выполняются до тех пор, пока не будет удовлетворено условие (3.5.17).

Аналогичный подход предлагается использовать при реше нии системы уравнений (3.5.1).

При его реализации формируется следующая система син гулярно возмущенных уравнений [110]:

dx j µ = F j (t, x1, x2,..., xn ), j = (1, n), t [, tk ] (3.5.19) dt с начальными условиями:

. (3.5.20) Для получения решения этой системы, приближенно представляющего точ ные решения системы уравнений (3.5.1), выбираются начальные значения, х0j,, удовлетворяющие неравенствам вида:

, (3.5.21) где – частная производная функции по переменной xj.

Здесь, как и выше, для выполнения этих неравенств в пра вой части уравнений (3.5.19) некоторые функции F j, могут быть умножены на (–1).

Величины, кроме выполнения условий, (3.5.21) должны удовлетворять неравенству:

F j (x1 (t1 ), x2 (t1 ),..., xn (t1 ), t1 ), n * * * (3.5.22) j = обеспечивающему решение системы уравнений (3.5.1) при с заданной точностью.

При решении этой системы методом малого параметра вы полняется указанная выше последовательность действий, в ко торых используются выражения (3.5.19)-(3.5.22). В п.п. 3 и этой последовательности используется представление системы дифференциальных уравнений (3.5.19) в следующей форме:

dx j = F j (t, x1, x2,..., xn ), j = (1, n), t [, t k ].

dt µ Основным достоинством применения метода малого пара метра для решения линейных и нелинейных параметрических уравнений является возможность снизить трудоемкость про цессов формирования искомых решений с привлечением чис ленных методов интегрирования дифференциальных уравнений (3.5.13) или (3.5.19).

3.6. Методы решения краевых задач для обыкновенных дифференциальных уравнений Задачи данного класса возникают при использовании на практике методов вариационного исчисления и оптимального управления, приведенных в Разд. 2.4 и 2.5.

Численные методы решения краевых задач применяются в тех случаях, когда невозможно получить аналитические (фор мульные) решения дифференциальных уравнений, получаемых этими методами.

В настоящее время для решения краевых задач широко при меняются следующие основные численные методы [25, 104]:

• метод «пристрелки»;

• конечно-разностный метод.

Проиллюстрируем расчетные схемы этих методов на при мерах краевых задач, описываемых выражениями (2.1.35) (2.1.38).

Метод «пристрелки». Суть этого метода состоит в много кратном решении уравнения (2.1.35) с изменяющимися началь ными условиями до тех пор, пока не будет удовлетворено с за данной точностью условие «попадания» интегральной кривой в заданную граничную точку (Рис. 3.6).

( 0) y y (1) y (k ) y у = у(х) y y 0 a b x Рис. 3. Физическая сущность метода «пристрелки», основанная на процессе корректировки артиллерийского огня, приведена в ра боте [104].

Пусть на заданном интервале требуется решить краевую задачу (2.1.35), (2.1.36). Для ее решения рассматриваемым ме тодом сформируем задачу Коши с дифференциальным уравне нием (2.1.35) и начальными условиями вида:

у(а) = у0 ;

у'(a) =, (3.6.1) где – вспомогательная переменная, определяющая тангенс угла наклона касательной к интегральной кривой в точке х = а.

Заметим, что изменение значения этой переменной позво ляет получить семейство интегральных кривых, представлен ное на Рис. 3.6.

Придадим некоторое начальное значение 0 этой перемен ной и решим одним из численных методов задачу Коши:

Обозначим через значение решения этой задачи в точке x = b.

Решая задачу Коши (2.1.35), (3.6.1) для другого значения = 1, получим значение Используя теорему о непрерывной зависимости решений дифференциального уравнения от значений используемых на чальных условий [20], можно утверждать, что решение задачи Коши (2.1.35), (3.6.1) в общем случае описывается функцией:

. (3.6.2) При x = b и заданном значении у 0 это соотношение будет функцией одной переменной.

Следуя этому, сформируем функцию вида:

, (3.6.3) которая описывает величину отклонения решения (3.6.2) в точке x = b от граничного значения у1, заданного выражением (2.1.36).

Отсюда возникает задача определения значения * такого, что (*) = 0.

Эта задача представляет собой задачу нахождения корня уравнения:

() = 0. (3.6.4) Особенностью этого уравнения является его алгоритмиче ский характер, который определяется тем, что функция (3.6.3) задается путем численного решения задачи Коши (2.1.35), (3.6.1) при конкретном значении вспомогательной переменной.

В связи с невозможностью вычисления производной функ ции () при решении уравнения (3.6.4) будем использовать расчетную схему (3.3.8) метода Ньютона, в которой производ ная приближенно представлена ее разностным аналогом [26], вычисляемым по двум приближениям k и k+1 искомого корня.

Вследствие этого расчетная схема решения уравнения (3.6.4) будет иметь вид:

k +1 k k + 2 = k +1 (k +1 ), k = 0,1, 2,... (3.6.5) (k +1 ) (k ) где в соответствии с выражением (3.6.3) имеем, что (k ) = y1 k ) y1, k = 0,1, 2,...

( (3.6.6) Отметим, что для применения этой схемы необходимо за даться начальными значениями 0 и 1.

Расчеты по формуле (3.6.5) выполняются до некоторого шага с номером N, при котором выполняется следующее усло вие останова итерационного процесса:

, (3.6.7) где – требуемая точность выполнения граничных условий (2.1.36).

Пример 3.4.

Пусть на интервале [0, 1] с точностью = 0,0001 требуется решить краевую задачу:

(3.6.8) y(0) = 1;

y(1) = 2, (3.6.9) которая является конкретизацией выражений (2.1.35), (2.1.36).

Представим дифференциальное уравнение (3.6.8) в виде системы (3.1.10):

y ' = z;

(3.6.10) z ' = e x + sin y.

Начальные условия для этой системы будут иметь следую щий вид:

y(0) = 1, z(0) =. (3.6.11) Задачу Коши (3.6.10), (3.6.11) будем решать с шагом интег рирования h = 0,1 с применением расчетной схемы (3.1.14), (3.1.15) метода Рунге-Кутта (см. Разд. 3.1.2).

Для использования расчетной схемы (3.6.5), (3.6.6) зада димся значениями 0 = 1,0 и 1 = 0,8.

Результаты применения метода «пристрелки» к решению краевой задачи (3.6.8), (3.6.9) приведены в табл. 3.2 и 3.3.

Таблица 3. ( k ) k 0 +1,00000 3,16889 1, 1 +0,80000 2,97448 0, 2 –0,20466 1,95376 0, 3 –0,15916 2,00179 0, 4 –0,16086 2,00003 0, Из этой таблицы видно, что условие (3.6.7) при заданном значении = 0,0001 выполняется на четвертом шаге (N = 4) итерационного процесса (3.6.5), (3.6.6).

Результаты решения рассматриваемой краевой задачи пред ставлены в табл. 3.3.

Таблица 3. xi 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1, y i 0,1 0,993 1,006 1,039 1,094 1,174 1,279 1,412 1,575 1,777 2, Данная таблица значений функции у = у(х) получена путем решения задачи Коши (3.6.10), (3.6.11) при = 4 = –0,16086.

Конечно-разностный метод. Пусть требуется решить на отрезке [a,b] краевую задачу для линейного неоднородного дифференциального уравнения второго порядка с переменными коэффициентами:

(3.6.12) и граничными условиями 1-го рода:

(3.6.13) Введем на отрезке [a,b] сетку значений x0 = a, x1, x2, …, xN = b, где x k = a + hk, k = 0, 1, …, N, h = (b – a) / N. Решение задачи (3.6.12), (3.6.13) будем определять в виде множества значений )}, где yk = y(xk). Представим производные, входя {yk, k = ( щие в уравнение (3.6.12), их разностной аппроксимацией вида [25]:

y yk 1 y 2 yk + yk = k +1 = k +1. (3.6.14) ;

2h h Подставляя выражения (3.6.14) в уравнение (3.6.12), полу чим систему уравнений для нахождения значений y k :

Здесь предполагается, что при использовании граничных условий (3.6.13) неизвестные y0 и y N считаются заданными.

Приводя подобные члены, получим следующую систему ли нейных алгебраических уравнений относительно неизвестных ( 2 + h2q( x1))y1 + 1 + p( x1)h y2 = p ( x1 )h = h 2 f ( x1 ) 1 ya ;

.............................................................................

( ) p ( xk ) h 1 yk 1 + 2 + h q ( xk ) yk + 2 (3.6.15) p ( xk ) h + 1 + yk +1 = h f ( xk );

k = (2, N 2);

.............................................................................

( ) p ( x N 1 )h 1 y N 1 + 2 + h q ( x N 1 ) y N 1 = p ( x N 1 )h = h 2 f ( x N 1 ) 1 + yb.

Для решения системы (3.6.15) при достаточно малых значе ниях шага сетки h и q(x k ) 0 применяется один из численных методов, приведенных в Разд. 3.2.

В случае использования граничных условий 2-го и 3-го ро да, описываемых выражениями (2.1.37), (2.1.38), аппроксима ция производных, входящих в соответствующие краевые зада чи в граничных точках, проводится с помощью следующих од носторонних разностей 1-го и 2-го порядков [25]:

y y0 y y N =N =1 ;

;

(3.6.16) h h 3 y0 + 4 y1 y 2 4 y N 1 + 3 y N y = N = ;

;

(3.6.17) 2h 2h В случае использования формул (3.6.16) линейная алгеб раическая система (3.6.15) аппроксимирует краевую задачу с первым порядком точности. При использовании формул (3.6.17) получается второй порядок точности решения краевой задачи.

Пример 3. Решить на интервале [0, 1] с шагом h = 0,2 краевую задачу:

y" – xy' – y = с граничными условиями:

y(0) = 1;

y'(1) + 2y(1) = 0.

Данное уравнение является частным случаем уравнения (3.6.12) при следующих значениях его коэффициентов и правой части:

p(x) = x;

q(x) = 1;

f(x) = 0.

Параметры решаемой задачи имеют следующие значения:

N = 5;

x0 = 0;

x 1 = 0,2;

x 2 = 0,4;

x 3 = 0,6;

x4 = 0,8;

x 5 = 1,0.

Во всех внутренних узлах интервала [0, 1] после замены производных их разностными аналогами (3.6.14) получим вы ражения вида:

Из заданных граничных условий следует, что на левой гра нице значение у0 = 1. На правой границе аппроксимируем про изводную y'(1) односторонней разностью 1-го порядка (3.6.16).

В этом случае второе граничное условие примет вид:

С помощью группировки слагаемых, приведения подобных членов, подстановки значений xk и с учетом того, что у0 = 1, полу чим следующую систему линейных алгебраических уравнений:

2,04 y1 + 1,02 y 2 = 0,98;

0,96 y1 2,04 y 2 + 1,04 y3 = 0;

0,94 y 2 2,04 y3 + 1,06 y 4 = 0;

0,92 y3 2,04 y4 + 1,08 y5 = 0;

y4 1,4 y5 = 0.

В результате решения этой системы методом простых ите раций (см. Разд.3.2.1) получим следующие значения искомых переменных:

y 1 = 0,7719;

y 2 = 0,5830;

y 3 = 0,4311;

y4 = 0,3126;

y 5 = 0,2233.

Решение краевой задачи приведено в табл. 3.4.

Таблица 3. k 0 1 2 3 4 xk 0 0,2 0,4 0,6 0,8 1, yk 1,0 0,77191 0,58303 0,43111 0,31265 0, Пусть на интервале [a, b] требуется решить краевую задачу для системы дифференциальных уравнений (2.1.9).

Будем считать, что граничные условия 3-го рода этой зада чи в точках x = a и x = b заданы в общем случае нелинейными функциями:

1 (a, y1a, y2 a, y2 a,..., yna ) = 0;

(3.6.18) 2 (b, y1b, y2b, y2b,..., ynb ) = 0.

Выберем величину шага h и введем сетку значений по правилу:

где N = (b – a)/h.

Решение краевой задачи (2.1.9), (3.6.18) конечно разностным методом, как и выше, заключается в определении множества значений yjk, где yjk = y j (x k ),.

Заменим производные в системе (2.1.9) их разностной ап проксимацией вида (3.6.14).

Проводя несложные преобразования, получаем с учетом вы ражений (3.6.18) следующую систему нелинейных уравнений:

1 ( x0, y10, y20,..., yn 0 ) = 0;

........................................

y j k +1 y j k 2hf j ( xk, y1k, y2 k,..., ynk ) = 0, k = (0, N 1);

(3.6.19).........................................

2 ( xN, y1N, y2 N,..., ynN ) = 0.

Решение этой системы уравнений может быть получено од ним из численных методов, приведенных в Разд. 3.4.

Система (3.6.19) значительным образом упрощается, если решения системы (2.1.9) должны удовлетворять граничным ус ловиям 1-го рода:

(3.6.20) Другие численные методы решения краевых задач для сис тем уравнений вида (2.1.9) и (2.1.11) приведены в работе [104].

3.7. Методы решения задач вариационного исчисления Приведенные в Разд. 2.4 дифференциальные уравнения Эй лера, с помощью которых осуществляется решение вариацион ных задач при сколь-нибудь сложном практическом примене нии, только в исключительных случаях интегрируются в анали тическом виде.

Как было показано в этом разделе, любая вариационная за дача является краевой задачей с различными видами граничных условий, рассматриваемых в Разд. 2.1. Поэтому для решения полученных уравнений Эйлера используются численные мето ды решения краевых задач, описанные в Разд. 3.6. В частности, для решения задачи (2.4.7), (2.4.8) применяется метод «при стрелки». Вариационные задачи с подвижными границами (см.

Разд. 2.4.2), которые сводятся к краевой задаче (2.4.4), (2.4.21), могут быть решены конечно-разностным методом, где в каче стве граничных условий (2.1.39) используются применяемые условия трансверсальности.

Практически одновременно с разработкой теоретических вопросов вариационного исчисления были начаты работы по разработке численных (приближенных) методов построения экстремалей, получившие название прямых методов решения вариационных задач [17, 20].

Характерной особенностью таких методов является отсутст вие необходимости формирования и решения уравнений Эйлера.

Для иллюстрации сути прямых методов рассмотрим два наиболее распространенных численных метода этого класса.

Конечно-разностный метод Эйлера. Пусть дана вариаци онная задача вида:

t J [x(t )] = F (t, x, x ) dt max ;

(3.7.1) t x(t0 ) = a ;

x(t1 ) = b.

Построим на интервале [t 0, t 1] сетку значений t 0 + t ;

t 0 + 2t ;

…;

t 0 + (n 1)t аргумента t функции x(t), где t = (t1 t 0 ) n.

Каждому узлу этой сетки поставим в соответствие значения ординаты функции x(t):

x1 = x(t 0 + t );

x2 = x(t 0 + 2t );

;

xn 1 = x(t 0 + (n 1)t ).

Отметим, что ординаты a и b являются заданными гранич ными условиями.

В этом случае искомая экстремаль x = x(t) заменяется лома ной Эйлера, представленной на Рис. 3.7.

x xn– x0 = a x1 x2 xn = b t0+t t0+2t t0 t0+(n–1)t t1 t Рис. 3. На такой ломаной функционал J[x(t)] превращается в функ цию J = f(x 1, x 2,…, x n–1 ) ординат x 1, x 2, …, x n–1 ее вершин.

Эти ординаты выбираются из условия достижения макси мума функции f (x 1, x 2,…, x n–1 ).

Для решения этой задачи будем использовать уравнения вида (2.2.4), которые в данном случае принимают вид:

J = 0, j = (1, n 1). (3.7.2) x j Для построения функции J заменим интеграл (3.7.1) его приближенным представлением, используемым в методе пря моугольников [8]:

t F (t, x, x ) dt F (t j, x j, x j )t.

n j = t Для производной будем использовать приближенное представление:

x j +1 x j xj.

t Тогда имеем, что x j +1 x j n J [x(t )] f ( x1, x2,, xn 1 ) = Ft j, x j, t.

t j =0 В этой сумме от x j зависят только (j – 1)-е и j-е слагаемые:

x j x j 1 x xj t и F t j, x j, j + F t j 1, x j 1, t.

t t Поэтому уравнения (3.7.1) могут быть представлены как:

x xj x x j t j, x j, j +1 t + Fx t j, x j, j +1 t + Fx t t t x x j 1 t = 0, j = (1, n 1).

t j 1, x j 1, j + Fx t t Проводя несложные преобразования, получим окончатель ный вид системы уравнений для определения значений x 1, x2, …, x n–1 искомых ординат экстремали x(t):

x j +1 x j x x j t + Fx t j 1, x j 1, j Fx t j, x j, t t (3.7.3) x j +1 x j = 0, j = (1, n 1).

Fx t j, x j, t В общем случае эта система нелинейных уравнений реша ется одним из численных методов, приведенных в Разд. 3.4.

Пример 3. Запишем систему уравнений (3.7.3) для решения конечно разностным методом Эйлера следующей вариационной задачи:

( ) J [x(t )] = 3tx + 2 x 2 dt max ;

x(0 ) = 1, x(5) = 4.

При F (t, x, x ) = 3 t x + 2 x 2 имеем:

Fx (t, x, x ) = 3t ;

Fx (t, x, x ) = 4 x.

Зададимся значением n = 100 и определим шаг сетки значе ний t:

t = = 0,05.

Узлы сетки вычислим по формуле:

j = (1,99 ).

t j = 0 + jt = 0,05 j, Производные функции x (t) в узлах сетки представим как:


x j +1 x j x(t j ) = 20( x j +1 x j ).

t Тогда система уравнений (3.7.3) примет следующий вид:

( x j x j 1 ) ( x j +1 x j ) = 0, j = (1,99 ).

3t j + 0, Проводя несложные преобразования, получим окончатель ный вид формируемой системы:

x j 1 + 2 x j x j +1 = 0,000938 j, j = (1,99 ).

Таким образом, ломаная Эйлера может быть построена пу тем численного решения этой системы линейных алгебраиче ских уравнений методами, изложенными в Разд. 3.2.

Аналогичным образом могут быть построены системы уравнений типа (3.7.3) для других видов вариационных задач.

Для получения с приемлемой для практики точностью ре шения вариационных задач методом Эйлера необходимо вы брать достаточно малую величину шага t. Это приводит к большой размерности системы уравнений (3.7.3) и, как следст вие, значительной трудоемкости ее численного решения.

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

Метод Ритца является наиболее часто применяемым методом приближенного решения различных видов вариационных задач.

Пусть задан функционал (2.4.4) и требуется найти такую функцию у(х), принимающую согласно выражениям (2.4.8) в точках х 0 и х1 заданные значения у 0 = у(х 0 ) и у1 = у(х 1 ), на ко торой функционал J[у(х)] будет достигать экстремума. Значе ния исследуемого на экстремум функционала (2.4.4) рассмат риваются не на всех допустимых в данной задаче функциях у(х), а на функциях, описываемых линейной комбинацией функций некоторой выбранной системы :

n y n ( x) = i i ( x), (3.7.4) i = где i – постоянные коэффициенты.

Отметим, что от удачного выбора этой системы функций зависит эффективность применения метода к решению кон кретных вариационных задач. Необходимым условием выбора подходящей системы функций является требование, что бы функции удовлетворяли условиям и для всех значений параметров i. Для выполнения этих условий в работе [20] предлагается использовать следую щую зависимость:

где функция 0 (х) имеет вид:

а функции i (х) выбираются так, чтобы i (х 0 ) = i (х 1 ) = 0.

Подставляя выражение (3.7.4) в функционал (2.4.4), имеем:

x J n = F ( x, yn ( x), y 'n ( x))dx = x n x1 n = F x, i i ( x), i 'i ( x) dx (3.7.5) x0 i =1 i = x = ( x, 1, 2,..., n )dx.

x Вычисляя по соответствующим правилам [8] последний ин теграл, получим функцию п переменных вида:

Jn = Jn. (3.7.6) Неизвестные значения коэффициентов будем определять из условия достижения экстремума функции (3.7.6).

Применяя к этой функции необходимые условия экстрему ма вида (2.2.4), запишем следующую систему уравнений:

J n = 0, i = (1, n). (3.7.7) i Решая эту систему определенным численным методом из Разд. 3.2 или Разд. 3.4, получаем значения искомых коэффици ентов и приближенное решение вариационной задачи вида:

n = i 0)i ( x), x [ x0, x1 ].

( yn0) ( x) ( (3.7.8) i = Пример 3. Пусть требуется решить следующую вариационную задачу:

J = ( y ' 2 2 xy )dx min;

(3.7.9) у(0) = 0;

у(1) = 0. (3.7.10) В качестве функции (3.7.4) выберем зависимость вида:

n y n ( x) = i (1 x) x i. (3.7.11) i = Нетрудно убедиться, что для такой функции граничные ус ловия (3.7.10) выполняются.

Будем решать задачу (3.7.9), (3.7.10) при п = 2. В этом слу чае выражение (3.7.11) конкретизируется как:

Для выполнения требуемых вычислений перепишем эту функцию в следующей форме:

. (3.7.12) Производная от функции (3.7.12) будет иметь вид:

Возводя это выражение в квадрат, имеем:

(3.7.13) Подставляя соотношения (3.7.12) и (3.7.13) в функционал (3.7.9), получим выражение вида (3.7.5):

1 J 2 = ( y2 2 xy2 )dx = [1 (1 4 x + 4 x 2 ) + 2 0 Вычисляя этот интеграл, сформируем функцию вида (3.7.6):

J2 = J С использованием этого выражения система уравнений (3.7.7) конкретизируется как:

J J Решая полученную систему линейных алгебраических уравнений, получим следующие значения:

С учетом этих значений и формул (3.7.12) решение задачи (3.7.9), (3.7.10) в виде выражения (3.7.8) конкретизируется как:

Вариационная задача (3.7.9), (3.7.10) имеет точное решение:

которое получается при интегрировании уравнения Эйлера:

у" + x = с граничными условиями (3.7.10). Абсолютное отклонение точ ного и приближенного решений в точке x = 0,5 [0, 1] будет равно:

Точность метода Ритца предлагается оценивать с помощью следующего практического приема [20].

В решаемой задаче используют функции.

Значения этих функций сравнивают между собой в нескольких точках интервала Если в пределах заданной точности решения задачи их значения совпадают, то в качестве результа та принимают функцию, описываемую выражением (3.7.8). В противном случае строится функция и срав ниваются значения функций и. Этот процесс продолжается до тех пор, пока значения и не совпадут в пределах заданной точности решаемой вариацион ной задачи.

Изложенный выше метод Ритца также может быть исполь зован при решении изопериметрических вариационных задач вида (2.4.14), (2.4.34).

Проиллюстрируем его применение на примере задачи с функционалом (2.4.4) и ограничением:

x F1 ( x, y, y ' ) = l1, (3.7.14) x которое является частным случаем выражений (2.4.34).

Для ее решения выберем приближенное представление экс тремали у(х) в форме зависимости (3.7.4).

Подставляя ее в функцию F 1 и вычисляя интеграл, входя щий в ограничение (3.7.14), получаем выражение вида:

(3.7.15) Искомые значения коэффициентов определя ются из решения задачи на условный экстремум целевой функ ции (3.7.6) при выполнении ограничений (3.7.15).

Данная задача решается методом Лагранжа, описываемого соотношениями (2.2.6), (2.2.7).

Для решения вариационных задач на условный экстремум вида (2.4.14), (2.4.28) и (2.4.14), (2.4.34) может быть применен метод малого параметра, предложенный в Разд. 3.4.

Суть такого подхода покажем на задаче из Примера 2.13.

Используя выражения (3.5.13), (3.5.14), заменим ограниче ние (x,y 1,y 2 ) = 0 сингулярным дифференциальным уравнением относительно множителя Лагранжа (х).

Тогда искомые функции у 1 (х), у2 (х), (х) решаемой задачи могут быть получены из решения следующей системы диффе ренциальных уравнений:

d y' = 0;

(3.7.16) y 2 dx 1 + y '1 + y ' с начальными условиями:

(3.7.17) Здесь µ – малый параметр;

x 0 – начальная точка процес са решения задачи Коши (3.7.16), (3.7.17).

Начальные значения в этой задаче вы бираются, исходя из того, что ее решения при х = х 0 должны удовлетворять условиям вида:

где – требуемая точность решения рассматриваемой вариаци онной задачи.

Для численного решения системы уравнений (3.7.16) пер вые два дифференциальных уравнения 2-го порядка преобра зуются с помощью замены:

к четырем уравнениям 1-го порядка.

Полученная система уравнений 5-го порядка решается од ним из методов, изложенных в Разд. 3.1.

3.8. Градиентные методы решения задач оптимального управления Как показала практика, применение принципа максимума Л.С. Понтрягина при решении реальных задач редко позволяет в получать вид оптимального управления аналитической форме. Основной причиной этого являются дос таточно сложные математические модели управляемого движе ния вида (2.5.1), (2.5.2), (2.5.5) и необходимость доопределения «свободных» параметров таких моделей.

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

Для иллюстрации содержания этого вида численных мето дов рассмотрим метод скорейшего спуска для решения задач оптимального управления со свободными границами [11]. Ис пользование термина «скорейший спуск» объясняется тем, что на каждом этапе «улучшения» управления каждое предыдущее управление корректируется в направлении наибыстрейшего приближения к минимальному (максимальному) значению функционала J, описываемого вектором его градиента [17].

Пусть требуется решить задачу выбора оптимального управления (2.5.1), (2.5.2), (2.5.5), (2.5.11) при следующих до пущениях:

1) моменты времени t 0 и t 1 являются заданными (фиксиро ванными) величинами, 2) начальные условия (2.5.2) определены для всех фазовых координат x1(t), x2(t), …, xn(t) объекта, входящих в модель (2.5.1), 3) функционал (2.5.11) зависит только от неизвестных («свободных») граничных значений фазовых координат объек та в момент времени t = t 1 и имеет вид:

J = f 0 ( x11, x21,, xn1 ) max, (3.8.1) где x j1 = x j (t 1 ),.

Зададим произвольное управление, удовлетворяющее ограничениям (2.5.5).

Поставим задачу формирования управления:

u r = u r + u r, r = (1, k ), (3.8.2) которому соответствует максимальное значение приращения функционала (3.8.1), вычисляемого как:

J = J J. (3.8.3) Здесь u r = ur (t) – искомое приращение функции u r (t) в момент времени t [t 0, t 1 ];

J* – значение функ ционала (3.8.1) при использовании управления вида (3.8.2).

Будем считать, что корректирующие поправки u r,, являются достаточно малыми по абсолютной вели чине. Это означает, что приращение J, описываемое выраже нием (3.8.3), достаточно точно представляется вариацией J функционала (3.8.1).

Управления, подставленные в уравнения, (2.5.1), после их решения при начальных условиях:

x j (t 0 ) = x j 0, j = (1, n ), (3.8.4) определяют траекторию x1 = x1 (t ), x2 = x2 (t ), …, xn = xn (t ) движения объекта на интервале времени [t 0, t 1 ].

В связи с малостью u r,, можно считать, что j = (1, n ), x = x j + x j, j где x j = x j (t) – непрерывные функции, достаточно малые по абсолютной величине в любой момент времени t [t 0, t 1 ],.


При этом:

x j (t0 ) = 0, j = (1, n ).

Учитывая это допущение, представим вариацию функцио нала (3.8.1) как:

n f J = 0 x j1, (3.8.5) j =1x j где x j1 = xj (t 1 ),.

Из доказательства принципа максимума Л.С. Понтрягина [11] будем использовать следующее соотношение:

t1 k n j1x j1 = u u r dt, (3.8.6) t0 r = j =1 r компоненты которого определяются выражениями (2.5.12) и (2.5.13).

Из второй группы условий (2.5.14) следует, что f j1 = 0, j = (1, n ). (3.8.7) x j Тогда с учетом этого равенства и выражений (3.8.5) и (3.8.6) получаем формулу вида:

t1 k J = u r dt, r =1 u r t отражающую непосредственную связь вариации функционала J с поправками управления u r,. Знаки этих поправок, которые должны доставлять максимум функционалу J, долж ны определяться из следующих соотношений:

u, r = (1, k ), (3.8.8) u r = µ, sign u r = sign u r r где – заданная константа, sign() – символ знака величины ().

Формируя достаточно малые значения поправок u r = ur (t), t [t 0, t 1 ] и конкретизируя их знаки из условий (3.8.8), находим «улучшенное» управление по формуле (3.8.2).

Так как управления не должны нарушать ограничений (2.5.5), то выбор поправок u r (t) необходимо проводить с уче том следующих неравенств:

u1r ur (t ) ur (t ) u2 r ur (t ), r = (1, k ). (3.8.9) Алгоритм метода включает в себя следующие этапы:

1°. Выбор шага t и формирование на интервале [t0, t 1 ] сет ки значений t s,.

2°. Задание в каждой точке сетки управлений u r (t s),,.

3°. Численное интегрирование с шагом t системы уравне ний (2.5.1) с начальными условиями (3.8.4).

4°. Вычисление значения функционала J и производных с помощью выражения (3.8.1) с использо, ванием полученных значений xj (t s ),,.

5°. Численное интегрирование системы уравнений (2.5.13) с начальными условиями (3.8.1) в направлении от момента вре мени t = t 1 к моменту t = t 0 с шагом t.

6°. Нахождение значений производных в, каждый момент времени t s, для вычисленных значе ний сопряженных функций j (t s ),.

7°. Формирование значения поправок u r (t s ),, и конкретизация их знаков с использованием условий (3.8.8) для выбранного значения параметра µ.

8°. Определение управлений по формуле (3.8.2).

9°. Проверка выполнения условия (3.8.9) для каждого мо мента времени t = t s,. При их нарушении – переход к п. 7° с уменьшением значения параметра µ.

10°.Полагая, переходим к п. 3° до тех пор, пока не будет выполнено условие |J|, где J вычисляется по фор муле (3.8.3), а величина определяет требуемую точность ре шения задачи.

Примечание.

Если на некотором шаге будет получено значение J 0, то осуществляется возврат к п. 7° с уменьшением поправок |u r |,.

В п. 5° производится интегрирование уравнений (2.5.13) в «обратном» времени, тогда как в п. 3° система (3.6.1) интегри руется в «прямом» времени от t = t0 к t = t 1.

Для выполнения общего интегрирования систем уравнений (2.5.1) и (2.5.13) в «прямом» времени функции j = j (t), предлагается вычислять по формулам:

j = (1, n ), n c ( ), j = j = где ( ) = ( ) (t ) – вспомогательные сопряженные функции, оп j j ределяемые из уравнений:

( ) ( ) = H j x j с начальными условиями ( ) (t ) =. j j Здесь функции () определяются с использованием правых частей уравнений (2.5.1) как:

= (1, n ).

n ( ) j, H ( ) = j j = Применяемый в начальных условиях параметр j (символ Кронекера) имеет следующие значения:

1, если = j, j = 0, если j, = (1, n ), j = (1, n ).

Постоянные c 1, c 2, …, c, …, c n определяются из решения системы линейных алгебраических уравнений вида:

f c ( ) (t1 ) = 0, j = (1, n ), n j x j = одним из численных методов, приведенных в Разд. 3.2.

Метод скорейшего спуска для численного решения задачи оптимального управления с фиксированными граничными ус ловиями вида (2.5.22) подробно изложен в работе [11]. Обзор существующих и перспективных численных методов решения задач оптимального управления с обширным списком исполь зованной литературы приведен в статье [53].

Следует отметить перспективность применения в сложных задачах оптимального управления БЛА генетических алгорит мов формирования искомого управляющего вектора, которые позволяют обеспечить достижение глобального максимума (минимума) используемых в них функционалов.

В заключение главы отметим, что все приведенные в ней численные методы прошли многолетнюю практическую апро бацию при решении сложных научно-технических задач, вклю чая задачи динамики полета и управления ЛА [52]. Эти методы в последующем реализуются в составе функционального про граммного обеспечения АРМ математика – системного про граммиста и используются персоналом БАК для решения задач формирования управлений БЛА при выполнении ими конкрет ных полетных заданий.

Глава 4. ВЕТРОВЫЕ ВОЗМУЩЕНИЯ ТРАЕКТОРИЙ ДВИЖЕНИЯ БЛА Необходимость учета ветровых возмущений при формиро вании программного управления БЛА определяется следую щими факторами, отмеченными в Главе 1:

1) БЛА имеют малую полетную массу и небольшие скоро сти полета по сравнению с современными пилотируемыми ЛА, 2) одним из преимуществ БЛА по сравнению с такими ЛА является возможность их применения в ночное время и при не благоприятных метеоусловиях, в частности при наличии ветра.

Как показала практика эксплуатации пилотируемых ЛА, действие ветра вызывает изменение величины перегрузок, а так же линейные колебание его центра масс (ЦМ) [48]. Последнее весьма отрицательно сказывается в процессе эксплуатации БЛА РМ на результаты решения задачи анализа оператором целевой нагрузки получаемого изображения наземной поверхности.

Анализ литературы по БЛА [1, 2-4, 7, 15, 31] показал пол ное отсутствие рассмотрения вопросов влияния ветровых воз мущений на параметры и характеристики их полетов, и, глав ное, на выбор соответствующих управлений, позволяющих эф фективно выполнять поставленные перед БЛА целевые задачи.

Учет влияния действия ветра на различных частных режи мах полета пилотируемых ЛА рассматривался в работах [13, 19, 23, 27-30, 39, 48]. Основным недостатком этих работ является отсутствие комплексного подхода, как к описанию произволь ных ветровых возмущений, так и к их влиянию на совокуп ность параметров и характеристик полета ЛА.

В связи с приведенными выше особенностями беспилотной авиационной техники учет ветровых возмущений является од ной из важных задач прикладной теории управления БЛА.

4.1. Общая характеристика ветровых возмущений Следует заметить, что в литературе [23, 30, 39, 40, 48], по священной достаточно подробному изложению этих вопросов, отсутствует единый подход к определению понятия ветра и описанию его характеристик.

В работе [48] отмечается, что на силу тяги и аэродинамиче ские силы существенное влияние оказывает распределение по траектории движения ЛА значений таких метеорологических параметров окружающей среды как ветер, давление, темпера тура и влажность воздуха. При этом наибольшее прямое воз действие на динамику полета ЛА оказывает ветер, под которым в общем случае понимается любое движение воздуха относи тельно земной поверхности.

Воздушная среда, в которой осуществляется полет ЛА, практически всегда находится в непрерывном движении, вызы ваемом неравномерным распределением температур, плотно стей и давлений воздуха в атмосфере. Вследствие этого движе ние воздуха в общем случае носит хаотический случайный ха рактер.

Как показала практика, некоторые составляющие этих по токов имеют постоянную скорость. Такие потоки воздуха по лучили название струйных течений.

Струйные течения распространяются горизонтально земной поверхности и имеют четкий максимум скорости на определен ной высоте. При этом отмечается, что средние направления та ких течений в вертикальной плоскости практически совпадают в каждой точке с горизонтальным направлением. Длина и ши рина струйных течений составляет величины порядка несколь ких сотен километров [48].

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

В работе [19] под турбулентной атмосферой понимается со стояние неспокойной атмосферы при наличии в ней вертикаль ных потоков воздуха. Турбулентность воздушных масс на ма лых высотах в основном возникает в теплое время года вслед ствие неравномерного нагрева земной поверхности (пашни, ле са, водоемы и т.п.). На средних высотах турбулентность появ ляется на границах холодных и теплых фронтов воздуха, а так же в кучевых и мощнокучевых облаках. На больших высотах (h = 11000-13000 м) наблюдаются горизонтальные течения воз душных масс с различными распределениями скоростей по вы соте. Здесь при большом перепаде скоростей образуется значи тельная турбулентность, вызывающая болтанку самолета. В этой работе отмечается, что скорость восходящих потоков в турбулентной атмосфере у земли, на средних и больших высо тах достигает 10-17 м/с, а в грозовом фронте – 25-47 м/с. Ско рости нисходящих потоков воздуха в такой атмосфере обычно несколько меньше скоростей восходящих потоков.

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

Следует отметить, что в настоящее время полные характе ристики случайного поля скоростей ветра для любой точки земной атмосферы отсутствуют. Поэтому в практических рас четах пользуются понятием средней скорости ветра.

Современные ЛА имеют эксплуатационные высоты полета до 10-15 км. С точки зрения характера поля скоростей ветра этот диапазон высот разбивается на два участка [40]:

1) пограничный слой до высот 1-10 км, 2) свободная атмосфера с высотами более 10 км.

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

В свободной атмосфере значения вертикальной составляю щей ветра обычно малы и его направление более стабильное.

Для струйных течений максимальная средняя скорость вет ра порядка 90 м/с. достигается на высоте около 10 км. Типовой интервал изменения этой скорости равен 30-50 м/с. с высотой струйных течений 5-17 км [48]. При этом отмечается, что в зимних условиях эта скорость увеличивается.

В приземном слое скорость ветра, начиная с высот 100 200 м. резко уменьшается с высотой. Для вычисления средней скорости ветра до высот h = 300-500 м. предлагается эмпириче ская формула вида:

n h wср (h) = w0. (4.1.1) h Здесь w 0 – скорость ветра на базовой высоте h 0 (например, h 0 = 10 м.);

n – показатель градиента ветра по высоте, который обычно имеет значения 0,2-0,5, но при значительных темпера турных изменениях может достигать значений 0,7 и выше.

Величина средней скорости ветра зависит от времени су ток. Так при h 50 м. ночью эта скорость выше, чем в дневные часы. Установлено, что максимальное значение такой скорости в Северо-западном регионе России составляет величину поряд ка 8-10 м/с. [48]. В этой же работе отмечается, что длина разбе га самолета среднего класса по ВВП при встречном ветре ско рости 10 м/с. почти в 2 раза меньше, чем при таком же попут ном ветре.

Наиболее часто в пограничном слое атмосферы возникают турбулентные течения. Установлено, что в таких течениях вих ревые порывы ветра не имеют предпочтительного направления при размерах вихрей от нескольких сантиметров до нескольких сот метров [48].

Процессы возникновения, развития и завершения турбулент ности являются достаточно длительными и при высокой скоро сти ЛА могут рассматриваться как установившиеся процессы.

Установлено [48], что толщина зон турбулентности менее 1000 м. встречается в северных широтах в 90% случаев, в уме ренных – в 85% случаев, в южных – в 70% случаев. Макси мальная толщина таких зон – не более 2000 м. Длина турбу лентных течений в 20% случаев менее 20 км, а в 70% случаев – менее 100 км. В зонах турбулентности вертикальные и горизон тальные порывы ветра длиной несколько километров могут достигать значений ±(20-40) м/с.

В работе [28] приводятся следующие количественные ха рактеристики ветровых возмущений.

Горизонтальная составляющая ветра может достигать сред ней скорости 12-18 м/с на умеренных высотах и до 20-30 м/с в стратосфере. В струйных течениях эта скорость может равнять ся 50-80 м/с. На малых высотах изменение горизонтальной со ставляющей скорости ветра по высоте имеет значение порядка 0,2-0,3 м/с на метр.

Вертикальная составляющая скорости ветра при отдельных порывах может достигать значений 18-20 м/с.

Среднее значения такой составляющей в турбулентных по токах воздуха не превышает обычно 3-5 м/с при длине волны 100-500 м.

В качестве общей количественной характеристики дейст вующего ветра будем рассматривать вектор его скорости (w x,wy,w z ). Компоненты этого вектора описывают состав ляющие скорости ветра, действующие в направлении соответ ствующих осей земной СК.

Непосредственное значение скорости действующего ветра вычисляется при известных значениях его компонент по фор муле вида [17]:

W = wx + w 2 + w z.

2 (4.1.2) y Как было отмечено выше, для струйных течений воздуха считается, что величина W является постоянной и равной сред ней скорости действующего ветра.

В этом случае имеем, что:

w x = w x ср = const;

wy = w y ср = const;

w z = wz ср = const,(4.1.3) Турбулентные течения воздуха предлагается также описы вать вектором, который представляется как сумма двух век торов [48]:

, (4.1.4) где – вектор постоянной составляющей средней скорости действующего ветра;

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

Компоненты вектора будем описывать выражениями вида (4.1.3). Как отмечается в работе [48], модуль этого векто ра, то есть среднее значение постоянной составляющей скоро сти турбулентного течения воздуха обычно имеет порядок не скольких метров в секунду. Для описания турбулентных дви жений воздуха используются следующие подходы [48]:

1) метод дискретных порывов, 2) метод непрерывных случайных порывов.

Первый метод (Рис. 4.1) применяется для анализа свойств ЛА и предполагает, что форма порыва ветра фиксирована. При этом ЛА перемещается под его воздействием вертикально без изменения угла тангажа, и динамические свойства всех ЛА, по павших в порыв ветра, одинаковы.

ЛА ЛА Рис. 4. Основной недостаток этого метода – неизвестная форма ре ального порыва ветра. Поэтому предлагается использовать второй метод, который применяется при следующих предположениях:

а) поле скоростей ветра на определенных участках турбу лентной атмосферы является однородным и изотропным, б) для движущегося ЛА поле скоростей ветра является «за мороженным», то есть не изменяющимся во времени (гипотеза Тейлора).

Однородность поля скоростей ветра означает, что его веро ятностные характеристики одинаковы для всего поля.

В работе [48] на основе полученных в экспериментальных полетах данных утверждается, что случайные значения компо распределены по нормальному закону с прак нент вектора тически нулевыми математическими ожиданиями и средне квадратическими отклонениями (СКО) wx 0, wy 0, wz 0.

При этом установлено, что w x ср w у ср и wx wz. Среднеквад ратические отклонения с уменьшением высоты полета с 500 м.

до 330 м. увеличиваются с 1,186 м/с. до 2,20 м/с. На трассе по лета 15 км величины этих отклонений изменяются не более чем на 25%. Для грозовых условий погоды значения СКО равны (4,5-5,0) м/с. [48].

Для описания в рамках данного метода характеристик тур булентной атмосферы будем использовать известное в теории вероятности правило «трех сигм» [17]. Это правило утверждает, что реализации нормально распределенной случайной величины Х с вероятностью не менее 0,995 лежат в интервале [m x – 3 x ;

m x + 3 x ], где m x и x – соответственно математическое ожида ние и СКО случайной величины Х.

Для получения гарантированных результатов при анализе влияния на движение БЛА турбулентных ветровых возмущений компоненты вектора (4.1.4) предлагается вычислять по сле дующим формулам:

w x = w x пост + 3 wx ;

w y = w y пост + 3 wy ;

wz = w z пост + 3 wz. (4.1.5) Классификация видов ветровых возмущений, действующих при движении БЛА, приведена на Рис. 4.2.

Отметим, что при учете первого вида возмущений исполь зуются выражения вида (4.1.3). Второй вид возмущений описы вается формулами (4.1.5). Данные, входящие в эти соотноше ния для различных районов полетов БЛА должны храниться в базе данных АРМ метеоролога БАЭ (см. Рис. 1.7).

Ветровые возмущения, действующие на БЛА Струйные ветровые Турбулентные вет возмущения ровые возмущения Рис. 4. На наш взгляд, наиболее точно учет ветровых возмущений проводится при использовании измеренных перед полетами БЛА компонент вектора скорости ветра, зависящих от времени.

Для района старта (взлета) и посадки БЛА метод расчета этих компонент иллюстрирует Рис. 4.3.

На этом рисунке представлена плоскость 0ст x ст z ст стартовой СК (см. Рис. 1.9). Угол, являющийся углом азимута ветра [82], определяет направление действия горизонтальной состав ляющей (w x,w z) пространственного вектора, измерен ной с помощью соответствующей аппаратуры.

zст wz 0ст wx xст Рис. 4. Обозначим через W xz (t) и (t) значения скорости этой со ставляющей и направления ее действия, измеренные в момент времени t.

Тогда изменение во времени значений соответствующих компонент вектора в районе старта (взлета) и посадки БЛА вычисляются по формуле вида:

w x (t) = w xz (t) cos (t);

wz (t) = w xz (t) sin (t). (4.1.6) Распределение значений wх (t) и w z (t) по высоте может быть получено с использованием выражения (4.1.1). При этом верти кальная компонента wу (t) скорости ветра измеряется отдельно.

Использование такого подхода предполагает наличие авто матической измерительной аппаратуры, устанавливаемой в рассматриваемом районе, данные от которой поступают по ра диоканалам в АРМ метеоролога БАЭ для их анализа и после дующей обработки.

Для определения зависимостей w x (t), w y (t), w z (t) компонент вектора скорости ветра, действующего в районах выполне ния БЛА полетных заданий, предлагается использовать специ альные метеорологические БЛА (МБЛА).

Измерения параметров полета таких БЛА позволяют вы числить параметры действующего ветра по определенным рас четным формулам, которые в общем случае имеют вид:

w x (t) = f 1 (p1 (t), p2 (t),…, p k (t));

w y (t) = f 2 (p1 (t), p2 (t),…, p k (t));

(4.1.7) w z (t) = f 3 (p1 (t), p 2(t),…, p k (t)), где p i (t) – фактические значения кинематических и динамиче ских параметров полета МБЛА в момент времени t,.

Отметим, что для достоверного измерения этих параметров МБЛА должен быть оснащен дополнительными датчиками.

В работе [48] приводятся следующие зависимости вида (4.1.7), использованные в летных экспериментах по измерению скорости ветра:

wx (t ) = (V Vзад ) g (nx )dt ;

wy (t ) = V (ф ) + g (n y 1)dt + l x ;

wx (t ) = V ( ф ) + g (nz + )dt + l y l x.



Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 14 |
 





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

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