|
Текст подпрограммы и версий de00r_c.zip , de00d_c.zip , de84r_c.zip , de84d_c.zip |
Тексты тестовых примеров tde00r_c.zip , tde00d_c.zip , tde84r_c.zip , tde84d_c.zip |
Выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = F(X, Y) ,
Y = ( y1,..., yM ) , F = ( f1( X, y1,..., yM ), ... , fM( X, y1,..., yM ) )
Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Предполагается также, что система (1) на [ xn, xn + H ] может быть достаточно хорошо аппроксимирована линейной однородной системой с постоянными коэффициентами, т.е. системой вида
Y ' = AY ,
где A - некоторая постоянная матрица. Вычисление производится по методу Лоусона. Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции y (x) на [ xn, xn + H ] по формуле
y(x) = exp [ ( x - xn ) A0 ] Z(x) ,
где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X) :
(2) Z ' (x) = F1(x, z) = exp [ - ( x - xn ) A0 ] { F(x, exp [ ( x - xn ) A0 ] Z(x)) -
- A0 exp [ ( x - xn ) A0 ] Z(x) }
Матрица Якоби ∂ F1 / ∂Z системы (2) и матрица Якоби ∂F / ∂y системы (1) связаны между собой соотношением
∂ F1 / ∂Z = exp [ - (x - xn) A0] { ∂F / ∂y - A0 } exp [ ( x - xn ) A0 ]
Указанное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу ∂F / ∂y , вычисленную либо в точке ( xn, yn ), либо в некоторой предыдущей точке ( xn - k, yn - k ). Если шаг H достаточно мал, то данное преобразование позволяет уменьшить характеристические корни матрицы ∂ F1 / ∂Z по сравнению с характеристическими корнями матрицы ∂F / ∂y . А это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1).
Для решения системы (2) применяются классический метод Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).
Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
J.Douglas Lowson, Generalized Runge - Kutta processes for stable systems with large lipshitz constants, SIAM Journal on Numerical Analisys, 1967, Vol.4, No 3.
int de00r_c (real *f, S_fp fj, integer *m, integer *jstart,
real *hmin, real *eps, real *p, real *yx, real *x, real *h,
logical *bul, real *xp, real *yp, real *a, real *e1, real *e2,
real *e3, real *e4, real *r, real *r0, real *r1, real *r2,
real *r3, real *yd, integer *ierr)
Параметры
| f - |
подпрограмма вычисления правой части системы в
точке x . Первый оператор подпрограммы должен
иметь вид: int f (float *x, float *y, float *dy, int *m). Здесь x и y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. m ≠ 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y и dy: вещественный); |
| fj - |
подпрограмма вычисления матрицы Якоби правой
части системы. Первый оператор подпрограммы должен иметь вид: int fj (float *x, float *y, float *df, int *m). Здесь x и y - значения независимой и зависимой переменных, соответственно, причем y представляет одномерный массив длины m. df - двумерный массив размера m*m. Вычисленное значение матрицы Якоби ∂F / ∂y должно быть помещено в массив df, при этом частная производная от правой части i - ого уравнения F по j - ой переменной y (j) запоминается в dF (i, j), т.е. dF (i, j) = ∂Fi / ∂yj (тип параметров x, y и df: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart; при этом значение параметра в качестве матрицы A0 принимается матрица ∂F / ∂y в точке [ x = xn, y = yn ]; вычисленное значение матрицы A0 будет использоваться в дальнейшем для значений jstart = 1, 2, 3; |
| 1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, заданных параметрами yx и x соответственно; при этом само значение величины шага интегрирования H, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней, а матрица A0 остается той же, что и на предыдущем шаге; |
| 2 - | то же, что и для jstart = 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме; матрица A0 остается той же, что и на предыдущем шаге; |
| 3 - | то же, что и для jstart = 1, но только с той разницей, что шаг интегрирования может принимать произвольное значение, матрица A0 остается той же, что и на предыдущем шаге. |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров h и/или hmin; при этом значении jstart в качестве матрицы A0 принимается матрица ∂F / ∂y в точке [ x = xp, y = yp ] таким образом, до тех пор, пока jstart принимает значения 1, 2 или 3, сохраняется одно и то же значение матрицы A0 , равное значению матрицы Якоби, вычисленному при последнем обращении к подпрограмме с параметром jstart = 0 или jstart = - 1; если требуется изменить матрицу A0 , то следует обратиться к подпрограмме со значением jstart = 0 или - 1; при выходе из подпрограммы параметр jstart полагается равным 2, если выполнив заданный в h шаг, подпрограмма рекомендовала для интегрирования вдвое большеe его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной h; |
| hmin - | минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| yx, x - | заданные вещественные значения решения и соответствующее ему значение аргумента; в результате работы подпрограммы в x получается новое значение аргумента, в yx - соответствующее значение решения; в случае системы уравнений т.е. когда m ≠ 1, yx задается одномерным массивом длины m; |
| h - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность eps; на выходе из подпрограммы h содержит рекомендуемое подпрограммой значение следующего шага, определяемое ею с целью достижения более экономного способа интегрирования; |
| bul - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в H шаг выводит в конец интервала интегрирования, и FALSE_ в противном случае; в результате работы подпрограммы bul равно FALSE_, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в h шаг, значение параметра bul не меняется; |
| xp, yp - | вещественная рабочая переменная и одномерный рабочий массив длины m соответственно; значения параметров xp, yp на выходе из подпрограммы равны тем значениям, которые имели параметры x, yx при входе в нее (т.е. предыдущий узел и решение в нем); |
|
a, e1, e2 - e3, e4 | вещественные двумерные рабочие массивы размера m*m; |
|
r, r0 - r1, r2 r3, yd | вещественные одномерные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1. |
Версии
| de00d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме f и x, y, df в подпрограмме fj должны иметь тип double; |
| de84r_c - |
выполнение одного шага численного интегрирования
устойчивой системы обыкновенных дифференциальных
уравнений первого порядка с большой константой Липшица
методом Лоусона. В отличии от подпрограммы de00r_c
первый оператор подпрограммы de84r_c не содержит в
списке формальных параметров имени подпрограммы
вычисления матрицы Якоби и имеет следующий вид: int de84r_c (real *f, integer *m, integer *jstart, real *hmin, real *eps, real *p, real *yx, real *x, real *h, logical *bul, real *xp, real *yp, real *a, real *e1, real *e2, real *e3, real *e4, real *r, real *r0, real *r1, real *r2, real *r3, real *yd, integer *ierr). При использовании этой подпрограммы пользователь не составляет подпрограмму вычисления матрицы Якоби. Все параметры подпрограммы de84r_c имеют тот же смысл, что и одноименные параметры подпрограммы de00r_c; |
| de84d_c - | выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de84r_c, при этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| ame2r_c - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de00r_c; |
| ame2d_c - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de00d_c; |
| utde20_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de00r_c; |
| utde21_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de00d_c; |
| Кроме того, при работе подпрограмм de00r_c, de84r_c и de00d_c, de84d_c вызываются рабочие подпрограммы de00rs_c, de04rm_c, de00rd_c и de00ds_c, de04dm_c, de00dd_c соответственно. | |
| При работе подпрограмм de84r_c и de84d_c дополнительно вызываются также рабочие подпрограммы de84rj_c и de84dj_c соответственно. |
Замечания по использованию
|
Данная подпрограмма предназначена для интегрирования таких
нелинейных систем, которые могут быть достаточно хорошо
аппроксимированы линейной однородной системой с постоянными
коэффициентами, т.е. системой вида Y ' = AY , где A - некоторая постоянная матрица. В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров m, hmin, eps, p сохраняются. Между последовательными обращениями к подпрограмме со значениями параметра jstart = 1, 2, 3 пользователь не должен изменять содержимое массивов a, e1, e2. При работе подпрограмм f и fj значения параметров m, x, y не должны изменяться. При обращении к подпрограмме со значением jstart = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров xp и yp соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart. После работы подпрограммы в массиве r1 содержится значение оценки абсолютной погрешности на шаге, вычисленной по правилу Рунге. Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (1). |
y1' = (- 500 + 2 x) y1 + (1 - x / 10 )y2 + (1 + 1 / ( x+1 ) ) sin y1 ,
y1( = 49
y2' = (- 1000 - 10 x) y1 + (1 + x / 10 ) y2 + (2 + 1 / ( x+2 ) ) cos y2 ,
y2( = 50
1 ≤ x ≤ 6
Приводятся подпрограммы вычисления правой части системы и матрицы Якоби правой части, а также фрагмент вызывающей программы и результаты счета.
int main(void)
{
/* Local variables */
extern int de00r_c(U_fp, U_fp, int *, int *, float *, float *,
float *, float *, float *, float *, logical *,
float *, float *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, int *);
static float hmin;
static int ierr;
static float a[4] /* was [2][2] */;
extern int f_c();
static float h__;
static int m;
static float p, r__[2], x, e1[4] /* was [2][2] */, e2[4] /* was [2][2] */,
e3[4] /* was [2][2] */, e4[4] /* was [2][2] */,
r0[2], r1[2], r2[2], r3[2];
extern int fj_c();
static int ih;
static float yd[2], xp, yp[2], yx[2];
static int jstart;
static logical bul;
static float eps;
m = 2;
x = 1.f;
yx[0] = 49.f;
yx[1] = 50.f;
hmin = 1e-10f;
eps = 1e-8f;
p = 100.f;
jstart = 0;
h__ = .01f;
ih = 0;
l100:
++ih;
de00r_c((U_fp)f_c, (U_fp)fj_c, &m, &jstart, &hmin, &eps, &p, yx, &x, &h__,
&bul, &xp, yp, a, e1, e2, e3, e4, r__, r0, r1, r2, r3, yd,
&ierr);
printf("\n %16.7e %16.7e %16.7e \n", x, yx[0], yx[1]);
printf("\n %16.7e %16.7e %16.7e \n\n", h__, r1[0], r1[1]);
switch (ih) {
case 1: goto l101;
case 2: goto l102;
case 3: goto l103;
case 4: goto l104;
case 5: goto l105;
case 6: goto l106;
}
l101:
h__ = 1.f;
jstart = 3;
goto l100;
l102:
eps = .01f;
goto l100;
l103:
jstart = -1;
eps = 1e-8f;
goto l100;
l104:
eps = .01f;
h__ = 1e-8f;
jstart = 3;
goto l100;
l105:
h__ = -.001f;
eps = 1e-8f;
jstart = 0;
goto l100;
l106:
return 0;
} /* main */
int f_c(float *x, float *y, float *dy, int *m)
{
/* Builtin functions */
double sin(double), cos(double);
/* Parameter adjustments */
--dy;
--y;
/* Function Body */
dy[1] = (*x * 2.f - 500.f) * y[1] + (1.f - *x * .1f) * y[2] +
(1.f / (*x + 1.f) + 1.f) * (float)sin(y[1]);
dy[2] = (-1e3f - *x * 10.f) * y[1] + (*x * .1f + 1.f) * y[2] +
(1.f / (*x + 2.f) + 2.f) * (float)cos(y[2]);
return 0;
} /* f_c */
int fj_c(float *x, float *y, float *df, int *m)
{
/* Builtin functions */
double cos(double), sin(double);
#define df_ref(a_1,a_2) df[(a_2)*2 + a_1]
/* Parameter adjustments */
df -= 3;
--y;
/* Function Body */
df_ref(1, 1) = *x * 2.f - 500.f +
(1.f / (*x + 1.f) + 1.f) * (float)cos(y[1]);
df_ref(1, 2) = 1.f - *x * .1f;
df_ref(2, 1) = -1e3f - *x * 10.f;
df_ref(2, 2) = *x * .1f + 1.f -
(1.f / (*x + 2.f) + 2.f) * (float)sin(y[2]);
return 0;
} /* fj_c */
Результаты:
После первого обращения к подпрограмме
x yx(1) yx(2)
1.000039062499+00 4.895766071531+00 4.808757932455+00
h__ r1(1) r1(2)
3.906250000002-05 9.597251282538-10 -7.233271996176-09
После второго обращения к подпрограмме
x yx(1) yx(2)
1.000069580077+00 4.733408100787+01 4.661901490984+01
h__ r1(1) r1(2)
3.051757812500-05 1.319373647371-10 6.946114202335-02
После третьего обращения к подпрограмме
x yx(1) yx(2)
1.000100097655+00 4.662140584405+01 4.517257692956+01
h__ r1(1) r1(2)
6.103515625000-05 -6.208817164108-11 1.830362084864-09
После четвертого обращения к подпрограмме
x yx(1) yx(2)
1.000100097655+00 4.662140584411+01 4.517257692973+01
h__ r1(1) r1(2)
3.051757812500-05 -5.820766091346-11 1.851003617047-09
После пятого обращения к подпрограмме
x yx(1) yx(2)
1.000100107654+00 4.662117408356+01 4.517210655950+01
h__ r1(1) r1(2)
2.000000000004-08 0.000000000000 + 00 0.000000000000 + 00
После шестого обращения к подпрограмме
x yx(1) yx(2)
1.000068857653+00 4.735108317091+01 4.665352269082+01
h__ r1(1) r1(2)
-3.125000000004-05 3.880510727560-11 -2.250696221985-09