Текст подпрограммы и версий de23r_c.zip , de23d_c.zip , de25r_c.zip , de25d_c.zip |
Тексты тестовых примеров tde23r_c.zip , tde23d_c.zip , tde25r_c.zip , tde25d_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 случае, когда система уравнений является жесткой, интегрирование осуществляется специальным методом, основанном на методе типа Адамса и использующим якобиан ( ∂F/∂Y ) системы, который вычисляется подпрограммой по формулам численного дифференцирования. При интегрировании данной системы уравнений численное решение проверяется на точность; считается что значение решения в узле xn вычислено с требуемой точностью ЕРS, если выполняется следующее соотношение:
M ( ∑ δI2 )1/2 ≤ EPS , I=1
где δI - одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I - й компоненты решения в узле xn подразумевается отношение абсолютной погрешности eI этого значения в узле xn к абсолютной величине значения I - й компоненты в предыдущем узле xn - 1, т.е. eI / | yIn - 1|, а под стандартной погрешностью - отношение eI / YPM (I), где
YPM(I) = max ( | yI0 |, | yI1 |,..., | yIn-1 | )
Tип погрешности специфицируется пользователем при обращении к подпрограмме.
Gear C.W. The automatic integration of ordinary differential equations. Communicatuons of the ACM, 14, 3 (March 1971), 176-179.
Gear C.W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice - Hall, Englewood Cliffs, N.J., 1971.
Gear C.W., The automatic integration of stiff ordinary differential equations. Information Processing 68, A.J.H.
int de23r_c (real *f, integer *m, real *xn, real *yn, real *xk, real *hmin, real *hmax, real *eps, integer *istifj, integer *iorder, integer *iu, real *h, real *y, real *ypm, real *delty, real *rab, real *yp, 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ое разрешается использовать при интегрировании данной системы уравнений; это значение должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром h (тип: вещественный); |
hmax - | максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
eps - | допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметpа iu (тип: вещественный); |
istifj - | целый указатель метода численного интегрирования: |
istifj=0 - | интегрирование системы ведется методом Адамса; |
istifj=1 - | интегрирование ведется специальным методом, предназначенным для жестких систем; |
iorder - | целая переменная, указывающая максимальный допустимый порядок метода; iorder должен быть не больше 7 для метода Адамса и не больше 6 для метода интегрирования жестких систем; |
iu - | целый указатель типа погрешности численного решения: |
iu = 1 - | для стандартной погрешности; |
iu = 2 - | для относительной погрешности; |
iu = 3 - | для абсолютной погрешности; |
h - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления итегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn, или без такого учета в виде абсолютной величины; |
y - | искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента xk; для системы уравнений (когда m ≠ 1) задается массивом длиной m; в случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный); |
ypm - delty | одномерные вещественные рабочие массивы длиной m; |
rab - | одномерный вещественный рабочий массив; при интегрировании нежесткой системы уравнений rab имеет размер 17*m, при интегрировании жесткой системы - m*(m + 17); |
yp - | двумерный вещественный рабочий массив размеpа 8*m; |
ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
ierr= 1 - | когда неправильно задан параметр iorder, а именно, когда iorder превосходит максимальный допустимый порядок метода; в этом случае интегрирование системы ведется методом Гира порядка не выше 7 для нежесткой системы, и не выше 6 для жесткой; |
ierr=65 - | когда решение системы не может быть вычислено с требуемой точностью eps при заданных начальном шаге h, его минимальном значении hmin и порядке метода iorder; |
ierr=66 - | когда приближенное значение решения не может быть вычислено, т.к. итерационный процесс его определения не сходится для шагов интегрирования h, больших заданного минимального значения hmin; |
ierr=67 - | когда требуемая точность eps вычисления приближенного решения меньше той, которая может быть достигнута для данной задачи при тех размерах шага интегрирования, начальное значение которого задано параметром h; |
при ierr = 65, 66, 67 интегрирование системы прекращается; при желании интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и iorder; |
ierr=68 - | когда приближенное значение решения для жесткой системы не может быть вычислено с заданной точностью; для достижения тебуемой точности следует воспользоваться версиями подпрограммы de23d_c, de25r_c или de25d_c. |
Версии
de23d_c - | вычисление решения задачи Коши для нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с повышенной точностью. При этом параметры xn, yn, xk, hmin, hmax, eps, h, y, ypm, delty, rab, yp и параметры x, y и dy в подпрограмме f должны иметь тип double. |
de25r_c - |
вычисление решения задачи Коши для жесткой
системы обыкновенных дифференциальных уравнений
первого порядка в конце интервала
интегрирования методом Гира с автоматическим выбором шага.
Первый оператор подпрограммы имеет вид: int de25r_c(real *f, real *fj, integer *m, real *xn, real *yn, real *xk, real *hmin, real *hmax, real *eps, integer *iorder, integer *iu, real *h, real *y, real *ypm, real *delty, real *rab, real *yp, integer *ierr) Здесь: fj - имя подпрограммы вычисления якобиана правой части системы; первый оператор этой подпрограммы имеет вид: int fj (float *x, float *y, float *z, int *m) Здесь: x, y - значения независимой и зависимой переменных, соответственно, причем y представляет одномерный массив длины m; вычисленное значение якобиана должно быть помещено в двумерный массив z размера m*m, при этом частная производная от правой части i - ого уравнения по j - ой переменной y (j) запоминается в элементе z (i, j) (тип параметров x, y и z: вещественный). Остальные параметры подпрограммы de25r_c имеют тот же смысл, что и одноименные параметры подпрограммы de23r_c. |
de25d_c - | вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de25r_c; при этом параметры xn, yn, xk, hmin, hmax, eps, h, y, ypm, delty, rab, yp и параметры x, y, dy и z в подпрограммах f и fj должны иметь тип double. |
Вызываемые подпрограммы
de21r_c - | выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы de23r_c. |
de21d_c - | выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью; вызывается при работе подпрограммы de23d_c. |
de24r_c - | выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы de25r_c. |
de24d_c - | выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью; вызывается при работе подпрограммы de25d_c. |
utde12_c - | подпрограмма выдачи диагностических сообщений при работе подпрограмм de21r_c, de23r_c, de24r_c, de25r_c. |
utde13_c - | подпрограмма выдачи диагностических сообщений при работе подпрограмм de21d_c, de23d_c, de24d_c, de25d_c. |
Замечания по использованию
B общем случае заданая точность eps не гарантируется. При работе подпрограммы и ее версий значения параметров m, xn, yn, xk, hmin, hmax, eps, istifj, iorder, iu сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить. Значение hmin должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром h при обращении к подпрограмме и ее версиям, т.к. интегрирование системы начинается с метода первого порядка. |
y1' = -20 y1 + y2 y2' = 19 y1 - 2 y2 , 0 ≤ x ≤ 1 y1(0) = 2 , y2(0) = 18 Собственные значения якобиевой матрицы системы λ1 = -21 , λ2 = -1 . Поэтому эту систему можно рассматривать как жесткую. int main(void) { /* Local variables */ extern int de23r_c(U_fp, int *, float *, float *, float *, float *, float *, float *, int *, int *, int *, float *, float *, float *, float *, float *, float *, int *); static float hmin, hmax; static int ierr; extern int f_c(); static float h__; static int m; static float x, y[2], delty[2]; static int iu; static float xk, xn, yn[2], yp[16]; static int iorder, istifj; static float rab[38], eps, ypm[2]; iu = 3; eps = 1e-5f; iorder = 6; hmin = 1e-10f; m = 2; xn = 0.f; yn[0] = 2.f; yn[1] = 18.f; xk = 1.f; istifj = 1; hmax = .1f; h__ = .01f; de23r_c((U_fp)f_c, &m, &xn, yn, &xk, &hmin, &hmax, &eps, &istifj, &iorder, &iu, &h__, y, ypm, delty, rab, yp, &ierr); printf("\n %16.7e \n", x); printf("\n %16.7e %16.7e \n", y[0], y[1]); return 0; } /* main */ int f_c(float *x, float *y, float *dy, int *m) { /* Parameter adjustments */ --dy; --y; /* Function Body */ dy[1] = y[1] * -20.f + y[2]; dy[2] = y[1] * 19.f - y[2] * 2.f; return 0; } /* f_c */ Результаты: ierr = 0 y(1) = 3.678796482-01 y(2) = 6.989709436