Scientific journal
Scientific Review. Technical science
ISSN 2500-0799
ПИ №ФС77-57440

IDENTIFICATION OF EXTREMUMS OF THE NORMS OF SOLUTIONS OF ORDINARY DIFFERENTIAL EQUATIONS WITH AN APPENDIX TO THE ESTIMATION OF STABILITY

Zaika I.V. 1
1 Taganrog Institute named after A.P. Chekhov (branch) Rostov State University of Economics
A computer method of numerical optimization based on sorting is proposed, which helps identify extremums and zeros of functions by values and indexes of data. The presented scheme allows us to programmatically calculate the extrema of an arbitrary function of the form in question in the domain of the function definition and can be transferred to the case of a function of three or more variables, but this requires an increase in the expenditure of computer time. The method differs from the known by the parallelization possibility, and without fundamental changes it is transferred to the case of solving systems of nonlinear equations of general form. In this case, the numerical optimization algorithm is invariant with respect to the dimensionality of the system and the type of problem. The method is used to search for extrema and zeros of difference solutions of systems of ordinary differential equations (ODE). The presented algorithm is applied to the search for extrema of the norms of the ODE solutions. Based on the optimization algorithm and computer stability estimation schemes, it is possible to estimate locally and globally the extreme deviation of the system from a stable state with variation of parameters. The work includes the texts of programs and the description of a numerical experiment, in particular, the application of the scheme to the analysis of perturbations of high-power power systems under perturbation of parameters is considered.
identification of extremums of solutions of differential equations
еstimation of the stability of difference solutions of differential equations

Предлагаются алгоритмы для определения по индексам входных данных экстремальных значений норм решений систем ОДУ. Схема модифицируется для поиска экстремумов норм [1, 2] преобразованных разностных решений систем уравнений при вариации параметров. Анализируется экстремальное отклонение решения системы от устойчивого состояния, и на основе этих данных алгоритмы применимы к оценке устойчивости, в частности при возмущении параметров.

Вычисление экстремальных значений функции двух действительных переменных

В начале рассматривается функция одной действительной переменной y = f(x). Пусть функция определена и задана на промежутке zaik01.wmf. На сетке считываются значения функции f(x), zaik02.wmf, zaik03.wmf.

Значения c[i] сортируются. Процедура сортировки в программе обозначена как sort [3, 4]. К выходу процедуры сортировки подсоединяется специально сконструированный оператор: zaik04.wmf,
zaik05.wmf. С помощью данного оператора производится автоматическая идентификация каждого минимума в промежутке ε среди элементов c[i].

Смысл условия в том, что в ε-промежутках входных значений функции с индексом e[k] нет значений уже в отсортированном массиве c1, больших по значению, чем элемент с индексом [5]. Это значит, что значение функции y в узле zaik07.wmf не превосходит значений во всех точках заданной ε-окрестности.

Если оператор zaik08.wmf
заменить на другой оператор, с помощью которого локализуется максимум zaik09.wmf, zaik10.wmf, то осуществляется вычисление максимумов функции.

Изложенная схема применяется для вычисления всех экстремумов функции двух переменных zaik11.wmf [6]. Задаются границы zaik12.wmf и zaik13.wmf,
где строится область: zaik14.wmf,
zaik15.wmf, zaik16.wmf, zaik17.wmf, zaik18.wmf.

Для нахождения минимумов функции zaik19.wmf перебираются значения вдоль j-го столбца по ординате OY, в результате чего находится наименьшее по строкам значение zaik20.wmf, zaik21.wmf, которое поступает на вход сортировки с индексом j.

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

Оператор для текущего узла e[k] находит каждый узел zaik22.wmf, в проекции ε, на абсциссу OX среди которых нет точек, предшествующим значениям в отсортированном массиве.

Значение точки минимума zaik23.wmf фиксируется, после чего локализуется ордината, в которой вычисляется приближенное минимальное значение функции двух действительных переменных.

Представленная схема позволяет программно вычислять экстремумы произвольной функции рассматриваемого вида. Алгоритм может использоваться для идентификации экстремумов функций нескольких переменных [7].

Решение нелинейных систем уравнений

Схема применяется к приближенному решению систем нелинейных уравнений. Пусть исходная система преобразована к виду zaik24.wmf, где левые части являются функциями от n действительных переменных. Совокупность функций zaik25.wmf образует n-мерную вектор-функцию аргументов zaik26.wmf. Пусть требуется приближенно найти все решения системы zaik27.wmf в заданной области zaik28.wmf, zaik29.wmf, …,
zaik30.wmf. Строится сетка zaik31.wmf, zaik32.wmf, zaik33.wmf, zaik34.wmf. Во всех ее узлах вычисляются значения нормы zaik35.wmf рассматриваемой функции, которые принимаются за элементы n-мерного массива: zaik36.wmf,
где индексы соответствуют узлам многомерной сетки. Левая часть zaik37.wmf интерпретируется как дискретная функция n переменных, ее минимумы можно идентифицировать программно по изложенной схеме.

Пример. Приближенное решение системы

zaik38.wmf

в области zaik39.wmf, zaik40.wmf, zaik41.wmf вычисляется по программе Program UravnNelin, представленной в листинге [2]:

label 21,22,23; const eps0=0.1; h=eps0/23; x00=0.1;

x11=1; y00=0; y11=1; z00=0; z11=1;

n00=trunc(6/h); nn=n00+round(n00/2)+1;

type vect1=array [1..2*nn+nn] of double;

vect2=array [1..2*nn+nn] of longint;

var i,j,k,k1,r,rx,ry,rz,nn0,kx,ky,kz : integer;

c,a1,az: vect1; e, ex, ey,ez: vect2;

x,x0,x1,xk,minx,minz: double;

y,y0,y1,yk,z0,z1,z,zk: double;

function ff1 (x,y,z: double): double;

begin ff1:= x*x+y*y-z ; end;

function ff2 (x,y,z: double): double;

begin ff2:= x*x+y*y-z*z; end;

function ff3 (x,y,z:double): double ;

begin ff3:= ln(x)-sqrt(y)+0.8; end;

function func (x,y,z: double): extended;

var p: double; i1: integer;begin

func:= abs(ff1(x,y,z))+abs(ff2(x,y,z))+abs(ff3(x,y,z)); end;

procedure minyz (var x,y,z,minx: double);

begin minx:=func(x,y,z); i:=1;

while i<= trunc((y11-y00)/h) do begin

j:=1; while j<= trunc((z11-z00)/h) do

begin y:=y0+i*h; z:=z0+j*h;

IF minx > func(x,y,z) then minx:=func(x,y,z);

j:=j+1; end; i:=i+1; end; end;

procedure minyy (var x,y,z,minz:extended);

begin minz:=func(x,y,z); i:=1;

while i<=trunc((y11-y00)/h) do begin

begin y:=y0+i*h; if minz > func(x,y,z) then

minz:=func(x,y,z) end; i:=i+1;end; end;

procedure sort(var nn0:

longint; var c: vect1; var e: vect2);

{ Oписание процедуры sort }

begin x0:=x00; x1:=x11;

y0:=y00; y1:=y11;z0:=z00; z1:=z11; nn0:=n00;

for rx:=1 to nn0 do begin

x:=x0+rx*h; minyz (x,y,z,minx); a1[rx]:=minx; end;

sort( nn0, a1, ex); kx:=1;

while kx<= nn0 do begin

for rx := 1 to kx-1 do

IF abs( ex[kx]-ex[kx-rx] ) <= eps0/h then goto 21;

xk:= x0+ex[kx]*h; for rz:=1 to nn0 do begin

z:=z0+rz*h; x:=xk;

minyy (x,y,z,minz); az[rz]:=minz; end;

sort( nn0, az, ez); kz:=1;

while kz<= nn0 do begin

for rz := 1 to kz-1 do

IF abs(ez[kz]-ez[kz-rz]) <=eps0/h then goto 22;

zk:= z0+ez[kz]*h;

for ry:=1 to nn0 do begin

y:=y0+ry*h; a1[ry]:=func(xk,y,zk) end;

sort( nn0, a1, ey); ky:=1;

while ky<= nn0 do begin

for ry := 1 to ky-1 do

IF abs(ey[ky]-ey[ky-ry]) <=eps0/h then goto 23;

yk:= y0+ey[ky]*h;

writeln (' ', xk,' ');

writeln (' ', yk,' ');

writeln ( zk,' ', ff1(xk,yk,zk):3:5,' ',ff2(xk,yk,zk):3:5,'',ff3(xk,yk,zk):3:5);

23: ky:=ky+1; end;

22: kz:=kz+1; end;

21: kx:=kx+1; end; end.

Результат вычислений:

Значения аргументов

Значения компонент

X

Y

z

0,8870

0,4609

1,0000

f1 = –0,0009

f2 = –0,0009

f3 = –0,0012

 

Представленный алгоритм может использоваться для решения трансцендентных уравнений [2] и по способу построения отличается от аналогов [8–10].

Идентификация экстремумов и нулей разностных решений систем ОДУ

Пусть рассматривается задача Коши для уравнения

zaik42.wmf (1)

и выполняются условия существования и единственности решения. Для решения (1) с малым шагом h на промежутке zaik43.wmf используется метод Эйлера:

zaik44.wmf (2)

Разностные значения из (2) рассматриваются как элементы массива

zaik45.wmf, zaik46.wmf. (3)

После сортировки массива (3) применяется оператор локализации экстремумов. Схема локализации экстремальных элементов в массиве (3) такая же, как и для массива (1) функции одной переменной, за исключением того, что в случае разностной схемы нельзя выполнить спуск. Погрешность вычисления экстремумов решения уравнения (1) полностью определяется погрешностью разностного приближения. С точностью до этой оговорки имеют место устойчивость локализации и вычисления экстремумов решения (2) [11, 12].

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

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

zaik47.wmf (4)

где

zaik48.wmf,

zaik49.wmf,

zaik50.wmf,

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

zaik52.wmf

zaik53.wmf zaik54.wmf

zaik55.wmf zaik56.wmf.

На выходе метода идентифицируются все локальные экстремумы каждой переменной yk, в виде последовательности ck[i], в окрестностях малого промежутка на всем отрезке zaik57.wmf.

Нули разностных приближений компонент получатся, если на вход подавать модули значений zaik58.wmf с оценкой на выходе малости минимумов. По-прежнему необходимо смещение текущего отрезка по всему промежутку поиска экстремумов с проверкой границ смещаемого отрезка на экстремум. Вместо метода Эйлера могут применяться разностные методы высшего порядка [14], численный эксперимент показывает возможность идентификации экстремумов с точностью до формата представления данных в языке программирования.

Под устойчивостью понимается устойчивость по Ляпунову [1]. Пусть рассматривается система нелинейных ОДУ:

zaik59.wmf

где zaik60.wmf определяются по методу Эйлера, zaik61.wmf; zaik62.wmf – варьируемые числовые параметры в диапазоне zaik63.wmf, zaik64.wmf, zaik65.wmf.
Возмущенные начальные данные обозначаются zaik66.wmf соответственное им возмущенное решение записывается в виде zaik67.wmf. Требуется найти все экстремальные отклонения от нуля нормы разности между возмущенным и невозмущенным решениями системы dY/dt при вариации числовых параметров.

Способ оценки отклонений опирается на схему идентификации экстремумов [4] дискретно представленных функций четырех действительных переменных. На вход алгоритма подается функция одной независимой переменной t, роль трех других играют варьируемые параметры. Роль функции играет норма разности вычисляемых по разностной схеме значений вектор-функций zaik68.wmf. При выборе нормы:

zaik69.wmf,

zaik70.wmf.

На вход метода поступает

zaik71.wmf,

zaik72.wmf.

При этом решения zaik73.wmf вычисляются по разностной схеме отдельно для каждого набора дискретизированных значений трех варьируемых параметров. Выбранные дискретизированные значения левой части c[i] задают трехмерный массив, к которому добавляется еще одно измерение по независимой переменной t. Получится четырехмерный массив с элементами zaik74.wmf, где zaik75.wmfиз zaik76.wmf, i = 1, 2, ..., N
при значениях параметров zaik84.wmf,
zaik85.wmf, индексы которых указывают номера шагов дискретизации.

Найденные экстремумы характеризуют меру отклонения возмущенного решения от невозмущенного (меру возмущения решения). Значение глобального максимума позволит найти наибольшее значение возмущения на отрезке zaik79.wmf при всех дискретных значениях трех параметров.

Соответственные программы приводятся в [2, 6]. В [2] детально изложено применение подхода к анализу возмущений энергетических систем при возмущении параметров. Изложенный метод применяется к оценке устойчивости на основе анализа поиска нулей характеристического полинома. Для системы [2] характеристическое уравнение имеет вид zaik80.wmf или zaik81.wmf.

Нули полинома и трансцендентной функции в левой части последнего уравнения можно найти как минимум модуля функции zaik82.wmf при вариации параметра zaik83.wmf. Данный метод переносится на поиск решений систем нелинейных уравнений, которые могут содержать трансцендентные функции, при этом вместо минимума модуля ищется минимум нормы компонент зависимых переменных на многомерной сетке.

Заключение

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

Работа выполнена при финансовой поддержке РФФИ (проект 16-07-00100).