|
Текст подпрограммы и версий de08r_c.zip , de08d_c.zip |
Тексты тестовых примеров tde08r_c.zip , tde08d_c.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Инглэнда.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
Y ' = F(X, Y) ,
Y = ( y1, ... , yM ) , F = ( f1(X, y1,..., yM), ... , fM(X, y1,..., yM) )
методом Инглэнда пятого порядка точности. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Все компоненты решения вычисляются с контролем точности по мере погрешности следующим образом.
Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. В качестве абсолютной погрешности решения принимается оценка главного члена асимптотического разложения погрешности метода на шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на одном шаге интегрирования для определения решения и погрешности используются всего шесть вычислений правой части системы.
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
| 1. | Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987. |
| 2. | England R. Error estimates for Runge - Kutta type solutions to systems of ordinary differential equations. The Computer Journal. 1969 - v.12, No.2. |
int de08r_c (S_fp f, integer *m, integer *jstart, real *hmin,
real *eps, real *p, real *yx, real *x, real *h, logical *bul,
real *xp, real *yp, real *r1, real *r2, real *r3, integer *ierr)
Параметры
| f - |
имя подпрограммы вычисления значений правой части
системы дифференциальных уравнений. Первый оператор
подпрограммы должен иметь вид int f (float *x, float *y, float *dy, int *m). Здесь x, y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в dy. В случае ситемы уравнений, т.е. когда m ≠ 1, параметры y и dy представляют одномерные массивы длины m (тип параметров x, y и dy: вещественный); |
| 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, соответственно, при входе в нее (т.е. предыдущий узел и решение в нем); |
| r1 - | двумерный вещественный рабочий массив размера m*5; |
| r2, r3 - | одномерные вещественные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1 |
Версии
| de08d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Инглэнда с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, r1, r2, r3 и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde20_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de08r_c. |
| utde21_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de08d_c. |
Замечания по использованию
|
В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров m, hmin, eps, p сохраняются. При работе подпрограммы счета правой части f значения параметров x, y и m не должны изменяться. При обращении к подпрограмме со значением jstart = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров xp и yp, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart. На выходе из подпрограммы в массиве r3 содержится оценка локальной погрешности метода на выполненном шаге. |
y1' = - y1 - 5y2 , y1(0) = 1
y2' = y1 + y2 , y2(0) = 1
Точное решение системы:
y1 = cos2x - 3 sin2x ,
y2 = cos2x + sin2x
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета.
int main(void)
{
/* Builtin functions */
double sin(double), cos(double);
/* Local variables */
extern int de08r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, logical *, float *,
float *, float *, float *, float *, int *);
static float hmin;
static int ierr;
static float h__;
static int m;
static float p, x;
extern int fengl_c();
static float r1[10] /* was [2][5] */, r2[2], r3[2], y1, y2;
static int ih;
static float xp, yp[2], yx[2];
static int jstart;
static logical bul;
static float eps;
m = 2;
p = 100.f;
hmin = 1e-12f;
h__ = .01f;
eps = 1e-5f;
jstart = 1;
x = 0.f;
yx[0] = 1.f;
yx[1] = 1.f;
ih = 0;
bul = FALSE_;
l20:
de08r_c((U_fp)fengl_c, &m, &jstart, &hmin, &eps, &p, yx, &x, &h__, &bul,
&xp, yp, r1, r2, r3, &ierr);
y1 = (float)cos((float)(x * 2.f)) - (float)sin((float)(x * 2.f)) * 3.f;
y2 = (float)cos((float)(x * 2.f)) + (float)sin((float)(x * 2.f));
++ih;
printf("\n %15.7e %15.7e %15.7e \n", x, yx[0], yx[1]);
printf("\n %15.7e %15.7e %15.7e \n", h__, y1, y2);
printf("\n %15.7e %15.7e \n\n", r3[0], r3[1]);
switch (ih) {
case 1: goto l20;
case 2: goto l20;
case 3: goto l103;
case 4: goto l104;
case 5: goto l105;
case 6: goto l20;
case 7: goto l100;
}
l103:
jstart = -1;
h__ = -.005f;
goto l20;
l104:
jstart = -1;
h__ = .02f;
goto l20;
l105:
eps = 1e-8f;
h__ = .01f;
goto l20;
l100:
return 0;
} /* main */
int fengl_c(float *x, float *y, float *z, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = -y[1] - y[2] * 5.f;
z__[2] = y[1] + y[2];
return 0;
} /* fengl_c */
Результаты:
после первого обращения к подпрограмме -
x yx(1) yx(2)
1.000000000001-02 9.398040065844-01 1.019798673355 + 00
h__ y1 y2
1.000000000001-02 9.398040065835-01 1.019798673356 + 00
контрольный член: -7.983658179000-11 2.673061771929-11
после второго обращения к подпрограмме -
x yx(1) yx(2)
2.000000000001-02 8.792321040974-01 1.039189440839 + 00
h__ y1 y2
1.000000000001-02 8.792321040992-01 1.039189440846 + 00
контрольный член: -8.037659426918-11 2.623323780426-11
после третьего обращения к подпрограмме -
x yx(1) yx(2)
3.00000000000-02 8.183085204928-01 1.058164546403 + 00
h__ y1 y2
1.00000000001-02 8.183085204937-01 1.058164546412 + 00
контрольный член: -8.088818503892-11 2.565069276093-11
после четвертого обращения к подпрограмме -
x yx(1) yx(2)
1.500000000000-02 9.095635331360-01 1.029545533938 + 00
h__ y1 y2
-1.000000000001-02 9.095635331387-01 1.029545533947 + 00
контрольный член: 2.458477865729-12 -8.100187187665-13
после пятого обращения к подпрограмме -
x yx(1) yx(2)
3.000000000000-02 8.183085204928-01 1.058164546403 + 00
h__ y1 y2
1.000000000001-02 8.183085204937-01 1.058164546412 + 00
контрольный член: -8.088818503892-11 2.565059276093-11
после шестого обращения к подпрограмме -
x yx(1) yx(2)
3.249999999997-02 8.030255271615-01 1.062842482492+00
h__ y1 y2
2.500000000001-03 8.030255271660-01 1.062842482504+00
контрольный член: -1.136868377216-13 1.953992523340-14
после седьмого обращения к подпрограмме -
x yx(1) yx(2)
3.499999999997-02 7.877224582326-01 1.067493847573+00
h__ y1 y2
2.500000000001-03 7.877224582389-01 1.067493847588+00
контрольный член: -9.237055564881-14 1.065814103640-14