|
Текст подпрограммы и версий de43r_c.zip , de43d_c.zip , de49r_c.zip , de49d_c.zip |
Тексты тестовых примеров tde43r_c.zip , tde43d_c.zip , tde49r_c.zip , tde49d_c.zip |
Вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений второго порядка в конце интервала интегрирования методом Штермера с контролем точности.
Решается задача Коши для системы M обыкновенных дифференциальных уравнений второго порядка
Y '' = F (X, Y) ,
Y = ( y1, ..., yM ) ,
F = ( f1 ( X, y1, ..., yM ), ..., fM ( X, y1, ..., yM ) )
с начальными условиями, заданными в точке XN:
Y (XN) = YN , YN = ( y10, ..., yM0 ) ,
Y ' (XN) = DYN , DYN = ( y10', ..., yM0' ) , -
многошаговым предсказывающе - исправляющим методом пятого порядка точности.
Суть метода состоит в том, что сначала по явной формуле Адамса - Штермера строятся предсказанные значения приближенного решения, которые затем исправляются по неявной формуле Адамса - Штермера.
Решение вычисляется в одной точке XK, которая является концом интервала интегрирования.
Bсе компоненты решения проверяются на точность по меpе погрешности, т.е. если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Березин И.С., Жидков Н.П. Методы вычислений, т.2, Физматгиз, М., 1960.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб."Вычислительные методы и программирование", вып.1, Изд - во МГУ, 1962.
int de43r_c (real *f, integer *m, real *xn, real *yn,
real *dyn, real *xk, real *hmin, real *hmax, real *eps, real *p,
real *h, real *y, real *dy, real *delty, real *df, real *rfn,
real *rf, real *yp, real *dyp, 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 - | количество уравнений в системе (тип: целый); |
|
xn, yn - dyn | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. m ≠ 1) yn и dyn представляют одномерные массивы длины m (тип вещественный); |
| xk - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или pавно xn (тип: вещественный); |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| hmax - | максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| h - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn, или без такого учета в виде абсолютной величины; |
| y - | искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента xk; для системы уравнений задается одномерным массивом длины m; в случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный); |
|
dy - delty rfn, rf yp, dyp | одномерные вещественные рабочие массивы длины m; |
| df - | двумерный вещественный рабочий массив размеpа m * 5; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| ierr=65 | и ierr=66 - когда какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps при заданных начальном шаге h и его минимальном значении hmin, при этом ierr = 66 указывает, что тpебуемая точность не может быть достигнута при разгоне, а ierr=65 - после разгона; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров h и hmin. |
Версии
| de43d_c - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений второго порядка в конце интервала интегрирования методом Штермера с повышенной точностью. При этом параметры xn, yn, dyn, xk, hmin, hmax, eps, p, h, y, dy, delty, df, rfn, rf, yp, dyp и параметры x, y и dy в подпрограмме f должны иметь тип double. |
| de49r_c - | назначение этой подпрограммы такое же, как и подпрограммы de43r_c, но схема организации вычислений в ней несколько иная, чем в de43r_c, а именно: во - первых, исправленное и предсказанное значения решения вычисляются через значения некоторой промежуточной переменной, введение которой позволяет уменьшить вычислительную погрешность в приближенном значении решения; во - вторых, для вычисления приближенного решения на одном шаге используется на одно вычисление правой части уравнения больше, чем в подпрограмме de43r_c. Таким образом, обеспечивается возможность вычислить решение уравнения более точно, чем в de43r_c. |
| de49d_c - | то же, что и de49r_c, но все вычисления ведутся с удвоенным числом значащих цифр; при этом параметры xn, yn, dyn, xk, hmin, hmax, eps, p, h, y, dy, delty, df, rfn, rf, yp, dyp и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| de42r_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с контролем точности. |
| de42d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с повышенной точностью. |
| de48r_c - de48d_c | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера, преобразованным к виду, где влияние вычислительной погрешности меньше. |
|
utde12_c - utde13_c utde16_c utde17_c | подпрограммы выдачи диагностических сообщений. |
|
Подпрограммы de42r_c и utde12_c вызываются при работе подпрограммы de43r_c, а подпрограммы de42d_c и utde13_c - при работе подпрограммы de43d_c. Подпрограммы de48r_c и utde16_c вызываются при работе подпрограммы de49r_c, а подпрограммы de48d_c и utde17_c - при работе подпрограммы de49d_c. |
Замечания по использованию
|
B общем случае заданная точность не гарантируется. При работе подпрограммы и ее версии значения параметров m, xn, yn, dyn, xk, hmin, hmax, eps, p сохраняются. Если после работы подпрограммы нет необходимости иметь начальные значения решения yn и / или его производной dyn, то параметры yn и y, а в случае производной параметры dyn и dy, при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением ierr = 65 или 66, значения параметров yn и dyn будут испорчены. B том случае, когда важнее сократить время вычисления, целесообразно использовать подпрограмму de43r_c; если требуется поточнее вычислить решение, то следует использовать подпрограмму de49r_c. Значение параметра hmax не должно быть слишком большим, т.к. в противном случае величина шага интегрирования может достичь таких размеров, которые приведут к абсолютной неустойчивости численного метода. |
y1'' = - 6y1 - 7y2 ,
y2'' = - 3y1 - 2y2 + 2x ;
y1 (0) = 0 , y2 (0) = 0 ,
y1' (0) = - 1 / 9 , y2' (0) = - 3 .
Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей программы, выполняющей дважды интегрирование данной системы, сначала на отрезке [0, 5], затем на отрезке [-5, 0], справа налево, а также результаты счета.
int main(void)
{
/* Builtin functions */
double exp(double), sin(double);
/* Local variables */
extern int de43r_c(U_fp, int *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, int *);
static float hmin, hmax;
static int ierr;
extern int f_c();
static float h__;
static int k, m;
static float p, y[2], delty[2], y1, y2, y3, df[10], rf[2], dy[2], xk, xn,
yn[2], yp[2], rfn[2], eps, dyn[2], dyp[2];
h__ = .01f;
hmin = 1e-10f;
m = 2;
xn = 0.f;
xk = 5.f;
yn[0] = 0.f;
yn[1] = 0.f;
dyn[0] = -.1111111111111111f;
dyn[1] = -3.f;
eps = 1e-4f;
y3 = ((float)exp(5.f) - (float)exp(-5.f)) / 3.f;
y1 = y3 - (float)sin(15.f) * 7.f / 9.f + 7.7777777777777777f;
y2 = -y3 - (float)sin(15.f) / 3.f - 6.666666666666667f;
p = 100.f;
hmax = .1f;
for (k = 1; k <= 2; ++k) {
de43r_c((U_fp)f_c, &m, &xn, yn, dyn, &xk, &hmin, &hmax, &eps, &p,
&h__, y, dy, delty, df, rfn, rf, yp, dyp, &ierr);
printf("\n %16.7e \n", xk);
printf("\n %16.7e %16.7e \n", y[0], y[1]);
printf("\n %16.7e \n", h__);
printf("\n %16.7e %16.7e \n", y1, y2);
if (ierr != 0) {
goto L20;
}
y1 = -y1;
y2 = -y2;
xk = -5.f;
h__ = .01f;
/* L200: */
}
L20:
return 0;
} /* main */
int f_c(float *t, float *y, float *z, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = y[1] * -6. - y[2] * 7.;
z__[2] = y[1] * -3. - y[2] * 2. + *t * 2.;
return 0;
} /* f_c */
Результаты:
после первого обращения к подпрограмме -
xk y(1) y(2)
5.000000000 + 00 5.674103844 + 01 - 5.655217994 + 01
h__ y1 y2
1.000000000 - 01 5.674080540 + 01 - 5.635223633 + 01
после второго обращения к подпрограмме -
xk y(1) y(2)
- 5.000000000 + 00 - 5.674103838 + 01 5.635217988 + 01
h__ y1 y2
- 1.000000000 - 01 - 5.674080540 + 01 5.635223633 + 01