В файле:
-
транспонированный двумерный массив ,
и соответственно
равны и ,
и соответственно
равны и ,
- размер
опоры ограничений начального согласованного опорного плана,
-
обратная матрица к опорной матрице,
- размер
опоры целевой функции начального согласованного опорного плана,
-
вектор-столбец, содержащий элементы верхнего треугольника матрицы, обратной к
опорной матрице целевой функции, записанные по строчкам,
-
вектор-столбец размерности , который
содержит следующую информацию: индексы опорных компонент начального плана,
вошедшие в опору ограничений (их количество);
индексы опорных компонент начального плана, вошедшие в опору целевой функции
(их количество );
индексы неопорных компонент начального плана (их количество ).
KOD(2.4) - код
завершения подпрограммы, принимает значение:
KOD(2.4)=0 -
полученный оптимальный план задачи,
KOD(2.4)=1 -
исчерпан ресурс итераций,
KOD(2.4)=2 -
получено направление бесконечного убывания целевой функции (массив R),
KOD(2.4)=3 -
вычислительный процесс остановлен из-за накопления ошибок округлений (точка ,
содержащаяся в массиве X, не удовлетворяет ограничениям задачи),
KOD(2.4)=4 -
вычислительный процесс остановлен из-за недостаточных размеров массива E (E -
двумерный массив содержит при в первых
строках
и первых столбцах
опорную матрицу), так как количество индексов, вошедших в опору ограничений,
превысило значение первого измерения, двухмерного массива E.
Выходной файл содержит следующую результирующую информацию:
значение параметра KOD(2.4).
Если он равен нулю, то построен оптимальный план задачи;
размерность опоры ограничений;
опору ограничений;
размерность опоры целевой функции;
опору целевой функции;
неопорные индексы;
оптимальный план задачи;
6. Примеры
Пример 1. Рассмотрим систему управления
(6.1)
Эта
система при выключенном управлении имеет
периодические движения для любых начальных состояний. Однако эти решения не
являются предельными циклами. Выделим из них одно движение
(6.2)
и
построим обратную связь (3. 3), с которой замкнутая система
(6.3)
имеет периодическое движение (6.2) в качестве асимптотически устойчивого
предельного цикла.
При построении обратной связи были выбраны следующие параметры
вспомогательной задачи оптимального управления (4.1):
Результаты решения задачи:
Рисунок 1 и 2 приведены для движения с началом в точке (-2,0)
Рисунок 1-Реализация управления
На
рисунке видно, что фунция стабилизируется
в 0.
Рисунок
2- Изменение функции Ляпунова
Рисунок
3 - Фазовые траектории замкнутой
системы.
С
началом движения в точках:
1
- (-2,0);
2
- (-0.5,0);
3
- (0.5,0);
4
- (2,0);
5
- траектория предельного цикла.
Пример 2. Рассмотрим систему управления
(6.4)
которая
при имеет
периодические движения, среди которых нет заданного предельного цикла
(6.5)
С помощью решения вспомогательной задачи (4.1) при
(управление, осуществляющее заданное движение) была построена
ограниченная обратная связь, после замыкания которой заданное движение (6.5)
стало асимптотически устойчивым предельным циклом.
Результаты решения задачи:
Рисунок 1 и 2 приведены для движения с началом в точке (-2,0)
Рисунок 4-Реализация управления
На
рисунке видно, что фунция стабилизируется
в -2.
Рисунок 5- Изменение функции Ляпунова
Рисунок 6 - Фазовые траектории замкнутой системы.
С началом движения в точках:
1 - (2,0);
3 - (2,2.5);
4 - (2,4);
5 - траектория предельного цикла.
Заключение
В работе показано, что методами оптимального управления с привлечением
современных средств вычислительной техники можно эффективно решить неизменно
актуальную задачу следящих систем. Были сформулированы понятия ограниченной
обратной связи и дискретного управления, а также доказаны свойства стартовой
обратной связи. Задача следящих систем была сведена к линейно-квадратичной
задаче оптимального управления. Для решения данной задачи был построен алгоритм
обратной связи и выбран двойственный метод. Приведена программа на языке С для
решения.
Список использованных источников
1 Айзерман, М.А. Лекции по теории автоматического
регулирования. М.: Физматгиз, 1958.
2 Габасов, Р., Кириллова, Ф.М., Ружицкая, Е.А. Решение
классической задачи регулирования методами оптимального управления // АиТ.
2001. М. С. 18-29.
3 Барбашип, Е.А. Введение в теорию устойчивости. М.: Наука,
1967.
Квакерпаак X., Сиван Р. Линейные оптимальные системы
управления.М.: Мир, 1977.
Малкин И.Г. Теория устойчивости движения. М.: Наука, 1966.
6 Поптрягин, Л.С., Болтянский, В.Г., Гамкрелидзе, Р.В.,
Мищенко, Е.Ф. Математическая теория оптимальных процессов. М.: Наука, 1969.
7 Габасов, Р., Кириллова, Ф.М., Костюкова, О.И. Построение
оптимальных управлений типа обратной связи в линейной задаче // Доклады АН
СССР, 1991. Т.320.
Ж.С.1294-1299.
8 Gabasov, R., Kirillova, F.M., Ruzhitskaya, Е.А. Stabilization of dynamical systems with the help of optimization
methods // Nonsmooth and discontinuous problems of control and optimization.
Proceedings volume from the IFAC Workshop, Chelyabinsk, Russia, 1998.P.35-41.
9 Габасов,
P., Кириллова, Ф.М., Тятюшкин, А.И. Конструктивныеметодыоптимизации. Часть 1. Линейные задачи. Минск:
Изд-во “Университетское”, 1984.
10 Данциг, Дж. Линейное программирование, его применения и
обобщения. М.: Прогресс, 1966.
11 Габасов, Р., Кириллова, Ф.М., Ружицкая, Е.А. Гашение
колебаний струны // ДНАН Беларуси. 1999. Т.43, №5, С.9-12.
12 Габасов, Р. Конструктивные методы оптимизации: в 5 ч. Ч.
1. Линейные задачи / Р. Габасов, Ф. М. Кириллова, А. И. Тятюшкин. - Мн.: БГУ,
1983. - 214 с.
Габасов, Р. Конструктивные методы оптимизации: в 5 ч. Ч. 2.
Задачи управления / Р. Габасов, Ф. М. Кириллова. - Мн.: Университетское, 1984.
- 204 с.
Конструктивные методы оптимизации: в 5 ч. Ч. 4. Выпуклые
задачи / Р. Габасов [и др.]. - Мн.: Университетское, 1987. - 223 c.
Ракецкий, В.М. Решение общей задачи выпуклого квадратичного
программирования прямым методом / В. М. Ракецкий // Програм. обеспечение ЭВМ. -
Мн.: Ин-т математики АН БССР, 1985. - Вып. 55. - С. 113-123.
Ракецкий, В.М. Решение общей задачи выпуклого квадратичного
программирования двойственным методом / В. М. Ракецкий // Програм. обеспечение
ЭВМ. - Мн.: Ин-т математики АН БССР, 1985. - Вып. 55. - С. 124-129.
Габасов, Р. Построение оптимальных управлений типа обратной
связи в линейной задаче / Р. Габасов, Ф. М. Кириллова, О. И. Костюкова // Докл.
АН СССР. - 1991. - Т. 320,№ 6. - С. 1294-1299.
Габасов, Р. К методам стабилизации динамических систем / Р.
Габасов, Ф. М. Кириллова, О. И. Костюкова // Изв. РАН. Техн. кибернетика. -
1994. -№ 3. - С. 67-77.
Габасов, Р. Стабилизация линейных динамических систем
оптимальными управлениями линейно-квадратичных задач / Р. Габасов, А. В.
Лубочкин // Прикл. матем. и механ. - 1998. - Т. 62, Вып. 4. - С. 556-565.
Лубочкин, А.В. Использование ограниченных оптимальных
обратных связей линейно-квадратичных задач для осуществления заданных движений
динамических систем / А.В. Лубочкин // IX Белорусская матем. конф.: тез. докл.: в 3 ч. Ч. 3
(2004;Гродно), 3-6 ноября 2004 г., г. Гродно. - Гродно: ГрГУ, 2004. - С.
124-125.
Приложения
Приложение А
Текст программы 1
#pragma hdrstop
# include <stdio.h>
# include <math.h>
# include <string.h>
# include <process.h>
# include <alloc.h>
# include <conio.h>
#define M 2
#define N 25
#define N_M 23
#define TBEG0.0
#define TETA1.0
#define L 4.0
#define NX 2tau,[NX] = { 0.5, 0.0 },0[NX] = { 1.0, 0.0 },
z[NX],// Точка на пред.цикле y=y(t):
// z(t) = y(t+TETA) = F(t+TETA)*z0[M],A[M][N],dn[N], dw[N],
eps1 = 0.00001, eps2=0.00001,[M][M], **Pr, Q[M][M], **D,v[N_M], **Ma,u[N],
Del[N],J, a1, a2;Jop[M],[N_M], kn,[N_M], koc,[N_M], knn,;callocPr ( int
);freePr ( int );puN1NIntU2 ( double );puIntU2 ( double );formParam ( double
);formElem ( void );procsogl( void );gs ( int, double**, /* Метод Гаусса */*,
int* );f_x ( double, double*,/* f1(t) = x2(t) */* ); /* f2(t) = -x1(t) + u(t)
*/FEHL(int, double,*, void (f_p)(double, double*, double*),*, double*,*,
double*,*, double*,*, double*);RKF45 (int, double, /* Метод
Рунге-Кутта-Фельдберга */, double*,(f_p)(double, double*, double*), double*,*,
double*,*, double*,*, double*,*, int*,*, int*,*, int*,*, double*,*,
double*);prifl0 (int*, double*, double*);*d;*r;*ru; /* Файл.dat : t, u(t)
*/*rJ; /* Файл.dat : t, J(t) */*rx; /* Файл.dat : t, x(t) */
#pragma argsusedmain(int argc, char* argv[])
{f_Name_u[13],/* Имя файла.dat : t, u(t) */_Name_J[13],/* Имя
файла.dat : t, J(t) */_Name_x[13];/* Имя файла.dat : t, x(t) */h, t, tout, ypx
[NX], hx=0, f1x [NX], f2x [NX], f3x [NX], f4x [NX], f5x [NX], rerr = 1.0e-12,
aerr = 0.0,= 0, savaex = 0;j, Ostanov, inte, iflx, jflx, kflx, nfex, kopx,
initx;= (double**)calloc( N_M, sizeof(double) );( j=0; j<N_M; j++ )[j] =
(double*)calloc( N_M, sizeof(double) );= (double**)calloc( N_M, sizeof(double)
);( j=0; j<N_M; j++ )[j] = (double*)calloc( N_M, sizeof(double) );
/* Задание возмущения: */
printf
( "Возмущение: a1 * sin ( a2 * t )\n" );
printf ( "a1=" ); scanf ( "%lf", &a1
);( "a2=" ); scanf ( "%lf", &a2 );
/* Открытие файлов: */( "\nFileName: " ); scanf (
"%s", f_Name_u );( f_Name_x, f_Name_u );( f_Name_J, f_Name_u );(
f_Name_u, "_u.dat" );( f_Name_J, "_J.dat" );( f_Name_x,
"_x.dat" );( ( ru = fopen ( f_Name_u, "wt" ) ) == NULL )
{( f_Name_u );(-1);
}( ( rJ = fopen ( f_Name_J, "wt" ) ) == NULL )
{( f_Name_J );(-1);
}( ( rx = fopen ( f_Name_x, "wt" ) ) == NULL )
{( f_Name_x );(-1);
}
tau = TBEG;
/* Печать в файл и на экран x0: */
fprintf ( rx, "%lf ", tau );( "\ntau = %lf
", tau );( "x = ( " );( j=0; j<NX; j++ )
{( rx, "%lf ", x[j] );( "%lf ", x[j] );
}( rx, "\n" ); ( ")\n"
);
/* Начальн. формирование некотор.параметров конечномерн.задачи: */
h = TETA/N; /* шаг(длина перем.)*/
printf ( "h = %lf\n", h );();
/* Работа регулятора: */= 0;( !Ostanov )
{
/* Программн. упpавл.
задачи с перем. структурой: */
puN1NIntU2 ( h );
/* Печ. в файл u(t), x(t) на [tau,tau+h]: */
fprintf ( ru, "%lf %lf\n", tau, u[0] );( ru,
"%lf %lf\n", tau+h, u[0] );( rJ, "%lf %lf\n", tau, J );=
tau;( "tau = %lf\n", tau );( tout < tau+h-0.000001 )
{= tout;( tau + h - tout > 0.1 )+= 0.1;= tau + h;(
"\nt = %lf", t );( " tout = %lf", tout );= 1;= 1;(inte)
{( NX, t, tout, x, f_x, &hx, ypx, f1x, f2x, f3x, f4x,
f5x,
&iflx, &nfex, &kopx, &initx, &jflx,
&kflx,
&rerr, &aerr, &savrex, &savaex);( iflx==2 )=
0;
{("x\n");(&iflx, &rerr, &aerr);
}
}( rx, "%lf ", tout );( "x = ( " );( j=0;
j<NX; j++ )
{( rx, "%lf ", x[j] );( "%lf ", x[j] );
}( rx, "\n" );( ") " );
}+= h;( "\n\n\n tau = %lf\n\n\n", tau );(1)
{( "|x1*x1 + x2*x2 - 1| = %lf\n",
fabs(x[0]*x[0]+x[1]*x[1]-1.0) ); ( "Остановить процесс ? ( 1 -- да, 0 -- нет )" );
scanf ( "%d", &Ostanov );( ( Ostanov==0 ) || (
Ostanov==1 ) );
}
} /* while( !Ostanov ) */();( j=0; j<N_M; j++ )( D[j] );(
D );( j=0; j<N_M; j++ )( Ma[j] );( Ma );0;
}callocPr ( int n )
{i;= (double**)calloc( n, sizeof(double) );( i=0; i<n; i++
)[i] = (double*)calloc( n, sizeof(double) );
}freePr ( int n )
{i;( i=0; i<n; i++ )( Pr[i] );( Pr );
}puN1NIntU2 ( double h )/* Прогр.упр.задачи пер.стр.*/
{
int i;
/* Программн. упр-ие миним. энергии: */
kod4 = 0;( h );( kod4==0 )
{= 0.0;( i=0; i<N; i++ )+= u[i]*u[i]/2.0;*= h;( "J =
%.12lf\n", J );
}
else
{
printf
( "Нет решения !!!\n"
);
getchar();
}
} /* puN1NIntU2 */
void puIntU2 ( double h )/* Прогр.упр.мин.эн. */
{ur,
**H;i, j, k,,,, UslOpt;= (double**)calloc( N_M,
sizeof(double) );( j=0; j<N_M; j++ )[j] = (double*)calloc( N_M,
sizeof(double) );
/* Формирование параметров конечномерн. задачи: */
formParam ( h );
/* Формирование элементов опорного псевдоплана: */
formElem ();
/* Пров-ка согласов-ти оп.пс-плана: */
UslSogl = 1;( j=0; j<knn; j++ )( ( ( Del[ Jnn[j] ] <
-eps1 ) &&
( fabs( u[ Jnn[j] ] - dn[ Jnn[j] ]) < eps2 ) ) ||
( ( Del[ Jnn[j] ] > eps1 ) &&
( fabs( u[ Jnn[j] ] - dw[ Jnn[j] ]) < eps2 ) ) )
{= 0;;
}( "UslSogl = %d\n", UslSogl );( !UslSogl )();= 1;(
j=0; j<M; j++ )( ( u[Jop[j]] < dn[Jop[j]] - eps2 ) ||
( u[Jop[j]] > dw[Jop[j]] + eps2 ) )
{= 0;;
}( !UslOpt )
{
/* Реш.конечномерн.з-чи дв. методом: */
/* Открытие файла для передачи в прогр. *.for (двойств.метод) */( ( d = fopen
( "dat", "wt" ) ) == NULL )
{( "dat" );(-1);
}
/* Выч.обр.м. к м.D */( koc );( j=0; j<koc; j++ )
{( i=0; i<koc; i++ )
{( k=0; k<koc; k++ )[i][k] = D[i][k];[i] = 0.0;
}[j] = 1.0;( koc, Pr, v, &kod );( kod )
{( "\nОпорная матрица (цел.ф-ии) вырождена !!!\n"
);();
}( i=0; i<koc; i++ )[i][j] = v[i];
}( koc );
/* Печать параметров задачи в файл */
fprintf
( d, "%d %d\n", M, N );/* размеры M,N */
for ( j=0; j<N; j++ ) /* матрица A */
{( i=0; i<M; i++ )( d, "%lf ", A[i][j] );( d,
"\n" );
}( i=0; i<M; i++ ) /* вектор b */( d, "%lf ",
b[i] );( d, "\n" );( i=0; i<N; i++ ) /* векторы dn, dw */( d,
"%lf %lf\n", dn[i], dw[i] ); ( d,
"%d\n", M );/*
размер оп.огр. M */
for ( i=0; i<M; i++ ) /* обр.оп.м. Q */
{( j=0; j<M; j++ )( d, "%lf ", Q[i][j] );( d,
"\n" );
}( d, "%d\n", koc ); /* разм.оп.ц.ф. koc */( j=0;
j<koc; j++ )/* верх.треуг.м. H */( i=0; i<=j; i++ )( d,
"%lf\n", H[i][j] );( i=0; i<M; i++ ) /* IX=(Jop,Joc,Jnn) */( d,
"%d\n", Jop[i] + 1 );( i=0; i<koc; i++ )( d, "%d\n",
Joc[i] + 1 );( i=0; i<knn; i++ )( d, "%d\n", Jnn[i] + 1 );(j=0;
j<N; j++ )( d, "%lf\n", u[j] );(d);( "errfile.exe" ); ( "dmqp.exe"
);
/* Открытие файла для чтения из прогр. *.for (двойств.метода) */( ( r =
fopen ( "res", "r" ) ) == NULL )
{( "res" );(-1);
}( r, "%d", &kod4 ); ( "Результаты:\nkod(4) = %d\n", kod4 );
if ( kod4==0 )
{( "Jop = (" );( r, "%d", &i );= i;(
j=0; j<i; j++ )
{( r, "%d", &k );[j] = k - 1;( "%d ",
Jop[j] );
}( " )\n" );( "Joc = (" );( r,
"%d", &i );+= i;( j=0; j<i; j++ )
{( r, "%d", &k );[j] = k - 1;( "%d ",
Joc[j] );
}( " )\n" );( "Jnn = (" );( j=0;
j<N-nn; j++ )
{( r, "%d", &k );[j] = k - 1;( "%d ",
Jnn[j] );
}( " )\n" );( k=0; k<N; k++ )
{( r, "%lf", &ur );[k] = ur;
}( "u =\n" );( i=0; i<10; i++ )(
"%7.1lf", u[i] );( "\n" );( i=10; i<20; i++ )( "%7.1lf",
u[i] );( "\n" );( i=20; i<N; i++ )( "%7.1lf", u[i] );(
"\n" );
} /* if ( kod4==0 ) */(r);
} /* if ( !UslOpt ) */( j=0; j<N_M; j++ )( H[j] );( H );
} /* puIntU2 */formParam ( double h )/*
Форм.парам.конечном.з-чи */
{t_tn, t_tw,, tauNh;i, j;( j=0; j<N; j++ )
{_tn = N*h - j*h;_tw = N*h - (j+1)*h;[0][j] = cos( t_tw ) -
cos( t_tn );[1][j] = sin( t_tn ) - sin( t_tw );
}= N*h;= tau + N*h;[0] = cos( tauNh )*z0[0] + sin( tauNh
)*z0[1];[1] = -sin( tauNh )*z0[0] + cos( tauNh )*z0[1];[0] = z[0] - cos( te
)*x[0] - sin( te )*x[1];[1] = z[1] + sin( te )*x[0] - cos( te )*x[1];( i=0;
i<N; i++ ) /* в-ры dn, dw */
{[i] = -L;[i] = L;
}
} /* formParam */formElem ( void ) /* Форм.эл-тов оп.пс-плана
*/
{i, j, k,,,;
/* Опора огр. Jop */[0] = 0;[1] = N-1;= N - M; /* неопорн.инд.
Jn */= 0;( j=0; j<N; j++ )
{= 0;( i=0; i<M; i++ )( j==Jop[i] )
{= 1;;
}( !PrOp )
{[k] = j;++;
}
}= 0;= kn; /* дв.неоп.инд. Jnn */( j=0; j<knn; j++ )[j] =
Jn[j];( i=0; i<M; i++ ) /* P=A(I,Jop) */( j=0; j<M; j++ )[i][j] =
A[i][Jop[j]];( M );( j=0; j<M; j++ ) /* м. Q -- обр.к P */
{( i=0; i<M; i++ )
{( k=0; k<M; k++ )[i][k] = P[i][k];[i] = 0.0;
}[j] = 1.0;( M, Pr, v, &kod );( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( i=0; i<M; i++ )[i][j] = v[i];
}( M );
/* псевдоплан u */( j=0; j<knn; j++ ) /* unn */[Jnn[j]] =
dn[Jnn[j]];( koc>0 ) /* uoc */
{
/* uoc */
}( koc>0 ) /* U=(*,Uoc,Unn) */( i=0; i<koc; i++
)[Joc[i]] = v[i];( M );( i=0; i<M; i++ ) /* uop */
{( j=0; j<M; j++ )[i][j] = P[i][j];[i] = b[i];( j=0;
j<kn; j++ )[i] -= A[i][Jn[j]] * u[Jn[j]];
}( M, Pr, v, &kod );( M ); ( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( i=0; i<M; i++ ) /* U=(Uop,Uoc,Unn) */[Jop[i]] = v[i];(
M );( i=0; i<M; i++ ) /* Вект.потенц. Up */
{( k=0; k<M; k++ )[i][k] = P[k][i];[i] = u[Jop[i]];
}( M, Pr, v, &kod ); ( M );
if ( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( j=0; j<N; j++ ) /* вект.оц. Delta */
{[j] = u[j];( i=0; i<M; i++ )[j] -= v[i] * A[i][j];
}
} /* formElem */procsogl( void )
{l[N],[M],, Thj;i, j, k,,,,;= 0;(1)
{( "iterS = %d\n", iterS + 1 );
/* l */( j=0; j<knn; j++ ) /* lnn */( Del[Jnn[j]] >
eps1 )[Jnn[j]] = dn[Jnn[j]] - u[Jnn[j]];( Del[Jnn[j]] < -eps1 )[Jnn[j]] =
dw[Jnn[j]] - u[Jnn[j]];[Jnn[j]] = 0.0;( koc>0 ) /* loc */
{( M );( j=0; j<koc; j++ ) /* Выч.матр. Q*Aoc */
{( i=0; i<M; i++ )
{( k=0; k<M; k++ )[i][k] = P[i][k];[i] = A[i][Joc[j]];
}( M, Pr, v, &kod );( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( i=0; i<M; i++ )[i][j] = v[i];
}( M );( i=0; i<koc; i++ ) /* D=E+Aoc'Q'QAoc */( j=0;
j<koc; j++ )
{[i][j] = 0.0;( k=0; k<M; k++ )[i][j] += Ma[k][i] *
Ma[k][j];
}( i=0; i<koc; i++ )[i][i] += 1.0;( M );( i=0; i<M; i++
) /* Qb=Q*Ann*lnn */
{( k=0; k<M; k++ )[i][k] = P[i][k];[i] = 0.0;( j=0;
j<knn; j++ )[i] += A[i][Jnn[j]] * l[Jnn[j]];
}( M, Pr, Qb, &kod ); ( M );
if ( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( i=0; i<koc; i++ ) /* boc=Aoc'Q'Qb */
{[i] = 0.0;( j=0; j<M; j++ )[i] -= Ma[j][i] * Qb[j];
}( koc );( i=0; i<koc; i++ ) /* Выч. loc */( k=0;
k<koc; k++ )[i][k] = D[i][k];( koc, Pr, v, &kod );( koc );( kod )
{( "\nОпорная матрица (цел.ф-ии) вырождена !!!\n"
);();
}
}( j=0; j<koc; j++ )/* l=(*,loc,lnn) */[Joc[j]] = v[j];( M
);( i=0; i<M; i++ )/* lop */
{( j=0; j<M; j++ )[i][j] = P[i][j];[i] = 0.0;( j=0;
j<kn; j++ )[i] -= A[i][Jn[j]] * l[Jn[j]];
}( M, Pr, v, &kod );( M ); ( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( j=0; j<M; j++ ) /* l=(lop,loc,lnn) */[Jop[j]] = v[j];
/* tnn */( M );( i=0; i<M; i++ ) /* Qb=Q'*lop */
{( k=0; k<M; k++ )[i][k] = P[k][i];[i] = l[Jop[i]];
}( M, Pr, Qb, &kod ); ( M );
if ( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( j=0; j<knn; j++ ) /* tnn */
{[j] = l[Jnn[j]];( i=0; i<M; i++ )[j] -= Qb[i] *
A[i][Jnn[j]];
}= 1; /* Theta */= -1;( j=0; j<koc; j++ ) /* Theta(Joc)
*/( l[Joc[j]] < -eps2 )
{= ( dn[Joc[j]] - u[Joc[j]] ) / l[Joc[j]];( Thj < Th )
{= Thj;= j; /* ?????? */= 0;
}
}( l[Joc[j]] > eps2 )
{= ( dw[Joc[j]] - u[Joc[j]] ) / l[Joc[j]];( Thj < Th )
{= Thj;= j;= 0;
}
}( j=0; j<knn; j++ ) /* Th(Jnn) */( Del[Jnn[j]] * v[j]
< -eps1 )
{= -Del[Jnn[j]] / v[j];( Thj < Th )
{= Thj;= j;= 1;
}
}( j=0; j<N; j++ ) /* u=u+Th*l */[j] += Th*l[j];( Th == 1
)
{;
}( prTh )
{[koc] = Jnn[j0];++;-;( i=j0; i<knn; i++ )[i] = Jnn[i+1];
}
{[knn] = Joc[j0];++;-;( i=j0; i<koc; i++ )[i] = Joc[i+1];
}( M );( i=0; i<M; i++ ) /* Вект.потенц. Up */
{( k=0; k<M; k++ )[i][k] = P[k][i];[i] = u[Jop[i]];
}( M, Pr, v, &kod ); ( M );
if ( kod )
{
printf
( "\nОпорная матрица (ограничений)
вырождена !!!\n" );
getchar();
}( j=0; j<N; j++ ) /* вект.оц. Delta */
{[j] = u[j];( i=0; i<M; i++ )[j] -= v[i] * A[i][j];
}++;
} /* while(1) */
}gs ( int m, double **pX,copX[], int *kod )
{tol = 1E-5,, pXpX,;i, j, k, ky, imax,, jx, i0, ib;( j=0;
j<m; j++ )
{= 0.;( i=j; i<m; i++ )
{= fabs ( pX [i] [j] );( fabs ( bigpX ) < pXpX )
{= pX [i] [j];= i;
}
}( fabs ( bigpX ) <= tol )
{
*kod = 1;1;
}( k=j; k<m; k++ )
{= pX [imax] [k];[imax] [k] = pX [j] [k];[j] [k] = sav /
bigpX;
}= copX [imax];[imax] = copX [j];[j] = sav / bigpX;( j==m-1
);( ix=j+1; ix < m; ix++ )
{( jx=j+1; jx < m; jx++ )[ix] [jx] -= pX [ix] [j] * pX [j]
[jx];[ix] -= copX [j] * pX [ix] [j];
}
}( ky=0; ky < m-1; ky++ )
{= m-2-ky;= m-1;( k=0; k <= ky; k++ )
{[ib] -= pX [ib] [i0] * copX [i0];--;
}
}
*kod=0;0;
}f_x ( double t, /* Выч. f1(t) = x2(t) */y[], /* f2(t) =
-x1(t) + u(t) */f[] )
{[0] = y[1];[1] = -y[0] + u[0] + a1*sin( a2*t );
} /* f_x() */
# define True 1
# define False 0FEHL(int neqn, double t,y[], void
(f_p)(double, double*, double*),*hc, double yp[],f1[], double f2[],f3[], double
f4[],f5[], double s[])
{ double h;ch,rab;k;(*fun)(double, double*, double*);=*hc;=
f_p;= h/4.0;(k=0; k<neqn; k++)[k] = y[k] + ch * yp[k];= t + ch;(rab,f5,f1);=
3.0 * h / 32.0;(k=0; k<neqn; k++)[k] = y[k] + ch * (yp[k] + 3.0 * f1[k]);= t
+ 3.0 * h / 8.0;(rab,f5,f2);= h / 2197.0;(k=0; k<neqn; k++)[k] = y[k] + ch *
(1932.0 * yp[k] + (7296.0 * f2[k] - 7200.0 * f1[k]));= t + 12.0 * h /
13.0;(rab,f5,f3);= h / 4104.0;(k=0; k<neqn; k++)[k] = y[k] + ch * ((8341.0 *
yp[k] - 845.0 * f3[k]) +
(29440.0 * f2[k] - 32832.0 * f1[k]));= t + h;(rab,f5,f4);= h
/ 20520.0;(k=0; k<neqn; k++)[k] = y[k] + ch * ((-6080.0 * yp[k] + (9295.0 *
f3[k] -
.0 * f4[k])) + (41040.0 * f1[k] - 28352.0 * f2[k]));= t + h /
2.0;(rab,f1,f5);
/*ВЫЧИСЛИТЬ ПРИБЛИЖЕННОЕ РЕШЕНИЕ В ТОЧКЕ t + h*/
ch = h / 7618050.0;(k=0; k<neqn; k++)[k] = y[k] + ch *
((902880.0 * yp[k] + (3855735.0 * f3[k] -
.0 * f4[k])) + (3953664.0 * f2[k] + 277020.0 * f5[k]));
*hc=h;;
}RKF45 (int neqn, double t,tout, double y[],(f_p)(double,
double*, double*), double *hc,yp[], double f1[],f2[], double f3[],f4[], double
f5[],*iflagc, int *nfec,*kopc, int *initc,*jflagc, int *kflagc,*relerrc, double
*abserrc,*savrec, double *savaec)
{iflag,nfe,kop,init,jflag,kflag;h,savre,savae,, abserr;maxnfe
= 30000;remin = 1.0e-15;hfaild,
output;a,ae,dt,ee,eeoet,esttol,et,hmin,eps,u26,,s,scale,tol,toln,ypk,rab1;k,
mflag;(*fun)(double, double*, double*);= f_p;=*iflagc; nfe=*nfec; kop=*kopc;
init=*initc; jflag=*jflagc;=*kflagc;=*hc; savre=*savrec; savae=*savaec;=
*relerrc;= *abserrc;
/*Проверить входные параметры*/(neqn < 1)m10;(relerr <
0.0 || abserr < 0.0)m10;=abs(iflag);(mflag == 0 || mflag > 8)m10; (mflag != 1)
goto m20;
/*Первый вызов, вычислить машинное эпсилон*/
eps = 1.0;(eps+1.0 > 1.0)= eps/2.0;26 = 26.0*eps;
goto m50;
/*Ошибки во входной информации*/
m10: iflag = 8;
goto mexit;
/*Проверить возможность продолжения*/
m20: if (t==tout && kflag != 3)m10;(mflag !=
2)m25;(kflag == 3 || init==0)m45;(kflag == 4)m40;(kflag == 5 && abserr
== 0.0)m30;(kflag == 6 && relerr <= savre && abserr <=
savae)m30;m50;: if (iflag == 3)m45;(iflag == 4)m40;(iflag == 5 &&
abserr > 0.0)m45;: goto mexit;: nfe = 0;(mflag == 2)m50;: iflag =
jflag;(kflag == 3)= abs(iflag);: jflag = iflag;= 0;= relerr;= abserr;=
2.0*eps+remin;(relerr >= rer)m55;
/*Заданная граница относительной погрешности слишком мала*/
relerr = rer;= 3;= 3;mexit;: dt =
tout-t;(mflag==1)m60;(init==0)m65;m80;: init = 0;= 0;= t;(a,y,yp);= 1;(t !=
tout)m65;= 2;mexit;: init = 1;= fabs(dt);= 0.0;(k=0; k<neqn; k++)
{ tol = relerr*fabs(y[k])+abserr;(tol > 0.0)
{ toln=tol;=fabs(yp[k]);=h*h*h*h*h;(ypk*rab1 >
tol)=exp(0.2*log(tol/ypk));
}
}(toln <= 0.0)=0.0;=u26*fabs(t);(fabs(t) <
fabs(dt))=u26*fabs(dt);(h < rab1)=rab1;=-2;(iflag > 0)=2;: if (dt <
0.0)=-h;(fabs(h) >= 2.0*fabs(dt))=kop+1;(kop != 100 )m85;
/*Много выходов*/=0;=7;mexit;: if (fabs(dt) > u26*fabs(t)) m95;
/*Близко к выходной точке, проэкстраполировать и вернуться по месту
вызова*/
for (k=0; k<neqn; k++)[k] = y[k] + dt *
yp[k];=tout;(a,y,yp);=nfe+1;m300;: output=False;=2.0/relerr;=scale*abserr;
/*Пошаговое интегрирование*/
m100: hfaild=False;
hmin=u26*fabs(t);=tout-t;(fabs(dt) >= 2.0*fabs(h))m200;(fabs(dt)
> fabs(h))m150;=True;=dt; m200;
m150: h=0.5*dt;
/*Внутренний одношаговый интегратор*/
m200: if(nfe <= maxnfe)m220;=4;
kflag=4;
goto mexit;
/*Продолжить приближенное решение на один шаг длины h*/
m220: FEHL(neqn, t, y, f_p, &h, yp, f1, f2, f3, f4, f5, f1);
nfe=nfe+5;=0.0;(k=0; k<neqn; k++)
{ et=fabs(y[k])+fabs(f1[k])+ae; (et
> 0.0)
goto m240;
/*Неправильная граница погрешности*/
iflag=5;mexit;:
ee=fabs((-2090.0*yp[k]+(21970.0*f3[k]-15048.0*f4[k]))+
(22528.0*f2[k]-27360.0*f5[k]));=ee/et;(eeoet < rab1)=rab1;
}=fabs(h)*eeoet*scale/752400.0; (esttol
<= 1.0)
goto m260;
/*Неудачный шаг; уменьшить его величину*/
hfaild=True;=False;=0.1;(esttol <
59049.0)=0.9/(exp(0.2*log(esttol)));=s*h;(fabs(h) > hmin) m200;
/*Неудача даже при минимальном шаге*/
iflag=6;
kflag=6;
goto mexit;
/*Успешный шаг*/
m260: t=t+h;(k=0; k<neqn;
k++)[k]=f1[k];=t;(a,y,yp);=nfe+1;=5.0;(esttol >
1.889568e-4)=0.9/(exp(0.2*log(esttol)));(hfaild)
{ if (s > 1.0)=1.0;
}=s*fabs(h);(rab1 < hmin)=hmin;(h >
0.0)=fabs(rab1);=-fabs(rab1);(output)m300;(iflag > 0)m100;
/*Интегрирование успешно завершено*/
iflag=-2;
m300: iflag=2;
mexit:
*iflagc=iflag; *nfec=nfe; *kopc=kop; *initc=init; *jflagc=jflag;
*kflagc=kflag;
*hc=h; *savrec=savre; *savaec=savae;
*relerrc = relerr;
*abserrc = abserr;;
}prifl0 ( int *iflc,*rerrc,*aerrc )
{;,;= *iflc;= *rerrc;= *aerrc;(ifl)
{3: break;4: printf ("много fp()\n ifl=%d\n",
ifl);;5: aerr = 1.0e-9;;6: rerr *= 10.0;("ifl=%d\n rerr=%.12lf
aerr=%.12lf\n", ifl, rerr, aerr);= 2;;7: printf ("много вых. точек\n
ifl=%d\n", ifl);= 2;;8: printf ("ifl=%d\n", ifl);;
}
*iflc = ifl;
*rerrc = rerr;
*aerrc = aerr;
}
Приложение Б
Текст программы 2
#pragma hdrstop
# include <stdio.h>
# include <math.h>
# include <string.h>
# include <process.h>
# include <alloc.h>
# include <conio.h>
#define M 2
#define N 25
#define N_M 23
#define TBEG0.0
#define TETA 5.0
#define L 3.0
#define Uf -2.0
#define NX 2tau,[NX] = { 2.0, 4.0 },[NX],[NX] = { 2.0, 3.0
},[M], A[M][N], dn[N], dw[N],= 0.00001, eps2=0.00001,[M][M], **Pr, Q[M][M],
**D,v[N_M], **Ma,u[N],[N], J, a1, a2;Jop[M],[N_M], kn,[N_M],
koc,[N_M], knn,;callocPr ( int );freePr ( int );puN1NIntU2 ( double );puIntU2 (
double );formParam ( double );formElem ( void );procsogl( void );gs ( int,
double**, /* Метод Гаусса */*, int* );f_x ( double, double*,/* f1(t) = x2(t) + u(t)
*/* ); /* f2(t) = -x1(t) - u(t) */FEHL(int, double,*, void (f_p)(double,
double*, double*),*, double*,*, double*,*, double*,*, double*);RKF45 (int,
double, /* Метод Рунге-Кутта-Фельдберга */, double*,(f_p)(double, double*,
double*), double*,*, double*,*, double*,*, double*,*, int*,*, int*,*, int*,*,
double*,*, double*);prifl0 (int*, double*, double*);*d;*r;*ru; /* Файл.dat : t, u(t) */*rJ; /* Файл.dat : t, J(t) */*rx; /* Файл.dat : t, x(t) */callocPr ( int n )
{i;= (double**)calloc( n, sizeof(double) );( i=0; i<n; i++
)[i] = (double*)calloc( n, sizeof(double) );
}freePr ( int n )
{i;( i=0; i<n; i++ )( Pr[i] );( Pr );
}puN1NIntU2 ( double h )/* Прогр.упр.задачи пер.стр.*/
{i;
/* Программн. упр-ие миним. энергии: */
kod4 = 0;( h );( kod4==0 )
{= 0.0;( i=0; i<N; i++ )+= u[i]*u[i]/2.0;*= h;( "J =
%.12lf\n", J );
}
else
{
printf
( "Нет решения !!!\n"
);
getchar();
}
} /* puN1NIntU2 */
void puIntU2 ( double h )/* Прогр.упр.мин.эн. */
{ur, **H;i, j, k, kod, nn, UslSogl, UslOpt;=
(double**)calloc( N_M, sizeof(double) );( j=0; j<N_M; j++ )[j] =
(double*)calloc( N_M, sizeof(double) );
/* Формирование параметров конечномерн. задачи: */( h );
/* Формирование элементов опорного псевдоплана: */();
/* Пров-ка согласов-ти оп.пс-плана: */
UslSogl = 1;( j=0; j<knn; j++ )( ( ( Del[ Jnn[j] ] <
-eps1 ) &&
( fabs( u[ Jnn[j] ] - dn[ Jnn[j] ]) < eps2 ) ) ||
( ( Del[ Jnn[j] ] > eps1 ) &&
( fabs( u[ Jnn[j] ] - dw[ Jnn[j] ]) < eps2 ) ) )
{= 0;;
}( "UslSogl = %d\n", UslSogl );( !UslSogl )();
/*printf ( "A * u - b = (" );= 1;( j=0; j<M; j++
)( ( u[Jop[j]] < dn[Jop[j]] - eps2 ) ||
( u[Jop[j]] > dw[Jop[j]] + eps2 ) )
{= 0;;
}( !UslOpt )
{
/* Реш.конечномерн.з-чи дв. методом: */
/* Открытие файла для передачи в прогр. *.for (двойств.метод) */( ( d = fopen (
"dat", "wt" ) ) == NULL )
{( "dat" );(-1);
}
/* Выч.обр.м. к м.D */( koc );( j=0; j<koc; j++ )
{( i=0; i<koc; i++ )
{( k=0; k<koc; k++ )[i][k] = D[i][k];[i] = 0.0;
}[j] = 1.0;( koc, Pr, v, &kod );( kod )
{( "\nОпорная матрица (цел.ф-ии) вырождена !!!\n" );();
}( i=0; i<koc; i++ )[i][j] = v[i];
}( koc );
/* Печать параметров задачи в файл */( d, "%d %d\n", M, N );/*
размеры M,N */
for ( j=0; j<N; j++ ) /* матрица A */
{( i=0; i<M; i++ )( d, "%lf ", A[i][j] );( d,
"\n" );
}( i=0; i<M; i++ ) /* вектор b */( d, "%lf ", b[i] );( d, "\n" );(
i=0; i<N; i++ ) /* векторы dn, dw */( d, "%lf
%lf\n", dn[i], dw[i] );
fprintf ( d, "%d\n", M );/* размер оп.огр. M */
for ( i=0; i<M; i++ ) /* обр.оп.м. Q */
{( j=0; j<M; j++ )( d, "%lf ", Q[i][j] );( d,
"\n" );
}( d, "%d\n", koc ); /* разм.оп.ц.ф. koc */( j=0; j<koc; j++ )/* верх.треуг.м. H */( i=0; i<=j; i++
)( d, "%lf\n", H[i][j] );( i=0; i<M; i++ ) /* IX=(Jop,Joc,Jnn) */(
d, "%d\n", Jop[i] + 1 );( i=0; i<koc; i++ )( d, "%d\n",
Joc[i] + 1 );( i=0; i<knn; i++ )( d, "%d\n", Jnn[i] + 1 );(j=0;
j<N; j++ )( d, "%lf\n", u[j] );(d);( "errfile.exe" ); ( "dmqp.exe"
);
/* Открытие файла для чтения из прогр. *.for (двойств.метода) */( ( r = fopen (
"res", "r" ) ) == NULL )
{( "res" );(-1);
}( r, "%d", &kod4 ); ( "Результаты:\nkod(4) = %d\n", kod4 );
if ( kod4==0 )
{( "Jop = (" );( r, "%d", &i );= i;(
j=0; j<i; j++ )
{( r, "%d", &k );[j] = k - 1;( "%d ",
Jop[j] );
}( " )\n" );( "Joc = (" );( r,
"%d", &i );+= i;( j=0; j<i; j++ )
{( r, "%d", &k );[j] = k - 1;( "%d ",
Joc[j] );
}( " )\n" );( "Jnn = (" );( j=0;
j<N-nn; j++ )
{( r, "%d", &k );[j] = k - 1;( "%d ",
Jnn[j] );
}( " )\n" );( k=0; k<N; k++ )
{( r, "%lf", &ur );[k] = ur;
}( "u =\n" );( i=0; i<10; i++ )(
"%7.1lf", u[i] );( "\n" );( i=10; i<20; i++ )(
"%7.1lf", u[i] );( "\n" );( i=20; i<N; i++ )(
"%7.1lf", u[i] );( "\n" );
} /* if ( kod4==0 ) */(r);
} /* if ( !UslOpt ) */( j=0; j<N_M; j++ )( H[j] );( H );
} /* puIntU2 */formParam ( double h )/* Форм.парам.конечном.з-чи */
{t_tn, t_tw,;i, j;( j=0; j<N; j++ )
{_tn = N*h - j*h;_tw = N*h - (j+1)*h;[0][j] = sin( t_tn ) -
sin( t_tw ) + cos( t_tn ) - cos( t_tw );[1][j] = cos( t_tn ) - cos( t_tw ) +
sin( t_tw ) - sin( t_tn );
}= N*h;[0] = x[0] - ( z0[1] + Uf )*sin( tau ) - ( z0[0] + Uf
)*cos( tau ) + Uf;[1] = x[1] + ( z0[0] + Uf )*sin( tau ) - ( z0[1] + Uf )*cos(
tau ) + Uf;[0] = -cos( te )*yx[0] - sin( te )*yx[1];[1] = sin( te )*yx[0] -
cos( te )*yx[1];( i=0; i<N; i++ ) /* в-ры dn, dw */
{[i] = -L-Uf;[i] = L-Uf;
}
} /* formParam */formElem ( void ) /* Форм.эл-тов оп.пс-плана */
{i, j, k, jr, kod, PrOp;
/* Опора огр. Jop */[0] = 0;[1] = 24;= N - M; /* неопорн.инд. Jn */= 0;( j=0; j<N; j++ )
{= 0;( i=0; i<M; i++ )( j==Jop[i] )
{= 1;;
}( !PrOp )
{[k] = j;++;
}
}= 0;= kn; /* дв.неоп.инд. Jnn */( j=0; j<knn; j++ )[j] =
Jn[j];( i=0; i<M; i++ ) /* P=A(I,Jop) */( j=0; j<M; j++ )[i][j] =
A[i][Jop[j]];( M );( j=0; j<M; j++ ) /* м. Q -- обр.к P */
{( i=0; i<M; i++ )
{( k=0; k<M; k++ )[i][k] = P[i][k];[i] = 0.0;
}[j] = 1.0;( M, Pr, v, &kod );( kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( i=0; i<M; i++ )[i][j] = v[i];
}( M );
/* псевдоплан u */( j=0; j<knn; j++
) /* unn */[Jnn[j]] = dn[Jnn[j]];( koc>0 ) /* uoc */
{
/* uoc */
}( koc>0 ) /* U=(*,Uoc,Unn) */( i=0; i<koc; i++
)[Joc[i]] = v[i];( M );( i=0; i<M; i++ ) /* uop */
{( j=0; j<M; j++ )[i][j] = P[i][j];[i] = b[i];( j=0;
j<kn; j++ )[i] -= A[i][Jn[j]] * u[Jn[j]];
}( M, Pr, v, &kod );( M );
if ( kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( i=0; i<M; i++ ) /* U=(Uop,Uoc,Unn) */[Jop[i]] = v[i];(
M );( i=0; i<M; i++ ) /* Вект.потенц. Up */
{( k=0; k<M; k++ )[i][k] = P[k][i];[i] = u[Jop[i]];
}( M, Pr, v, &kod ); ( M );(
kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( j=0; j<N; j++ ) /* вект.оц. Delta */
{[j] = u[j];( i=0; i<M; i++ )[j] -= v[i] * A[i][j];
}
} /* formElem */procsogl( void )
{l[N], Qb[M], Th, Thj;i, j, k, kod, iterS, prTh, j0;= 0;(1)
{( "iterS = %d\n", iterS + 1 );
/* l */( j=0; j<knn; j++ ) /* lnn */( Del[Jnn[j]] >
eps1 )[Jnn[j]] = dn[Jnn[j]] - u[Jnn[j]];( Del[Jnn[j]] < -eps1 )[Jnn[j]] =
dw[Jnn[j]] - u[Jnn[j]];[Jnn[j]] = 0.0;( koc>0 ) /* loc */
{( M );( j=0; j<koc; j++ ) /* Выч.матр. Q*Aoc */
{( i=0; i<M; i++ )
{( k=0; k<M; k++ )[i][k] = P[i][k];[i] = A[i][Joc[j]];
}( M, Pr, v, &kod );( kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( i=0; i<M; i++ )[i][j] = v[i];
}( M );( i=0; i<koc; i++ ) /* D=E+Aoc'Q'QAoc */( j=0;
j<koc; j++ )
{[i][j] = 0.0;( k=0; k<M; k++ )[i][j] += Ma[k][i] *
Ma[k][j];
}( i=0; i<koc; i++ )[i][i] += 1.0;( M );( i=0; i<M; i++
) /* Qb=Q*Ann*lnn */
{( k=0; k<M; k++ )[i][k] = P[i][k];[i] = 0.0;( j=0;
j<knn; j++ )[i] += A[i][Jnn[j]] * l[Jnn[j]];
}( M, Pr, Qb, &kod );
freePr ( M );( kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( i=0; i<koc; i++ ) /* boc=Aoc'Q'Qb */
{[i] = 0.0;( j=0; j<M; j++ )[i] -= Ma[j][i] * Qb[j];
}( koc );( i=0; i<koc; i++ ) /* Выч. loc */( k=0; k<koc; k++ )[i][k]
= D[i][k];( koc, Pr, v, &kod );( koc );( kod )
{( "\nОпорная матрица (цел.ф-ии) вырождена !!!\n" );();
}
}( j=0; j<koc; j++ )/* l=(*,loc,lnn) */[Joc[j]] = v[j];( M
);( i=0; i<M; i++ )/* lop */
{( j=0; j<M; j++ )[i][j] = P[i][j];[i] = 0.0;( j=0;
j<kn; j++ )[i] -= A[i][Jn[j]] * l[Jn[j]];
}( M, Pr, v, &kod );( M );
if ( kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( j=0; j<M; j++ ) /* l=(lop,loc,lnn) */[Jop[j]] = v[j];
/* tnn */( M );( i=0; i<M; i++ ) /* Qb=Q'*lop */
{( k=0; k<M; k++ )[i][k] = P[k][i];[i] = l[Jop[i]];
}( M, Pr, Qb, &kod );
freePr ( M );( kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( j=0; j<knn; j++ ) /* tnn */
{[j] = l[Jnn[j]];( i=0; i<M; i++ )[j] -= Qb[i] *
A[i][Jnn[j]];
}= 1; /* Theta */= -1;( j=0; j<koc; j++ ) /* Theta(Joc)
*/( l[Joc[j]] < -eps2 )
{= ( dn[Joc[j]] - u[Joc[j]] ) / l[Joc[j]];( Thj < Th )
{= Thj;= j; /* ?????? */= 0;
}
}( l[Joc[j]] > eps2 )
{= ( dw[Joc[j]] - u[Joc[j]] ) / l[Joc[j]];( Thj < Th )
{= Thj;= j;= 0;
}
}( j=0; j<knn; j++ ) /* Th(Jnn) */( Del[Jnn[j]] * v[j]
< -eps1 )
{= -Del[Jnn[j]] / v[j];( Thj < Th )
{= Thj;= j;= 1;
}
}( j=0; j<N; j++ ) /* u=u+Th*l */[j] += Th*l[j];( Th == 1
)
{;
}( prTh )
{[koc] = Jnn[j0];++;-;( i=j0; i<knn; i++ )[i] = Jnn[i+1];
}
{[knn] = Joc[j0];++;-;( i=j0; i<koc; i++ )[i] = Joc[i+1];
}( M );( i=0; i<M; i++ ) /* Вект.потенц. Up */
{( k=0; k<M; k++ )[i][k] = P[k][i];[i] = u[Jop[i]];
}( M, Pr, v, &kod ); ( M );(
kod )
{( "\nОпорная матрица (ограничений) вырождена !!!\n" );
getchar();
}( j=0; j<N; j++ ) /* вект.оц. Delta */
{[j] = u[j];( i=0; i<M; i++ )[j] -= v[i] * A[i][j];
}++;
} /* while(1) */
}gs ( int m, double **pX,copX[], int *kod )
{tol = 1E-5,, pXpX,;i, j, k, ky, imax,, jx, i0, ib;( j=0;
j<m; j++ )
{= 0.;( i=j; i<m; i++ )
{= fabs ( pX [i] [j] );( fabs ( bigpX ) < pXpX )
{= pX [i] [j];= i;
}
}( fabs ( bigpX ) <= tol )
{
*kod = 1;1;
}( k=j; k<m; k++ )
{= pX [imax] [k];[imax] [k] = pX [j] [k];[j] [k] = sav /
bigpX;
}= copX [imax];[imax] = copX [j];[j] = sav / bigpX;( j==m-1
);( ix=j+1; ix < m; ix++ )
{( jx=j+1; jx < m; jx++ )[ix] [jx] -= pX [ix] [j] * pX [j]
[jx];[ix] -= copX [j] * pX [ix] [j];
}
}( ky=0; ky < m-1; ky++ )
{= m-2-ky;= m-1;( k=0; k <= ky; k++ )
{[ib] -= pX [ib] [i0] * copX [i0];--;
}
}
*kod=0;0;
}f_x ( double t, /* Выч. f1(t) = x2(t) + u(t) */y[], /* f2(t) = -x1(t) - u(t) */f[]
)
{[0] = y[1] + ( u[0] + Uf );[1] = -y[0] - ( u[0] + Uf ) +
a1*sin( a2*t );
} /* f_x() */
# define True 1
# define False 0FEHL(int neqn, double t,y[], void
(f_p)(double, double*, double*),*hc, double yp[],f1[], double f2[],f3[], double
f4[],f5[], double s[])
{ double h;ch,rab;k;(*fun)(double, double*, double*);=*hc;=
f_p;= h/4.0;(k=0; k<neqn; k++)[k] = y[k] + ch * yp[k];= t + ch;(rab,f5,f1);=
3.0 * h / 32.0;(k=0; k<neqn; k++)[k] = y[k] + ch * (yp[k] + 3.0 * f1[k]);= t
+ 3.0 * h / 8.0;(rab,f5,f2);= h / 2197.0;(k=0; k<neqn; k++)[k] = y[k] + ch *
(1932.0 * yp[k] + (7296.0 * f2[k] - 7200.0 * f1[k]));= t + 12.0 * h /
13.0;(rab,f5,f3);= h / 4104.0;(k=0; k<neqn; k++)[k] = y[k] + ch * ((8341.0 *
yp[k] - 845.0 * f3[k]) +
(29440.0 * f2[k] - 32832.0 * f1[k]));= t + h;(rab,f5,f4);= h
/ 20520.0;(k=0; k<neqn; k++)[k] = y[k] + ch * ((-6080.0 * yp[k] + (9295.0 *
f3[k] -
.0 * f4[k])) + (41040.0 * f1[k] - 28352.0 * f2[k]));= t + h /
2.0;(rab,f1,f5);
/*ВЫЧИСЛИТЬ ПРИБЛИЖЕННОЕ РЕШЕНИЕ В ТОЧКЕ t + h*/
ch = h / 7618050.0;(k=0; k<neqn; k++)[k] = y[k] + ch *
((902880.0 * yp[k] + (3855735.0 * f3[k] -
.0 * f4[k])) + (3953664.0 * f2[k] + 277020.0 * f5[k]));
*hc=h;;
}RKF45 (int neqn, double t,tout, double y[],(f_p)(double,
double*, double*), double *hc,yp[], double f1[],f2[], double f3[],f4[], double
f5[],*iflagc, int *nfec,*kopc, int *initc,*jflagc, int *kflagc,*relerrc, double
*abserrc,*savrec, double *savaec)
{iflag,nfe,kop,init,jflag,kflag;h,savre,savae,, abserr;maxnfe
= 30000;remin = 1.0e-15;hfaild,
output;a,ae,dt,ee,eeoet,esttol,et,hmin,eps,u26,,s,scale,tol,toln,ypk,rab1;k,
mflag;(*fun)(double, double*, double*);= f_p;=*iflagc; nfe=*nfec; kop=*kopc;
init=*initc; jflag=*jflagc;=*kflagc;=*hc; savre=*savrec; savae=*savaec;=
*relerrc;= *abserrc;
/*Проверить входные
параметры*/(neqn < 1)m10;(relerr
< 0.0 || abserr < 0.0)m10;=abs(iflag);(mflag == 0 || mflag > 8)m10; (mflag != 1)m20;
/*Первый вызов, вычислить машинное эпсилон*/
eps = 1.0;(eps+1.0 > 1.0)= eps/2.0;26 = 26.0*eps;
goto m50;
/*Ошибки во входной информации*/: iflag = 8;mexit;
/*Проверить возможность продолжения*/
m20: if (t==tout && kflag != 3)m10;(mflag !=
2)m25;(kflag == 3 || init==0)m45;(kflag == 4)m40;(kflag == 5 && abserr
== 0.0)m30;(kflag == 6 && relerr <= savre && abserr <=
savae)m30;m50;: if (iflag == 3)m45;(iflag == 4)m40;(iflag == 5 &&
abserr > 0.0)m45;: goto mexit;: nfe = 0;(mflag == 2)m50;: iflag =
jflag;(kflag == 3)= abs(iflag);: jflag = iflag;= 0;= relerr;= abserr;=
2.0*eps+remin;(relerr >= rer)m55;
/*Заданная граница относительной погрешности слишком мала*/
relerr = rer;= 3;= 3;mexit;: dt =
tout-t;(mflag==1)m60;(init==0)m65;m80;: init = 0;= 0;= t;(a,y,yp);= 1;(t !=
tout)m65;= 2;mexit;: init = 1;= fabs(dt);= 0.0;(k=0; k<neqn; k++)
{ tol = relerr*fabs(y[k])+abserr;(tol > 0.0)
{ toln=tol;=fabs(yp[k]);=h*h*h*h*h;(ypk*rab1 >
tol)=exp(0.2*log(tol/ypk));
}
}(toln <= 0.0)=0.0;=u26*fabs(t);(fabs(t) <
fabs(dt))=u26*fabs(dt);(h < rab1)=rab1;=-2;(iflag > 0)=2;: if (dt <
0.0)=-h;(fabs(h) >= 2.0*fabs(dt))=kop+1;(kop != 100 )m85;
/*Много выходов*/=0;=7;mexit;: if (fabs(dt) > u26*fabs(t))
goto m95;
/*Близко к выходной точке, проэкстраполировать и вернуться по месту
вызова*/
for (k=0; k<neqn; k++)[k] = y[k] + dt *
yp[k];=tout;(a,y,yp);=nfe+1;m300;: output=False;=2.0/relerr;=scale*abserr;
/*Пошаговое интегрирование*/: hfaild=False;
hmin=u26*fabs(t);=tout-t;(fabs(dt) >=
2.0*fabs(h))m200;(fabs(dt) > fabs(h))m150;=True;=dt; m200;
m150: h=0.5*dt;
/*Внутренний одношаговый интегратор*/
m200: if(nfe <= maxnfe)m220;
iflag=4;=4;mexit;
/*Продолжить приближенное решение на один шаг длины h*/: FEHL(neqn, t, y,
f_p, &h, yp, f1, f2, f3, f4, f5, f1);
nfe=nfe+5;=0.0;(k=0; k<neqn; k++)
{ et=fabs(y[k])+fabs(f1[k])+ae;
if (et > 0.0)m240;
/*Неправильная граница погрешности*/
iflag=5;mexit;:
ee=fabs((-2090.0*yp[k]+(21970.0*f3[k]-15048.0*f4[k]))+
(22528.0*f2[k]-27360.0*f5[k]));=ee/et;(eeoet < rab1)=rab1;
}=fabs(h)*eeoet*scale/752400.0;
if (esttol <= 1.0)m260;
/*Неудачный шаг; уменьшить его величину*/
hfaild=True;=False;=0.1;(esttol < 59049.0)=0.9/(exp(0.2*log(esttol)));=s*h;(fabs(h)
> hmin) m200;
/*Неудача даже при минимальном шаге*/=6;=6;mexit;
/*Успешный шаг*/
m260: t=t+h;(k=0; k<neqn;
k++)[k]=f1[k];=t;(a,y,yp);=nfe+1;=5.0;(esttol >
1.889568e-4)=0.9/(exp(0.2*log(esttol)));(hfaild)
{ if (s > 1.0)=1.0;
}=s*fabs(h);(rab1 < hmin)=hmin;(h >
0.0)=fabs(rab1);=-fabs(rab1);(output)m300;(iflag > 0)m100;
/*Интегрирование успешно завершено*/
iflag=-2;
goto mexit;
m300: iflag=2;
mexit:
*iflagc=iflag; *nfec=nfe; *kopc=kop; *initc=init; *jflagc=jflag;
*kflagc=kflag;
*hc=h; *savrec=savre; *savaec=savae;
*relerrc = relerr;
*abserrc = abserr;;
}prifl0 ( int *iflc,*rerrc,*aerrc )
{ifl;rerr, aerr;= *iflc;= *rerrc;= *aerrc;(ifl)
{3: break;4: printf ("много fp()\n ifl=%d\n", ifl);;5: aerr = 1.0e-9;;6: rerr *=
10.0;("ifl=%d\n rerr=%.12lf aerr=%.12lf\n", ifl, rerr, aerr);= 2;;7:
printf ("много вых. точек\n ifl=%d\n", ifl);=
2;;8: printf ("ifl=%d\n", ifl);;
}
*iflc = ifl;
*rerrc = rerr;
*aerrc = aerr;
}
#pragma argsusedmain(int argc, char* argv[])
{ char f_Name_u[13],/* Имя файла.dat : t, u(t) */_Name_J[13],/* Имя файла.dat
: t, J(t) */_Name_x[13];/* Имя файла.dat : t, x(t) */h, t, tout, ypx
[NX], hx=0, f1x [NX], f2x [NX], f3x [NX], f4x [NX], f5x [NX], rerr = 1.0e-12,
aerr = 0.0,= 0, savaex = 0;j, Ostanov, inte, iflx, jflx, kflx, nfex, kopx,
initx;= (double**)calloc( N_M, sizeof(double) );( j=0; j<N_M; j++ )[j] =
(double*)calloc( N_M, sizeof(double) );= (double**)calloc( N_M, sizeof(double)
);( j=0; j<N_M; j++ )[j] = (double*)calloc( N_M, sizeof(double) );
/* Задание возмущения: */
printf
( "Возмущение: a1 * sin ( a2 * t )\n" );
printf ( "a1=" ); scanf ( "%lf", &a1
);( "a2=" ); scanf ( "%lf", &a2 );
/* Открытие файлов: */( "\nFileName: " ); scanf ( "%s",
f_Name_u );( f_Name_x, f_Name_u );( f_Name_J, f_Name_u );( f_Name_u,
"_u.dat" );( f_Name_J, "_J.dat" );( f_Name_x,
"_x.dat" );( ( ru = fopen ( f_Name_u, "wt" ) ) == NULL )
{( f_Name_u );(-1);
}( ( rJ = fopen ( f_Name_J, "wt" ) ) == NULL )
{( f_Name_J );(-1);
}( ( rx = fopen ( f_Name_x, "wt" ) ) == NULL )
{( f_Name_x );(-1);
}= TBEG;
/* Печать в файл и на экран x0: */
fprintf ( rx, "%lf ", tau );( "\ntau = %lf
", tau );( "x = ( " );( j=0; j<NX; j++ )
{( rx, "%lf ", x[j] );( "%lf ", x[j] );
}( rx, "\n" ); ( ")\n"
);
/* Начальн. формирование некотор.параметров конечномерн.задачи: */=
TETA/N; /* шаг(длина перем.)*/
printf ( "h = %lf\n", h );();
/* Работа регулятора: */= 0;( !Ostanov )
{
/* Программн. упpавл. задачи с перем. структурой: */NIntU2 ( h );
/* Печ. в файл u(t), x(t) на [tau,tau+h]: */
fprintf ( ru, "%lf %lf\n", tau, u[0] + Uf );( ru,
"%lf %lf\n", tau+h, u[0] + Uf );( rJ, "%lf %lf\n", tau, J
);= tau;( "tau = %lf\n", tau );( tout < tau+h-0.000001 )
{= tout;( tau + h - tout > 0.1 )+= 0.1;= tau + h;(
"\nt = %lf", t );( " tout = %lf", tout );= 1;= 1;(inte)
{( NX, t, tout, x, f_x, &hx, ypx, f1x, f2x, f3x, f4x,
f5x,
&iflx, &nfex, &kopx, &initx, &jflx,
&kflx,
&rerr, &aerr, &savrex, &savaex);( iflx==2 )=
0;
{("x\n");(&iflx, &rerr, &aerr);
}
}( rx, "%lf ", tout );( "x = ( " );( j=0;
j<NX; j++ )
{( rx, "%lf ", x[j] );( "%lf ", x[j] );
}( rx, "\n" );( ") " );
}+= h;( "\n\n\n tau = %lf\n\n\n", tau );(1)
{( "|(x1-2)^2 + (x2-2)^2 - 1| = %lf\n",
fabs( (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0)-1.0 ) );(
"Остановить процесс ? ( 1 -- да, 0 -- нет )" );
scanf
( "%d", &Ostanov );
if ( ( Ostanov==0 ) || ( Ostanov==1 ) );
}
} /* while( !Ostanov ) */();( j=0; j<N_M; j++ )( D[j] );(
D );( j=0; j<N_M; j++ )( Ma[j] );
free ( Ma );0;
}