Текст подпрограммы и версий de32r_p.zip , de32e_p.zip |
Тексты тестовых примеров tde32r_p.zip , tde32e_p.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 ) ,
с помощью модифицированного приема А.Н.Крылова отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Адамса значения решения Y в нескольких первых узлах и конечные разности правой части F(X,Y) системы до некоторого порядка в точке XN. Количество таких узлов и максимальный порядок получаемых разностей определяются в подпрограмме автоматически и соответствуют порядку точности используемого метода Адамса, который задается пользователем при обращении к подпрограмме.
Указанная процедура вычисления необходимых для счета начальных значений, или "фронта Адамса", называется разгоном. Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на экстраполяционную формулу Адамса, записанную в разностной форме.
B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы, при этом все компоненты решения проверяются на точность по меpе погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, называемой границей перехода, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып.3., Изд-во МГУ, 1965.
procedure DE32R(F :Proc_F_DE; M :Integer; IORDER :Integer; XN :Real; var YN :Array of Real; HMIN :Real; EPS :Real; P :Real; var H :Real; var DF :Array of Real; var X :Array of Real; var Y :Array of Real; var BUL :Boolean; var RS :Array of Real; var RF :Array of Real; var R :Array of Real; var IERR :Integer);
Параметры
F - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператоp подпрограммы должен иметь вид: procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer); Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный); |
M - | количество уравнений в системе (тип: целый); |
IORDER - | порядок точности того метода Адамса, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений; IORDER должен быть не больше 10 (тип: целый); |
XN, YN - | начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный); |
HMIN - | минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
EPS - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
P - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
H - | вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность при разгоне достигается, то именно он и реализуется на разгоне, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность; |
DF - | двумерный вещественный массив размера M*IORDER, в котоpом запоминаются значения правой части системы и ее разностей до порядка (IORDER - 1) включительно, вычисленные в точке XN, при этом элемент DF (I, 1) этого массива содержит значение правой части I - ого уравнения, а DF (I, J + 1) - ее J - ю разность, погрешность вычисления которой имеет (IОRDЕR + 1) - й порядок по H; |
X - | одномерный вещественный массив длины IORDER, содержащий на выходе из подпрограммы IORDER узлов, включая XN, в которых вычисляются при разгоне приближенные значения решения, при этом элемент X (J) этого массива содержит узел, равный XN + (J - 1) * H; |
Y - | двумерный вещественный массив размера M*IORDER, в котоpом запоминаются приближенные значения решения, вычисленные для значений аргумента, хранящихся в массиве X, а именно, значению аргумента в X (J) соответствует приближенное решение Y (I,J), при этом погрешность этого решения имеет порядок IORDER по H; |
BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE, если заданный в H шаг выводит в конец интервала интегрирования,т.е. узел X (IORDER) совпадает с концом интервала интегрирования, и FALSE в противном случае. B результате работы подпрограммы BUL pавно FALSE, если вместо исходного шага интегрирования при разгоне был использован меньший шаг; в противном случае значение переменной BUL не меняется; |
RS - RF, R | одномерные вещественные рабочие массивы длины M; |
IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
IERR= 1 - | когда неправильно задан параметр IORDER,т.е. IORDER > 10; в этом случае разгон выполняется для значения IORDER = 10; |
IERR=65 - | когда на разгоне не может быть достигнута требуемая точность EPS; в этом случае разгон можно начать сначала обращение к подпрограмме с новыми значениями параметров H, HMIN и IORDER. |
Версии
DE32E - | построение начальных значений с расширенной (Extended) точностью при интегрировании системы обыкновенных дифференциальных уравнений первого порядка методом Адамса с контролем точности. При этом параметры XN, YN, HMIN, EPS, P, H, DF, X, Y, RS, RF, R и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. |
Вызываемые подпрограммы
UTDE12 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE32R. |
UTDE13 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE32E. |
Kpоме того, подпрограммы DE32R и DE32E используют рабочие подпрограммы DE28RS и DE28ES, соответственно. |
Замечания по использованию
B общем случае требуемая точность при разгоне не гарантируется. При работе подпрограммы и ее версии значения параметров M, IORDER, XN, YN, HMIN, EPS, P сохраняются. Подпрограмма и ее версия могут использоваться как для разгона, так и непосредственно для численного решения системы уравнений на всем интервале интегрирования. |
y1' = y2 y2' = -y1 y1 (3/4 π) = √2 /2 y2 (3/4 π) = -√2 /2Unit TDE32R_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc, UtRes_p, FDE32R_p, DE32R_p; function TDE32R: String; implementation function TDE32R: String; var IORDER,M,J,I,IERR :Integer; P,EPS,HMIN,H,XN :Real; BUL :Вoolean; YN :Array [0..1] of Real; DF :Array [0..9] of Real; Y :Array [0..9] of Real; X :Array [0..4] of Real; RS :Array [0..1] of Real; RF :Array [0..1] of Real; R :Array [0..1] of Real; begin Result := ''; { результат функции } BUL := False; IORDER := 5; P := 100.0; EPS := 0.0001; HMIN := 1.E-10; M := 2; Н := 0.01; XN := 0.75*3.14159265359; YN[0] := Sqrt(2.0)/2.0; YN[1] := -YN[0]; DE32R(FDE32R,M,IORDER,XN,YN,HMIN,EPS,P,H,DF,X,Y,BUL,RS,RF,R,IERR); Result := Result + Format(' %20.16f %20.16f %20.16f ', [X[4],Y[8],Y[9]]) + #$0D#$0A; Result := Result + #$0D#$0A; for J:=1 to 5 do begin for I:=1 to 2 do begin Result := Result + Format('%20.16f ',[DF[(I-1)+(J-1)*2]]) + #$0D#$0A; end; end; Result := Result + #$0D#$0A; UtRes('TDE32R',Result); { вывод результатов в файл TDE32R.res } exit; end; end. Unit fde32r_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc; procedure fde32r(X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer); implementation procedure fde32r(X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer); begin DY[0] := Y[1]; DY[1] := -Y[0]; end; end. Результаты: IERR = 0 X(5) = 2.396194490 + 00 Y(1, 5) = 6.782644418-01 Y(2, 5) = -7.348179005-01 DF(1, 1) = -7.07106781186-03 DF(1, 2) = -7.10630480683-05 DF(1, 3) = 6.99988802921-07 DF(1, 4) = 7.18290493751-09 DF(1, 5) = -7.27808924239-11 DF(2, 1) = -7.07106781186-03 DF(2, 2) = 7.03559472512-05 DF(2, 3) = 7.14142103675-07 DF(2, 4) = -6.97064450605-09 DF(2, 5) = -6.86739554112-11