|
Текст подпрограммы и версий de04r_c.zip , de04d_c.zip |
Тексты тестовых примеров tde04r_c.zip , tde04d_c.zip |
Выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования линейной устойчивой системы M обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
(1) Y ' (X) = A(X) * Y(X) + φ(X) ,
Y = ( y1,..., yM ) ,
A(X) = ( ai j(X) ) , i, j = 1, ..., M ,
φ(X) = ( φ1(X), ... , φM(X) )
Предполагается, что среди характеристических корней матрицы A(X) имеются большие по модулю корни, а функция φ(x) является достаточно малой. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Вычисление производится по методу Лоусона.
Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции Y (x) на [ xn, xn + H ] по формуле
Y(x) = exp [ ( x - xn ) A0 ] Z(x),
где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X):
(2) Z ' (x) = A1(x) Z(x) + φ1(x) = exp [ - ( x - xn ) A0 ] { A(x) - A0 }
exp [ ( x - xn ) A0 ] Z(x) + exp [ ( - ( x - xn ) A0 ] φ(x)
Данное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу A0 = A (xn + H /2). Если шаг H достаточно мал, то преобразование позволяет уменьшить характеристические корни матрицы A1 (x) по сравнению с характеристическими корнями исходной матрицы A (x). Это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1). Для решения системы (2) применяются формулы классического метода Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме. Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
J.Douglas Louson. Generalized Runge - Kutta processes for stable systems with large Lipshitz constants, SIAM Journal on Numerical Analisys. - 1967. - Vol 4, No 3.
int de04r_c (S_fp fa, S_fp fi, 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 *w1, real *r1, real *r2, real *r3, real *r4,
real *r5, real *r6, real *yd, integer *ierr)
Параметры
| fa - |
подпрограмма вычисления матрицы системы A(x) в
точке x. Первый оператор подпрограммы должен
иметь вид: int fa(float *a, float *x, int *m). Здесь a - двумерный массив размера m*m, в котором помещается матрица системы, вычисленная при значении аргумента x (тип параметров a, x: вещественный); |
| fi - |
подпрограмма вычисления неоднородности правой части системы
φ (x) в любой точке x. Первый
оператор подпрограммы должен иметь вид: int fi(float *g, float *x, int *m). Здесь g - одномерный массив длины m, в который помещается неоднородность правой части системы, вычисленная при значении аргумента x (тип параметров g, x: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами x, yx, и h, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров h и/или hmin; на выходе из подпрограммы jstart полагается равным + 1; |
| 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, w1 | вещественные двумерные рабочие массивы размера m*m; |
|
r1, r2 - r3, r4, r5, r6, yd | вещественные одномерные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1 |
Версии
| de04d_c - | выполнение одного шага численного интегрирования линейной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, w1, r1, r2, r3, r4, r5, r6, yd и параметры a, g, x в подпрограммах fa и fi должны иметь тип double. |
Вызываемые подпрограммы
| ame2r_c - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de04r_c; |
| ame2d_c - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de04d_c; |
| utde20_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de04r_c; |
| utde21_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de04d_c |
|
Кроме того при работе подпрограмм de04r_c и de04d_c вызываются подпрограммы de04rd_c, de04rm_c, de04rp_c, de04rq_c и de04dd_c, de04dm_c, de04dp_c, de04dq_c соответственно. |
Замечания по использованию
|
Данная подпрограмма предназначена для интегрирования линейных систем, имеющих малую неоднородность φ (x). В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров m, hmin, eps, p сохраняются. При работе подпрограмм fa и fi значения параметров x и m не должны изменяться. При обращении к подпрограмме со значением jstart = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров xp, yp, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart. После работы подпрограммы в массиве r1 содержится значение оценки погрешности на шаге, вычисленной по правилу Рунге. Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (1). |
y1' = - ( 2 + x ) y1 /( 1 + x ) + 20 x y2 ,
y1(0) = 0 ,
y2' = -20 x y1 + ( 2 + x ) y2 /( 1 + x ) ,
y2(0) = 1 .
Приводятся подпрограммы вычисления матрицы системы и неоднородной части, фрагмент вызывающей программы, выполняющей несколько шагов из одной точки, а также результаты счета.
int main(void)
{
/* Local variables */
extern int de04r_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 *, float *, int *);
static float hmin;
static int ierr;
static float a[4] /* was [2][2] */, h__;
static int m;
static float p, x, e1[4] /* was [2][2] */, e2[4] /* was [2][2] */,
e3[4] /* was [2][2] */, r1[2], r2[2], r3[2], r4[2],
r5[2], w1[4] /* was [2][2] */, r6[2];
extern int fa_c(), fi_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 = 0.f;
yx[0] = 0.f;
yx[1] = 1.f;
hmin = 1e-10f;
eps = 1e-7f;
p = 100.f;
jstart = 1;
h__ = .01f;
ih = 0;
l1000:
++ih;
de04r_c((U_fp)fa_c, (U_fp)fi_c, &m, &jstart, &hmin, &eps, &p, yx, &x,
&h__, &bul, &xp, yp, a, e1, e2, e3, w1, r1, r2, r3, r4, r5, r6,
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;
case 7: goto l107;
}
l101:
eps = 1e-7f;
goto l1000;
l102:
eps = 1.f;
goto l1000;
l103:
jstart = -1;
h__ = .01f;
eps = 1e-5f;
goto l1000;
l104:
jstart = -1;
h__ = -1e-4f;
eps = .01f;
goto l1000;
l105:
h__ = 1.f;
eps = 1e-7f;
goto l1000;
l106:
eps = 1e-4f;
goto l1000;
l107:
return 0;
} /* main */
int fa_c(float *a, float *x, int *m)
{
#define a_ref(a_1,a_2) a[(a_2)*2 + a_1]
/* Parameter adjustments */
a -= 3;
/* Function Body */
a_ref(1, 1) = -(*x + 2.f) / (*x + 1.f);
a_ref(1, 2) = *x * 20.f;
a_ref(2, 1) = -a_ref(1, 2);
a_ref(2, 2) = a_ref(1, 1);
return 0;
} /* fa_c */
int fi_c(float *r1, float *x, int *m)
{
/* Parameter adjustments */
--r1;
/* Function Body */
r1[1] = 0.f;
r1[2] = 0.f;
return 0;
} /* fi_c */
Результаты:
После первого обращения к подпрограмме
x yx(1) yx(2)
1.000000000001-02 9.802471967628-04 9.802468700200-01
h__ r1(1) r1(2)
2.000000000001-02 0.000000000000+00 6.063298011814-14
после второго обращения к подпрограмме
x yx(1) yx(2)
3.000000000000-02 8.479506692396-03 9.421419716109-01
h__ r1(1) r1(2)
4.000000000002-02 1.231607408649-14 1.333925562599-12
после третьего обращения к подпрограмме
x yx(1) yx(2)
6.999999999994-02 4.268132414620-02 8.703501915807-01
h__ r1(1) r1(2)
8.000000000004-02 1.970571853839-12 3.565219230945-11
после четвертого обращения к подпрограмме
x yx(1) yx(2)
3.999999999996-02 1.478074532271-02 9.237177506848-01
h__ r1(1) r1(2)
2.000000000001-02 0.000000000000+00 6.063298011814-14
после пятого обращения к подпрограмме
x yx(1) yx(2)
2.990000000000-02 8.424732657545-03 9.423281849631-01
h__ r1(1) r1(2)
-2.000000000004-04 0.000000000000+00 6.063298011814-14
после шестого обращения к подпрограмме
x yx(1) yx(2)
6.115000000000-02 3.314040588361-02 8.858545422663-01
h__ r1(1) r1(2)
3.125000000000-02 4.395891058563-13 1.097456940137-11
после седьмого обращения к подпрограмме
x yx(1) yx(2)
9.240000000000-02 7.117143016819-02 8.315812963056-01
h__ r1(1) r1(2)
6.250000000000-02 7.80649190203-13 8.852415097246-12