Текст подпрограммы и версий de37r_p.zip , de37e_p.zip |
Тексты тестовых примеров tde37r_p.zip , tde37e_p.zip |
Вычисление решения задачи Коши для жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования неявным методом Рунге - Кутта.
Решается задача Коши для жесткой линейной системы M обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
Y ' (X) = A(X) * Y(X) + φ(X) , Y = ( y1, ..., yM ) , A(X) = ( ai j(X) ), i, j = 1, ..., M , φ(X) = ( φ1(X), ..., φM(X) ) .
с начальными условиями, заданными в точке XN:
Y (XN) = YN, YN = ( y10, ..., yM0 ) ,
трехстадийным А - устойчивым неявным методом Рунге - Кутта шестого порядка точности. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования.
Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
Современные численные методы решения обыкновенных дифференциальных уравнений. Ред. Дж.Холл и Дж.Уатт. "Мир", M., 1979.
Butcher J.C. Implicit Runge - Kutta processes. Math. Comp., 18, 50 - 64, 1964.
procedure DE37R(FA :Proc_F3_DE; FI :Proc_F3_DE; M :Integer; XN :Real; var YN :Array of Real; XK :Real; HMIN :Real; EPS :Real; P :Real; var H :Real; var Y :Array of Real; var RAB :Array of Real; var IR :Array of Integer; var IERR :Integer);
Параметры
FA - |
подпрограмма вычисления матрицы системы A (X) в
любой точке X. Первый оператор подпрограммы
должен иметь вид:
procedure FA (var A :Array of Real; X :Real; M :Integer); Здесь: A - двумерный массив размера M * M, в который помещается матрица системы, вычисленная при значении аргумента X (тип параметров A, X: вещественный); |
FI - |
подпрограмма вычисления неоднородности правой
части системы φ (X)
в любой точке X. Первый
оператор подпрограммы должен иметь вид:
procedure FI (var G :Array of Real; X :Real; M :Integer); Здесь G - одномерный массив длины M, в который помещается неоднородность правой части системы, вычисленная при значении аргумента X (тип параметров G, X: вещественный); |
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 - | одномерный рабочий массив вещественного типа длины (10*M*M + 9*M + 1); |
IR - | целый одномерный рабочий массив длины 3*M; |
IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая-нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. |
Версии
DE37E - | вычисление решения задачи Коши для жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования неявным методом Рунге - Кутта с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, RAB и параметры A, G, X в подпрограммах FA и FI должны иметь тип Extended. |
Вызываемые подпрограммы
DE36R - DE36E | выполнение одного шага численного интегрирования жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка неявным методом Рунге - Kутта. |
UTDE16 - UTDE17 | подпрограммы выдачи диагностических сообщений. |
Подпрограммы DE36R, UTDE16 вызываются при работе подпрограммы DE37R, а подпрограммы DE36E, UTDE17 - при pаботе DE37E. |
Замечания по использованию
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, P сохраняются. При работе подпрограмм FA и FI значения параметров X и M не должны изменяться. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65, значение параметра YN будет испорчено. Tак как при интегрировании уравнений с помощью подпрограмм DE37R и DE37E используются глобальные записи (структуры данных) с именами _COM36R и _COM36D, соответственно, то пользователю не рекомендуется использовать для своих целей глобальные записи (структуры данных) с указанными именами. |
y1' = - 20y1 + y2 , y1 = 2 , y2' = 19y1 - 2y2 , y2 = 18 , 0 ≤ x ≤ 5 . Точное решение системы: y1 = e - x + e - 21x , y2 = 19e - x - e - 21x .Unit TDE37R_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc, UtRes_p, FB1DE37R_p, FB2DE37R_p, DE37R_p; function TDE37R: String; implementation function TDE37R: String; var M,K,_i,IERR :Integer; XN,XK,HMIN,EPS,P,H,X,Y1,Y2 :Real; YN :Array [0..1] of Real; Y :Array [0..1] of Real; RАВ :Array [0..58] of Real; IR :Array [0..5] of Integer; label _200; begin Result := ''; { результат функции } M := 2; XN := 0.0; YN[0] := 2.0; YN[1] := 18.0; ХК := 5.0; HMIN := 1.E-5; EPS := 1.E-5; P := 100.0; for K:=1 to 2 do begin H := 0.01; DE37R(FB1DE37R,FB2DE37R,M,XN,YN,XK,HMIN,EPS,P,H,Y,RAB,IR,IERR); Result := Result + Format('%s',[' IERR=']); Result := Result + Format('%4d ',[IERR]) + #$0D#$0A; X := XK; Y1 := Exp(-X)+Exp(-21.0*X); Y2 := 19.0*Exp(-X)-Exp(-21.0*X); Result := Result + Format('%20.16f ',[XK]) + #$0D#$0A; Result := Result + #$0D#$0A; for _i:=0 to 1 do begin Result := Result + Format('%20.16f ',[Y[_i]]); if ( ((_i+1) mod 2)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + Format('%20.16f %20.16f %20.16f ', [H,Y1,Y2]) + #$0D#$0A; EPS := 0.5E-5; _200: end; UtRes('TDE37R',Result); { вывод результатов в файл TDE37R.res } exit; end; end. Unit fb1de37r_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc; procedure fb1de37r(var A :Array of Real; X :Real; M :Integer); implementation procedure fb1de37r(var A :Array of Real; X :Real; M :Integer); begin A[0] := -20.0; A[2] := 1.0; A[1] := 19.0; A[3] := -2.0; end; end. Unit fb2de37r_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc; procedure fb2de37r(var FI :Array of Real; X :Real; M :Integer); implementation procedure fb2de37r(var FI :Array of Real; X :Real; M :Integer); begin FI[0] := 0.0; FI[1] := 0.0; end; end. Результаты: после первого обращения к подпрограмме - XK Y(1) Y(2) 5.000000000000 + 00 6.737934370626 - 03 1.280207530417 - 01 H Y1 Y2 1.280000000001 + 00 6.737946999117 - 03 1.280209929830 - 01 после второго обращения к подпрограмме - XK Y(1) Y(2) 5.000000000000 + 00 6.737946721380 - 03 1.280209877059 - 01 H Y1 Y2 6.400000000003 - 01 6.737946999117 - 03 1.280209929830 - 01