|
Текст подпрограммы и версий de10r_c.zip , de10d_c.zip |
Тексты тестовых примеров tde10r_c.zip , tde10d_c.zip |
Вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона.
Решается задача Коши для системы M обыкновенных дифференциальных уравнений
Y ' = F (X, Y) ,
Y = ( y1, ... , yM ) ,
F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
с начальными условиями, заданными в точке XN :
Y(XN) = YN , YN = ( y10, ... , yM0 ) ,
методом Mеpсона.
Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа Р (это число называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.
Дж.Н.Ланс, Численные методы для быстродействующих вычислительных машин, Изд-во иностранной литературы, M., 1962.
int de10r_c (S_fp f, integer *m, real *xn, real *yn, real *xk,
real *hmin, real *eps, real *p, real *h, real *y, real *ra,
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 - | минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| h - | вещественная переменная, содержащая начальное значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn, или без всякого учета в виде абсолютной величины; |
| y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk. Для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m. В случае совпадения значений параметров xn и xk значение y полагается равным начадьному значению yn (тип: вещественный); |
| ra - | одномерный рабочий массив вещественного типа длины 4*m. |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая-нибудь компонента решения не может быть вычислена с требуемой точностью eps. B этом случае интегрирование системы прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров hmin и h. |
Версии
| de10d_c - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона с повышенной точностью. При этом параметры xn, yn, xk, hmin, eps, p, h, y, ra и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de10r_c. |
| utde11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de10d_c. |
Замечания по использованию
|
Подпрограммы de10r_c и de10d_c предназначены для численного решения дифференциальных уравнений и систем уравнений с правой частью, имеющей непрерывные частные производные вплоть до 5 порядка включительно. Они являются эффективными для нежестких уравнений и систем уравнений с несложными правыми частями (т.е. не являющимися трудоемкими для вычислений). Хотя заданная точность eps не гарантируется в общем случае, большой опыт эксплуатации данной подпрограммы убедительно показывает, что вычисляемое ею численное решение достаточно близко приближает точное решение. При работе подпрограммы значения параметров m, xn, yn, xk, hmin, eps и p сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить. При работе подпрограммы счета правой части f значения параметров x, y и dy не должны изменяться. |
Использование подпрограммы иллюстрируется на примере:
y '1 = y2 ,
y '2 = -y1 , 3π/4 ≤ x ≤ π
y1(3π/4) = √2 / 2 , y2(3π/4) = - √2 / 2
Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей (основной) программы, а также результаты счета.
int main(void)
{
/* Builtin functions */
double sqrt(double);
/* Local variables */
extern int de10r_c(U_fp, int *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, int *);
static float hmin;
static int ierr;
extern int f_c();
static float h__;
static int m;
static float p, y[2], ra[8], xk, xn, eps;
m = 2;
xn = 2.3561944901925003f;
h__ = .01f;
y[0] = (float)sqrt(2.) / 2.f;
y[1] = -y[0];
xk = 3.14159265359f;
hmin = 1e-4f;
eps = 1e-5f;
p = 1e-7f;
de10r_c((U_fp)f_c, &m, &xn, y, &xk, &hmin, &eps, &p, &h__, y, ra, &ierr);
printf("\n %16.7e \n", xk);
printf("\n %16.7e %16.7e \n", y[0], y[1]);
return 0;
} /* main */
int f_c(float *x, float *y, float *z, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = y[2];
z__[2] = -y[1];
return 0;
} /* f_c */
Результаты:
y (1) = 0.4495909707 * 10-8
y (2) = -0.9999999996
ierr = 0