Текст подпрограммы и версий de00r_c.zip , de00d_c.zip , de84r_c.zip , de84d_c.zip |
Тексты тестовых примеров tde00r_c.zip , tde00d_c.zip , tde84r_c.zip , tde84d_c.zip |
Выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = F(X, Y) , Y = ( y1,..., yM ) , F = ( f1( X, y1,..., yM ), ... , fM( X, y1,..., yM ) )
Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Предполагается также, что система (1) на [ xn, xn + H ] может быть достаточно хорошо аппроксимирована линейной однородной системой с постоянными коэффициентами, т.е. системой вида
Y ' = AY ,
где A - некоторая постоянная матрица. Вычисление производится по методу Лоусона. Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции y (x) на [ xn, xn + H ] по формуле
y(x) = exp [ ( x - xn ) A0 ] Z(x) ,
где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X) :
(2) Z ' (x) = F1(x, z) = exp [ - ( x - xn ) A0 ] { F(x, exp [ ( x - xn ) A0 ] Z(x)) - - A0 exp [ ( x - xn ) A0 ] Z(x) }
Матрица Якоби ∂ F1 / ∂Z системы (2) и матрица Якоби ∂F / ∂y системы (1) связаны между собой соотношением
∂ F1 / ∂Z = exp [ - (x - xn) A0] { ∂F / ∂y - A0 } exp [ ( x - xn ) A0 ]
Указанное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу ∂F / ∂y , вычисленную либо в точке ( xn, yn ), либо в некоторой предыдущей точке ( xn - k, yn - k ). Если шаг H достаточно мал, то данное преобразование позволяет уменьшить характеристические корни матрицы ∂ F1 / ∂Z по сравнению с характеристическими корнями матрицы ∂F / ∂y . А это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1).
Для решения системы (2) применяются классический метод Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).
Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
J.Douglas Lowson, Generalized Runge - Kutta processes for stable systems with large lipshitz constants, SIAM Journal on Numerical Analisys, 1967, Vol.4, No 3.
int de00r_c (real *f, S_fp fj, integer *m, integer *jstart, real *hmin, real *eps, real *p, real *yx, real *x, real *h, logical *bul, real *xp, real *yp, real *a, real *e1, real *e2, real *e3, real *e4, real *r, real *r0, real *r1, real *r2, real *r3, real *yd, integer *ierr)
Параметры
f - |
подпрограмма вычисления правой части системы в
точке x . Первый оператор подпрограммы должен
иметь вид: int f (float *x, float *y, float *dy, int *m). Здесь x и y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. m ≠ 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y и dy: вещественный); |
fj - |
подпрограмма вычисления матрицы Якоби правой
части системы. Первый оператор подпрограммы должен иметь вид: int fj (float *x, float *y, float *df, int *m). Здесь x и y - значения независимой и зависимой переменных, соответственно, причем y представляет одномерный массив длины m. df - двумерный массив размера m*m. Вычисленное значение матрицы Якоби ∂F / ∂y должно быть помещено в массив df, при этом частная производная от правой части i - ого уравнения F по j - ой переменной y (j) запоминается в dF (i, j), т.е. dF (i, j) = ∂Fi / ∂yj (тип параметров x, y и df: вещественный); |
m - | количество уравнений в системе (тип: целый); |
jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart; при этом значение параметра в качестве матрицы A0 принимается матрица ∂F / ∂y в точке [ x = xn, y = yn ]; вычисленное значение матрицы A0 будет использоваться в дальнейшем для значений jstart = 1, 2, 3; |
1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, заданных параметрами yx и x соответственно; при этом само значение величины шага интегрирования H, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней, а матрица A0 остается той же, что и на предыдущем шаге; |
2 - | то же, что и для jstart = 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме; матрица A0 остается той же, что и на предыдущем шаге; |
3 - | то же, что и для jstart = 1, но только с той разницей, что шаг интегрирования может принимать произвольное значение, матрица A0 остается той же, что и на предыдущем шаге. |
-1 - | повторить последний шаг интегрирования с новыми значениями параметров h и/или hmin; при этом значении jstart в качестве матрицы A0 принимается матрица ∂F / ∂y в точке [ x = xp, y = yp ] таким образом, до тех пор, пока jstart принимает значения 1, 2 или 3, сохраняется одно и то же значение матрицы A0 , равное значению матрицы Якоби, вычисленному при последнем обращении к подпрограмме с параметром jstart = 0 или jstart = - 1; если требуется изменить матрицу A0 , то следует обратиться к подпрограмме со значением jstart = 0 или - 1; при выходе из подпрограммы параметр jstart полагается равным 2, если выполнив заданный в h шаг, подпрограмма рекомендовала для интегрирования вдвое большеe его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной h; |
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 при входе в нее (т.е. предыдущий узел и решение в нем); |
a, e1, e2 - e3, e4 | вещественные двумерные рабочие массивы размера m*m; |
r, r0 - r1, r2 r3, yd | вещественные одномерные рабочие массивы длины m; |
ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1. |
Версии
de00d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме f и x, y, df в подпрограмме fj должны иметь тип double; |
de84r_c - |
выполнение одного шага численного интегрирования
устойчивой системы обыкновенных дифференциальных
уравнений первого порядка с большой константой Липшица
методом Лоусона. В отличии от подпрограммы de00r_c
первый оператор подпрограммы de84r_c не содержит в
списке формальных параметров имени подпрограммы
вычисления матрицы Якоби и имеет следующий вид: int de84r_c (real *f, integer *m, integer *jstart, real *hmin, real *eps, real *p, real *yx, real *x, real *h, logical *bul, real *xp, real *yp, real *a, real *e1, real *e2, real *e3, real *e4, real *r, real *r0, real *r1, real *r2, real *r3, real *yd, integer *ierr). При использовании этой подпрограммы пользователь не составляет подпрограмму вычисления матрицы Якоби. Все параметры подпрограммы de84r_c имеют тот же смысл, что и одноименные параметры подпрограммы de00r_c; |
de84d_c - | выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de84r_c, при этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
ame2r_c - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de00r_c; |
ame2d_c - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de00d_c; |
utde20_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de00r_c; |
utde21_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de00d_c; |
Кроме того, при работе подпрограмм de00r_c, de84r_c и de00d_c, de84d_c вызываются рабочие подпрограммы de00rs_c, de04rm_c, de00rd_c и de00ds_c, de04dm_c, de00dd_c соответственно. | |
При работе подпрограмм de84r_c и de84d_c дополнительно вызываются также рабочие подпрограммы de84rj_c и de84dj_c соответственно. |
Замечания по использованию
Данная подпрограмма предназначена для интегрирования таких
нелинейных систем, которые могут быть достаточно хорошо
аппроксимированы линейной однородной системой с постоянными
коэффициентами, т.е. системой вида Y ' = AY , где A - некоторая постоянная матрица. В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров m, hmin, eps, p сохраняются. Между последовательными обращениями к подпрограмме со значениями параметра jstart = 1, 2, 3 пользователь не должен изменять содержимое массивов a, e1, e2. При работе подпрограмм f и fj значения параметров m, x, y не должны изменяться. При обращении к подпрограмме со значением jstart = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров xp и yp соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart. После работы подпрограммы в массиве r1 содержится значение оценки абсолютной погрешности на шаге, вычисленной по правилу Рунге. Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (1). |
y1' = (- 500 + 2 x) y1 + (1 - x / 10 )y2 + (1 + 1 / ( x+1 ) ) sin y1 , y1( = 49 y2' = (- 1000 - 10 x) y1 + (1 + x / 10 ) y2 + (2 + 1 / ( x+2 ) ) cos y2 , y2( = 50 1 ≤ x ≤ 6
Приводятся подпрограммы вычисления правой части системы и матрицы Якоби правой части, а также фрагмент вызывающей программы и результаты счета.
int main(void) { /* Local variables */ extern int de00r_c(U_fp, U_fp, int *, int *, float *, float *, float *, float *, float *, float *, logical *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, int *); static float hmin; static int ierr; static float a[4] /* was [2][2] */; extern int f_c(); static float h__; static int m; static float p, r__[2], x, e1[4] /* was [2][2] */, e2[4] /* was [2][2] */, e3[4] /* was [2][2] */, e4[4] /* was [2][2] */, r0[2], r1[2], r2[2], r3[2]; extern int fj_c(); static int ih; static float yd[2], xp, yp[2], yx[2]; static int jstart; static logical bul; static float eps; m = 2; x = 1.f; yx[0] = 49.f; yx[1] = 50.f; hmin = 1e-10f; eps = 1e-8f; p = 100.f; jstart = 0; h__ = .01f; ih = 0; l100: ++ih; de00r_c((U_fp)f_c, (U_fp)fj_c, &m, &jstart, &hmin, &eps, &p, yx, &x, &h__, &bul, &xp, yp, a, e1, e2, e3, e4, r__, r0, r1, r2, r3, yd, &ierr); printf("\n %16.7e %16.7e %16.7e \n", x, yx[0], yx[1]); printf("\n %16.7e %16.7e %16.7e \n\n", h__, r1[0], r1[1]); switch (ih) { case 1: goto l101; case 2: goto l102; case 3: goto l103; case 4: goto l104; case 5: goto l105; case 6: goto l106; } l101: h__ = 1.f; jstart = 3; goto l100; l102: eps = .01f; goto l100; l103: jstart = -1; eps = 1e-8f; goto l100; l104: eps = .01f; h__ = 1e-8f; jstart = 3; goto l100; l105: h__ = -.001f; eps = 1e-8f; jstart = 0; goto l100; l106: return 0; } /* main */ int f_c(float *x, float *y, float *dy, int *m) { /* Builtin functions */ double sin(double), cos(double); /* Parameter adjustments */ --dy; --y; /* Function Body */ dy[1] = (*x * 2.f - 500.f) * y[1] + (1.f - *x * .1f) * y[2] + (1.f / (*x + 1.f) + 1.f) * (float)sin(y[1]); dy[2] = (-1e3f - *x * 10.f) * y[1] + (*x * .1f + 1.f) * y[2] + (1.f / (*x + 2.f) + 2.f) * (float)cos(y[2]); return 0; } /* f_c */ int fj_c(float *x, float *y, float *df, int *m) { /* Builtin functions */ double cos(double), sin(double); #define df_ref(a_1,a_2) df[(a_2)*2 + a_1] /* Parameter adjustments */ df -= 3; --y; /* Function Body */ df_ref(1, 1) = *x * 2.f - 500.f + (1.f / (*x + 1.f) + 1.f) * (float)cos(y[1]); df_ref(1, 2) = 1.f - *x * .1f; df_ref(2, 1) = -1e3f - *x * 10.f; df_ref(2, 2) = *x * .1f + 1.f - (1.f / (*x + 2.f) + 2.f) * (float)sin(y[2]); return 0; } /* fj_c */ Результаты: После первого обращения к подпрограмме x yx(1) yx(2) 1.000039062499+00 4.895766071531+00 4.808757932455+00 h__ r1(1) r1(2) 3.906250000002-05 9.597251282538-10 -7.233271996176-09 После второго обращения к подпрограмме x yx(1) yx(2) 1.000069580077+00 4.733408100787+01 4.661901490984+01 h__ r1(1) r1(2) 3.051757812500-05 1.319373647371-10 6.946114202335-02 После третьего обращения к подпрограмме x yx(1) yx(2) 1.000100097655+00 4.662140584405+01 4.517257692956+01 h__ r1(1) r1(2) 6.103515625000-05 -6.208817164108-11 1.830362084864-09 После четвертого обращения к подпрограмме x yx(1) yx(2) 1.000100097655+00 4.662140584411+01 4.517257692973+01 h__ r1(1) r1(2) 3.051757812500-05 -5.820766091346-11 1.851003617047-09 После пятого обращения к подпрограмме x yx(1) yx(2) 1.000100107654+00 4.662117408356+01 4.517210655950+01 h__ r1(1) r1(2) 2.000000000004-08 0.000000000000 + 00 0.000000000000 + 00 После шестого обращения к подпрограмме x yx(1) yx(2) 1.000068857653+00 4.735108317091+01 4.665352269082+01 h__ r1(1) r1(2) -3.125000000004-05 3.880510727560-11 -2.250696221985-09