|
Текст подпрограммы и версий de16r_c.zip , de16d_c.zip |
Тексты тестовых примеров tde16r_c.zip , tde16d_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 ) ,
методом Рунге - Кутта - Фельберга пятого порядка точности. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа Р, и по абсолютной погрешности на остальных участках, т.е. там, где модуль компоненты меньше этого числа. B качестве абсолютной погрешности решения используется оценка главного члена асимптотического разложения погрешности метода на одном шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на каждом шаге интегрирования для определения решения и его погрешности используются всего шесть вычислений правой части системы.
Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.
int de16r_c (real *f, integer *m, real *xn, real *yn, real *xk,
real *hmin, real *eps, real *p, real *h, real *y, real *rab,
integer *ierr)
Параметры
| f - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператоp подпрограммы должен иметь вид: int f (float *x, float *y, float *dy, int *m). Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. B случае системы уравнений, т.е. когда m ≠ 1 , параметры y и dy представляют массивы длины m (тип параметров x, y и dy: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; в случае системы уравнений (т.е. m ≠ 1) yn представляет одномерный массив длины m (тип: вещественный); |
| xk - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или pавно xn (тип: вещественный); |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| h - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если xn < xk, отрицательным, если xn > xk, или без всякого учета в виде абсолютной величины; на выходе из подпрограммы содержит значение последнего шага интегрирования; |
| y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk; для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m; в случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный); |
| rab - | одномерный рабочий массив вещественного типа длиной 6*m + 1; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps. B этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров hmin и h. |
Версии
| de16d_c - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Рунге - Кутта - Фельберга с удвоенным числом значащих цифр. При этом параметры xn, yn, xk, hmin, eps, p, h, y, rab и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
|
de15r_c - de15d_c | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Рунге - Кутта - Фельберга пятого порядка точности. |
|
utde16_c - utde17_c | подпрограммы выдачи диагностических сообщений. |
| Подпрограммы de15r_c и utde16_c вызываются при работе подпрограммы de16r_c, а подпрограммы de15d_c и utde17_c - при работе подпрограммы de16d_c. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров m, xn, yn, xk, hmin, eps, p сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить. При работе подпрограммы счета правой части f значения параметров x, y и m не должны изменяться. |
y1' = -y2 - 0.1 x - 0.9
y2' = y1 - 0.1 x - 1.1 0 ≤ x ≤ π
y1 (0) = 1 , y2 (0) = - 2
Точное решение системы:
y1 = 0.1 x + sin x + 1.
y2 = - 0.1 x - cos x - 1.
Приводятся подпрограммы вычисления значений правой части и фрагмент вызывающей программы, а также результаты счета.
int main(void)
{
/* Builtin functions */
double sin(double), cos(double);
/* Local variables */
extern int de16r_c(U_fp, int *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, int *);
static float hmin;
static int ierr;
static float h__;
static int m;
extern int ffelb_c();
static float p, y[2], y1, y2, xk, xn, yn[2], rab[13], eps;
m = 2;
xn = 0.f;
yn[0] = 1.f;
yn[1] = -2.f;
xk = 3.141592653589793238f;
hmin = 1e-14f;
eps = 1e-5f;
p = 100.f;
h__ = .01f;
de16r_c((U_fp)ffelb_c, &m, &xn, yn, &xk, &hmin, &eps, &p, &h__, y, rab,
&ierr);
printf("\n %5i \n", ierr);
y1 = xk * .1f + (float)sin(xk) + 1.f;
y2 = xk * -.1f - (float)cos(xk) - 1.f;
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);
return 0;
} /* main */
int ffelb_c(float *x, float *y, float *z, int *m)
{
static float r__;
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
r__ = *x * .1f;
z__[1] = -y[2] - r__ - .9f;
z__[2] = y[1] - r__ - 1.1f;
return 0;
} /* ffelb_c */
Результаты:
ierr = 0
y(1) = 1.314159265929 + 00
y(2) = -3.141592563520 - 01