|
Текст подпрограммы и версий de94r_p.zip, de94e_p.zip, de96r_p.zip, de96e_p.zip, de90r_p.zip, de90e_p.zip, de92r_p.zip, de92e_p.zip |
Тексты тестовых примеров tde94r_p.zip, tde94e_p.zip, tde96r_p.zip, tde96e_p.zip , tde90r_p.zip, tde90e_p.zip, tde92r_p.zip, tde92e_p.zip |
Выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
Y ' = F( X, Y ) ,
Y = ( y1, ..., yM ) , F = ( f1 ( X, y1, ..., yM ), ..., fM ( X, y1, ..., yM ) )
методом типа Розенброка четвертого порядка.
По заданному значению YX решения в узле Xn вычисляется значение этого решения в узле Xn + H. Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Артемьев С.С., Демидов Г.В. A - устойчивый метод типа Розенброка четвертого порядка точности решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Некоторые проблемы вычислительной и прикладной математики. Новосибирск: "Наука", 1975.
Современные численные методы решения обыкновенных дифференциальных уравнений. Под ред. Дж.Холла и Дж.Уатта. М.: "Мир", 1979.
procedure DE94R(F :Proc_F_DE; FJ :Proc_F_DE; FX :Proc_F_DE; M :Integer;
var JSTART :Integer; HMIN :Real; EPS :Real; P :Real;
var YX :Array of Real; var X :Real; var H :Real;
var BUL :Boolean; var XP :Real; var YP :Array of Real;
var A :Array of Real; var A1 :Array of Real;
var A2 :Array of Real; var RAB :Array of Real;
var R1 :Array of Real; var R2 :Array of Real;
var R3 :Array of Real; var R4 :Array of Real;
var YD :Array of Real; var IERR :Integer);
Параметры
| F - |
имя подпрограммы вычисления значений правой
части дифференциальных уравнений. Первый оператор
подпрограммы должен иметь вид: procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer); Здесь X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY; в случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный); |
| FJ - |
имя подпрограммы вычисления значений элементов матрицы Якоби
∂f / ∂y
правой части системы. Первый оператор подпрограммы должен иметь вид: procedure FJ (X :Real; var Y :Array of Real; var Z :Array of Real; M :Integer); Здесь X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметр Y представляет собой одномерный массив длины M, а параметр Z - двумерный массив размера M * M. Значения элементов матрицы Якоби ∂f / ∂y правой части системы должны быть помещены в массив Z, при этом частная производная от правой части I - го уравнения по J - ой переменной Y ( J ) запоминается в элементе Z ( I, J )(тип параметров X, Y и Z: вещественный); |
| FX - |
имя подпрограммы вычисления значений частных
производных по X правой части системы
∂f / ∂x.
Первый оператор подпрограммы должен иметь вид: procedure FX (X :Real; var Y :Array of Real; var Z1 :Array of Real; M :Integer); Здесь: X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметры Y и Z1 представляют собой одномерные массивы длиной M. Значения частных производных по X правой части системы ∂f / ∂x должны быть помещены в массив Z1; при этом частная производная от правой части I - го уравнения запоминается в элементе Z1 ( I ) (тип параметров X, Y и Z1: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN; |
| на выходе из подпрограммы JSTART равен + 1; | |
| 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 - A1, A2 | вещественные двумерные рабочие массивы размера M * M; |
| RAB - | вещественный одномерный рабочий массив длины 5 * M; на выходе из подпрограммы в первых M элементах этого массива содержится погрешность вычисленного приближенного решения, полученная по правилу Рунге; |
|
R1, R2 - R3, R4 YD | вещественные одномерные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE96R - |
выполнение одного шага численного интегрирования
жесткой системы обыкновенных дифференциальных уравнений
первого порядка методом типа Розенброка четвертого
порядка. Отличие подпрограммы DE96R от подпрограммы
DE94R состоит в том, что при обращении к подпрограмме
DE96R не требуется задавать подпрограмму вычисления
матрицы Якоби системы и подпрограмму вычисления
частных производных правой части системы по x . Все
частные производные от правой части вычисляются в
подпрограмме DE96R с помощью разностных отношений.
Первый оператор подпрограммы DE96R имеет вид:
procedure DE96R(F :Proc_F_DE; M :Integer; var JSTART :Integer;
HMIN :Real; EPS :Real; P :Real; var YX :Array of Real;
var X :Real; var H :Real; var BUL :Boolean;
var XP :Real; var YP :Array of Real;
var A :Array of Real; var A1 :Array of Real;
var A2 :Array of Real; var RAB :Array of Real;
var R1 :Array of Real; var R2 :Array of Real;
var R3 :Array of Real; var R4 :Array of Real;
var YD :Array of Real; var IERR :Integer);
Список формальных параметров подпрограммы DE96R
отличается от списка параметров подпрограммы DE94R
отсутствием параметров FJ и FX.
Параметры подпрограммы DE96R имеют тот же смысл, что
и одноименные параметры подпрограммы DE94R.
|
| DE90R - |
выполнение одного шага численного интегрирования
автономной жесткой системы обыкновенных дифференциальных
уравнений первого порядка методом типа Розенброка четвертого порядка.
Данная подпрограмма предназначена для интегрирования
систем уравнений, в правые части которых не входит
независимая переменная x ,
.
Первый оператор подпрограммы DE90R имеет вид:
procedure DE90R(F :Proc_F_DE; FJ :Proc_F_DE; M :Integer;
var JSTART :Integer; HMIN :Real; EPS :Real; P :Real;
var YX :Array of Real; var X :Real; var H :Real;
var BUL :Boolean; var XP :Real; var YP :Array of Real;
var A :Array of Real; var A1 :Array of Real;
var A2 :Array of Real; var RAB :Array of Real;
var R1 :Array of Real; var R2 :Array of Real;
var YD :Array of Real; var IERR :Integer);
Список формальных параметров подпрограммы DE90R
отличается от списка параметров DE94R отсутствием
параметров FX, R3, R4.
Параметры подпрограммы DE90R имеют тот же смысл, что
и одноименные параметры подпрограммы DE94R, кроме
параметра RAB. В подпрограмме DE90R параметр RAB
представляет одномерный вещественный рабочий массив длины
4 * M.
|
| DE92R - |
выполнение одного шага численного интегрирования
жесткой автономной системы обыкновенных дифференциальных
уравнений первого порядка методом типа Розенброка
четвертого порядка. Данная подпрограмма предназначена
для интегрирования систем уравнений, в правые части
которых не входит независимая переменная X, т.е.
.
Отличие подпрограммы DE92R от
подпрограммы DE94R состоит в том, что при обращении
к подпрограмме DE92R не требуется задавать подпрограмму
вычисления матрицы Якоби системы
¶f / ¶y
и подпрограмму вычисления частных производных
¶f / ¶x .
Первый оператор подпрограммы DE92R имеет вид:
procedure DE92R(F :Proc_F_DE; M :Integer; var JSTART :Integer;
HMIN :Real; EPS :Real; P :Real; var YX :Array of Real;
var X :Real; var H :Real; var BUL :Boolean;
var XP :Real; var YP :Array of Real;
var A :Array of Real; var A1 :Array of Real;
var A2 :Array of Real; var RAB :Array of Real;
var R1 :Array of Real; var R2 :Array of Real;
var YD :Array of Real; var IERR :Integer);
Список формальных параметров подпрограммы DE92R
отличается от списка параметров подпрограммы DE94R
отсутствием параметров FJ, FX, R3, R4.
Параметры подпрограммы имеют тот же смысл, что и в
подпрограмме DE94R, кроме параметра RAB. В подпрограмме
DE92R параметр RAB представляет одномерный вещественный
рабочий массив длины 4 * M.
|
| DE94E - | выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE94R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD и параметры X, Y, DY, Z, Z1 в подпрограммах F, FJ, FX должны иметь тип Extended; |
| DE96E - | выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE96R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD и параметры X, Y, DY в подпрограмме F должны иметь тип Extended. |
| DE90E - | выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE90R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, YD и параметры X, Y, DY, Z в подпрограммах F, FJ должны иметь тип Extended. |
| DE92E - | выполнение одного шага численного интегрирования жесткой автономной системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE92R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, YD и параметры X, Y, DY в подпрограмме F должны иметь тип Extended. |
Вызываемые подпрограммы
| AIG2R - | обращение вещественной матрицы методом Жордана с выбором ведущего элемента по всей матрице. Вызывается при работе подпрограмм DE90R, DE92R, DE94R, DE96R. |
| AIG2E - | обращение вещественной матрицы, заданной с расширенной (Extended) точностью, методом Жордана с выбором ведущего элемента по всей матрице. Вызывается при работе подпрограмм DE90E, DE92E, DE94E, DE96E. |
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограмм DE90R, DE92R, DE94R, DE96R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограмм DE90E, DE92E, DE94E, DE96E. |
| Кроме того, при работе подпрограммы DE90R вызываются рабочие подпрограммы DE92RS, DE92RK; при работе подпрограммы DE92R вызываются рабочие подпрограммы DE92RS, DE92RK, DE84RJ; при работе подпрограммы DE94R выэываются рабочие подпрограммы DE96RS, DE92RK; при работе подпрограммы DE96R вызываются рабочие подпрограммы DE96RS, DE92RK, DE84RJ, DE96RX; при работе подпрограммы DE90E вызываются рабочие подпрограммы DE92ES, DE92EK; при работе подпрограммы DE92E вызываются подпрограммы DE92ES, DE92EK, DE84EJ; при работе подпрограммы DE94E вызываются подпрограммы DE94ES, DE92EK; при работе подпрограммы DE96D вызываются подпрограммы DE96ES, DE92EK, DE84EJ, DE96EX. |
Замечания по использованию
|
В общем случае заданная точность не гарантируется. При работе подпрограммы и ее версий значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограмм F, FJ, FX значения параметров M, X, Y не должны изменяться. При выполнении из точки xn одного шага H метода типа Розенброка в подпрограмме DE94R и ее версиях вычисляются четыре значения правой части системы дифференциальных уравнений. Среди этих значений есть одно, которое должно вычисляться при значении независимой переменной x, равном xn - H. Это значение аргумента xn находится слева от точки xn при H > 0 и справа от точки xn при H < 0. Это следует учитывать при составлении подпрограммы F вычисления правой части системы. Если правая часть системы дифференциальных уравнений не определена для x < xn или для x < xn , то попытка вычислить правую часть для указанного значения аргумента x может привести к авосту.Если при решении системы дифференциальных уравнений не задаются подпрограммы FJ вычисления матрицы Якоби и FX вычисления частных производных по x правой части системы (т.е. производится , то все частные производные по y и по x в этих подпрограммах аппроксимируются центральными разностными отношениями. Замена точных значений производных разностными аппроксимациями может привести к росту погрешности приближенного решения системы дифференциальных уравнений по сравнению с тем, что было бы, если бы все частные производные вычислялись точно с помощью подпрограмм FJ и FX. Даже если будет мала погрешность аппроксимации частных производных (например, если правая часть системы является линейной функцией своих аргументов, то погрешность аппроксимации , может оказаться значительной вычислительная погрешность, особенно, когда вычисления выполняются с одинарной точностью. При этом вычислительная погрешность может даже превосходить верхний предел погрешности приближенного решения, заданный при обращении к подпрограмме параметром EPS. В этом случае целесообразно использовать не подпрограммы DE92R, DE96R, а их версии, выполняющие вычисления с удвоенным числом значащих цифр, т.е. подпрограммы DE92E, DE96E. Если частные производные по x от правой части системы вычисляются с помощью центральных разностных отношений (т.е. , то значения правой части системы, необходимые для аппроксимации этой производной, будут вычисляться с помощью подпрограммы F при значениях аргумента X, находящихся по разные стороны от точки Xn, а именно, при Xn - D, Xn + D где D - значение, задаваемое параметром EPS. Если правая часть системы дифференциальных уравнений не определена для x < xn при интегрировании слева направо или для x > xn при интегрировании справа налево, то попытка вычислить правую часть для указанного значения аргумента x может привести к авосту. |
Решается задача Коши
y1' = - 100y1
y2' = - 100y1 - 2y2 + 20e - 100x + 2e - x cos x
y3' = - 100y1 + 9998y2 - 9990y3 - 10y4 + 20e - 100x + 2e - x cos x
y4' = - 100y1 + 9988y2 + 20y3 - 10010y4 + 20e - 100x + 2e - x cos x
x = 0 , y1(0) = 10 , y2(0) = 11 , y3(0) = 111 , y4(0) = 111 .
Точное решение задачи имеет вид:
y1 = 10e - 100x
y2 = 10e - 100x + e - x cos x + e - x sin x
y3 = y2 + 100e - 10000x cos 10x
y4 = y3 + 100e - 10000x sin 10x
Матрица Якоби правой части системы имеет вид:
- 100 0 0 0
- 100 - 2 0 0
- 100 9998 - 9990 - 10
- 100 9988 20 - 10010
Частные производные по x правой части системы имеют вид:
0
- 2000e - 100x - 2e - x (cos x + sin x)
- 2000e - 100x - 2e - x (cos x + sin x)
- 2000e - 100x - 2e - x (cos x + sin x)
Ниже приводятся фрагмент вызывающей программы для DE94R, подпрограммы F, FJ, FX и результаты счета, полученные после нескольких обращений к подпрограмме DE94R.
Unit TDE94R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE94R_p, FJDE94R_p, FXDE94R_p, DE94R_p;
function TDE94R: String;
implementation
function TDE94R: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,Y1,Y2,Y3,Y4,ХР :Real;
BUL :Boolean;
YX :Array [0..3] of Real;
YP :Array [0..3] of Real;
A :Array [0..15] of Real;
A1 :Array [0..15] of Real;
A2 :Array [0..15] of Real;
RАВ :Array [0..19] of Real;
R1 :Array [0..3] of Real;
R2 :Array [0..3] of Real;
R3 :Array [0..3] of Real;
R4 :Array [0..3] of Real;
YD :Array [0..3] of Real;
label
_100,_101,_102,_103,_104,_105,_106;
begin
Result := ''; { результат функции }
M := 4;
X := 0.0;
YX[0] := 10.0;
YX[1] := 11.0;
YX[2] := 111.0;
YX[3] := 111.0;
HMIN := 1.E-10;
EPS := 1.E-4;
P := 1.E3;
JSTART := 0;
H := 0.01;
IH := 0;
_100:
IH := IH + 1;
DE94R (FDE94R,FJDE94R,FXDE94R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,
A,A1,A2,RAB,R1,R2,R3,R4,YD,IERR);
{ ВЫЧИСЛЕНИЕ ТОЧНОГО РЕШЕНИЯ CИCTEMЫ: }
Y1 := 10.0*Exp(-100.0*X);
Y2 := Y1+Exp(-X)*(Cos(X)+Sin(X));
Y3 := Y2+100.0*Exp(-1.E4*X)*Cos(10.0*X);
Y4 := Y3+100.0*Exp(-1.E4*X)*Sin(10.0*X);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y1,Y2]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y3,Y4]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[YX[0],YX[1]]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[YX[2],YX[3],X,H]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',
[' РАСПЕЧАТКА ПОГРЕШНОСТИ ПРИБЛИЖЕННОГО PEШEHИЯ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
{ РАСПЕЧАТКА ПОГРЕШНОСТИ ПРИБЛИЖЕННОГО РЕШЕНИЯ }
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[RAB[0],RAB[1],RAB[2],RAB[3]]) + #$0D#$0A;
case IH of
1: goto _101;
2: goto _102;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _106;
end;
_101:
H := 1.E-8;
goto _100;
_102:
H := -1.E-8;
goto _100;
_103:
JSTART := -1;
goto _100;
_104:
JSTART := -1;
H := -1.E-8;
goto _100;
_105:
H := 0.01;
goto _100;
_106:
UtRes('TDE94R',Result); { вывод результатов в файл TDE94R.res }
exit;
end;
end.
Unit fde94r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde94r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
implementation
procedure fde94r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
var
T1 :Real;
begin
DY[0] := -100.0*Y[0];
T1 := 20.0*Exp(-100.0*X)+2.0*Exp(-X)*Cos(X);
DY[1] := DY[0]-2.0*Y[1]+T1;
DY[2] := DY[0]+9998.0*Y[1]-9990.0*Y[2]-10.0*Y[3]+T1;
DY[3] := DY[0]+9988.0*Y[1]+20.0*Y[2]-10010.0*Y[3]+T1;
end;
end.
Unit fjde94r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fjde94r(X :Real; var Y :Array of Real; var DF :Array of Real;
M :Integer);
implementation
procedure fjde94r(X :Real; var Y :Array of Real; var DF :Array of Real;
M :Integer);
begin
DF[0] := -100.0;
DF[4] := 0.0;
DF[8] := 0.0;
DF[12] := 0.0;
DF[1] := -100.0;
DF[5] := -2.0;
DF[9] := 0.0;
DF[13] := 0.0;
DF[2] := -100.0;
DF[6] := 9998.0;
DF[10] := -9990.0;
DF[14] := -10.0;
DF[3] := -100.0;
DF[7] := 9988.0;
DF[11] := 20.0;
DF[15] := -10010.0;
end;
end.
Unit fxde94r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fxde94r(X :Real; var Y :Array of Real; var DX :Array of Real;
M :Integer);
implementation
procedure fxde94r(X :Real; var Y :Array of Real; var DX :Array of Real;
M :Integer);
begin
DX[0] := 0.0;
DX[1] := -2.E3*Exp(-1.E2*X)-2.0*Exp(-X)*(Cos(X)+Sin(X));
DX[2] := DX[1];
DX[3] := DX[1];
end;
end.
Результаты:
после первого обращения к подпрограмме -
Y1 Y2
9.997558891759 + 00 1.099755889174 + 01
Y3 Y4
1.085857138761 + 02 1.085880963993 + 02
YX(1) YX(2)
9.997558891700 + 00 1.099755889169 + 01
YX(3) YX(4)
1.085857138677 + 02 1.085880963909 + 02
X H
2.441406250001 - 06 2.441406250001 - 06
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 9.701276818918 - 13 - 2.910383045670 - 12
RAB(3) RAB(4)
7.823109626764 - 09 7.784304519497 - 09
после второго обращения к подпрограмме -
Y1 Y2
9.997548894156 + 00 1.099754889414 + 01
Y3 Y4
1.085759455506 + 02 1.085783375935 + 02
YX(1) YX(2)
9.997548894098 + 00 1.099754889408 + 01
YX(3) YX(4)
1.085759455420 + 02 1.085783375847 + 02
X H
2.451406249999 - 06 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
0.000000000000 + 00 0.000000000000 + 00
RAB(3) RAB(4)
- 7.761021455135 - 12 - 1.552204291027 - 11
после третьего обращения к подпрограмме -
Y1 Y2
9.997558891759 + 00 1.099755889174 + 01
Y3 Y4
1.085857138761 + 02 1.085880963993 + 02
YX(1) YX(2)
9.997558891613 + 00 1.099755889160 + 01
YX(3) YX(4)
1.085857138667 + 02 1.085880963901 + 02
X H
2.441406249998 - 06 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 9.701276818918 - 13 - 9.701276818918 - 13
RAB(3) RAB(4)
- 2.328306436536 - 11 - 1.552204291027 - 11
после четвертого обращения к подпрограмме -
Y1 Y2
9.997568889303 + 00 1.099756888929 + 01
Y3 Y4
1.085954831772 + 02 1.085978561790 + 02
YX(1) YX(2)
9.997568889157 + 00 1.099756888914 + 01
YX(3) YX(4)
1.085954831680 + 02 1.085978561698 + 02
X H
2.431406249996 - 06 - 4.000000000008 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.940255363783 - 12 - 1.940255363783 - 12
RAB(3) RAB(4)
- 1.552204291027 - 11 - 1.552204291027 - 11
после пятого обращения к подпрограмме -
Y1 Y2
9.997558891759 + 00 1.099755889174 + 01
Y3 Y4
1.085857138761 + 02 1.085880963993 + 02
YX(1) YX(2)
9.997558891613 + 00 1.099755889160 + 01
YX(3) YX(4)
1.085857138667 + 02 1.085880963901 + 02
X H
2.441406249998 - 06 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 9.701276818918 - 13 - 9.701276818918 - 13
RAB(3) RAB(4)
- 2.328306436536 - 11 - 1.552204291027 - 11
после шестого обращения к подпрограмме -
Y1 Y2
9.995118379389 + 00 1.099511837936 + 01
Y3 Y4
1.062995982552 + 02 1.062342483763 + 02
YX(1) YX(2)
9.995118379200 + 00 1.099511837916 + 01
YX(3) YX(4)
1.062995982379 + 02 1.062342483588 + 02
X H
4.882812499996 - 06 2.441406250001 - 06
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.940255363783 - 12 - 2.910387045670 - 12
RAB(3) RAB(4)
7.644606133297 - 09 7.590278983110 - 09
Пример 2.
Решается задача Коши из примера 1. Приводятся фрагмент вызывающей
программы для DE96R, подпрограмма F и результаты счета,
полученные после нескольких обращений к подпрограмме DE96R.
Unit tde96r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE96R_p, DE96R_p;
function tde96r: String;
implementation
function tde96r: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,Y1,Y2,Y3,Y4,XP :Real;
BUL :Boolean;
YX :Array [0..3] of Real;
YP :Array [0..3] of Real;
A :Array [0..15] of Real;
A1 :Array [0..15] of Real;
A2 :Array [0..15] of Real;
RAB :Array [0..19] of Real;
R1 :Array [0..3] of Real;
R2 :Array [0..3] of Real;
R3 :Array [0..3] of Real;
R4 :Array [0..3] of Real;
YD :Array [0..3] of Real;
label
_100,_101,_102,_103,_104,_105,_106;
begin
Result := ''; { результат функции }
M := 4;
X := 0.0;
YX[0] := 10.0;
YX[1] := 11.0;
YX[2] := 111.0;
YX[3] := 111.0;
HMIN := 1.E-10;
EPS := 1.E-4;
P := 1.E3;
JSTART := 0;
H := 0.01;
IH := 0;
_100:
IH := IH + 1;
DE96R (FDE96R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,
A,A1,A2,RAB,R1,R2,R3,R4,YD,IERR);
{ BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y1 := 10.0*Exp(-100.0*X);
Y2 := Y1+Exp(-X)*(Cos(X)+Sin(X));
Y3 := Y2+100.0*Exp(-1.E4*X)*Cos(10.0*X);
Y4 := Y3+100.0*Exp(-1.E4*X)*Sin(10.0*X);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y1,Y2]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y3,Y4]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[YX[0],YX[1]]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[YX[2],YX[3],X,H]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',
[' PACПEЧATKA ПOГPEШHOCTИ ПPИБЛИЖEHHOГO PEШEHИЯ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
{ PACПEЧATKA ПOГPEШHOCTИ ПPИБЛИЖEHHOГO PEШEHИЯ }
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[RAB[0],RAB[1],RAB[2],RAB[3]]) + #$0D#$0A;
case IH of
1: goto _101;
2: goto _102;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _106;
end;
_101:
H := 1.E-8;
goto _100;
_102:
H := -1.E-8;
goto _100;
_103:
JSTART := -1;
goto _100;
_104:
JSTART := -1;
H := -1.E-8;
goto _100;
_105:
H := 0.01;
goto _100;
_106:
UtRes('tde96r',Result); { вывод результатов в файл tde96r.res }
exit;
end;
end.
Unit fde96r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde96r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
implementation
procedure fde96r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
var
T1 :Real;
begin
DY[0] := -100.0*Y[0];
T1 := 20.0*Exp(-100.0*X)+2.0*Exp(-X)*Cos(X);
DY[1] := DY[0]-2.0*Y[1]+T1;
DY[2] := DY[0]+9998.0*Y[1]-9990.0*Y[2]-10.0*Y[3]+T1;
DY[3] := DY[0]+9988.0*Y[1]+20.0*Y[2]-10010.0*Y[3]+T1;
end;
end.
Результаты:
после первого обращения к подпрограмме -
Y1 Y2
9.999961853100 + 00 1.099996185309 + 01
Y3 Y4
1.109618221549 + 02 1.109618602874 + 02
YX(1) YX(2)
9.999961853042 + 00 1.099996185304 + 01
YX(3) YX(4)
1.109618221468 + 02 1.109618602339 + 02
X H
3.814697265626 - 08 3.814697265626 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.940255363783 - 12 - 1.940255363783 - 12
RAB(3) RAB(4)
- 1.777273913224 - 09 3.523503740628 - 09
после второго обращения к подпрограмме -
Y1 Y2
9.999951853111 + 00 1.099995185310 + 01
Y3 Y4
1.109518164690 + 02 1.109518645927 + 02
YX(1) YX(2)
9.999951853039 + 00 1.099995185304 + 01
YX(3) YX(4)
1.109518164597 + 02 1.109518645362 + 02
X H
4.814697265626 - 08 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 9.701276818918 - 13 - 9.701276818918 - 13
RAB(3) RAB(4)
1.862645149229 - 10 - 1.552204291027 - 11
после третьего обращения к подпрограмме -
Y1 Y2
9.999961853100 + 00 1.099996185309 + 01
Y3 Y4
1.109618221549 + 02 1.109618602874 + 02
YX(1) YX(2)
9.999961852940 + 00 1.099996185294 + 01
YX(3) YX(4)
1.109618221460 + 02 1.109618602293 + 02
X H
3.814697265621 - 08 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.940255363783 - 12 - 1.940255363783 - 12
RAB(3) RAB(4)
- 6.208817164108 - 11 6.984919309620 - 11
после четвертого обращения к подпрограмме -
Y1 Y2
9.999971853074 + 00 1.099997185306 + 01
Y3 Y4
1.109718288408 + 02 1.109718569798 + 02
YX(1) YX(2)
9.999971852914 + 00 1.099997185291 + 01
YX(3) YX(4)
1.109718288386 + 02 1.109718569178 + 02
X H
2.814697265622 - 08 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
0.000000000000 + 00 0.000000000000 + 00
RAB(3) RAB(4)
1.396983861924 - 10 3.182018796603 - 10
после пятого обращения к подпрограмме -
Y1 Y2
9.999961853100 + 00 1.099996185309 + 01
Y3 Y4
1.109618221549 + 02 1.109618602874 + 02
YX(1) YX(2)
9.999961852940 + 00 1.099996185294 + 01
YX(3) YX(4)
1.109618221460 + 02 1.109618602293 + 02
X H
3.814697265621 - 08 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.940255363783 - 12 - 1.940255363783 - 12
RAB(3) RAB(4)
- 6.208817164108 - 11 6.984919309620 - 11
после шестого обращения к подпрограмме -
Y1 Y2
9.999923706346 + 00 1.099992370633 + 01
Y3 Y4
1.109236588570 + 02 1.109237350927 + 02
YX(1) YX(2)
9.999923706127 + 00 1.099992370613 + 01
YX(3) YX(4)
1.109236588565 + 02 1.109237350506 + 02
X H
7.629394531243 - 08 3.814697265626 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.940255363783 - 12 - 1.940255363783 - 12
RAB(3) RAB(4)
4.020209113755 - 09 3.554547826448 - 09
Пример 3.
Решается задача Коши
y1' = - 10000y1 + 100y2 - 10y3 + y4
y2' = - 1000y2 + 10y3 -10y4
y3' = - y3 + 10y4
y4' = - 0.1y4
x = 0 , y1(0) = 1 , y2(0) = 1 , y3(0) = 1 , y4(0) = 1
Точное решение задачи имеет вид:
y1 = (1 - d - f - g) t3 + d t2 + f t1 + g y4
y2 = (1 - a - b) t2 + a t1 + b y4
y3 = (- 9.1 / 0.9) t1 + (10 / 0.9) y4
y4 = e - 0.1 x
где t = e - x , t = e - 1000x , t = e - 10000x ,
a = - 91 / (0.9*999) , b = 91 / (0.9*999.9) , d = (1 - a - b) / 90 ,
f = [100 a + (91/0.9)] / 9999 , g = [100 b - (100/0.9) + 1] / 9999.9 .
Матрица Якоби правой части системы имеет вид:
- 10000 100 - 10 1
0 - 1000 10 - 10
0 0 - 1 10
0 0 0 - 0.1
Ниже приводятся фрагмент вызывающей программы для DE90R, подпрограммы F и FJ и результаты счета, полученные после нескольких обращений к подпрограмме DE90R.
Unit tde90r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE90R_p, FJDE90R_p, DE90R_p;
function tde90r: String;
implementation
function tde90r: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,Y4,T1,Y3,T2,AA,BB,C,Y2,T3,DD,FF,GG,QQ,Y1,XP :Real;
BUL :Boolean;
YX :Array [0..3] of Real;
YP :Array [0..3] of Real;
A :Array [0..15] of Real;
A1 :Array [0..15] of Real;
A2 :Array [0..15] of Real;
RAB :Array [0..15] of Real;
R1 :Array [0..3] of Real;
R2 :Array [0..3] of Real;
YD :Array [0..3] of Real;
label
_100,_101,_102,_103,_104,_105,_106;
begin
Result := ''; { результат функции }
M := 4;
X := 0.0;
YX[0] := 1.0;
YX[1] := 1.0;
YX[2] := 1.0;
YX[3] := 1.0;
HMIN := 1.E-10;
EPS := 1.E-4;
P := 1.E3;
JSTART := 0;
H := 0.01;
IH := 0;
_100:
IH := IH + 1;
DE90R (FDE90R,FJDE90R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,A,A1,
A2,RAB,R1,R2,YD,IERR);
{ BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y4 := Exp(-0.1*X);
T1 := Exp(-X);
Y3 := -(9.1/0.9)*T1 + (10.0/0.9)*Y4;
T2 := Exp(-1.E3*X);
AA := -91.0/(0.9*999.0);
BB := 91.0/(0.9*999.9);
C := 1.0-AA-BB;
Y2 := C*T2+AA*T1+BB*Y4;
T3 := Exp(-1.E4*X);
DD := C/90.0;
FF := (100.0*AA+91.0/0.9)/9999.0;
GG := (100.0*BB-100.0/0.9+1.0)/9999.9;
QQ := 1.0-DD-FF-GG;
Y1 := QQ*T3+DD*T2+FF*T1+GG*Y4;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y1,Y2]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y3,Y4]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[YX[0],YX[1]]);
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[YX[2],YX[3],X,H]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',
[' PACПEЧATKA ПOГPEШHOCTИ ПPИБЛИЖEHHOГO PEШEHИЯ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
{ PACПEЧATKA ПOГPEШHOCTИ ПPИБЛИЖEHHOГO PEШEHИЯ }
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[RAB[0],RAB[1],RAB[2],RAB[3]]) + #$0D#$0A;
case IH of
1: goto _101;
2: goto _102;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _106;
end;
_101:
H := 1.E-8;
goto _100;
_102:
EPS := 1.E-2;
goto _100;
_103:
JSTART := -1;
EPS := 1.E-5;
goto _100;
_104:
JSTART := -1;
H := -1.E-8;
goto _100;
_105:
H := 0.01;
goto _100;
_106:
UtRes('tde90r',Result); { вывод результатов в файл tde90r.res }
exit;
end;
end.
Unit fde90r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde90r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
implementation
procedure fde90r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
begin
DY[0] := -1.E4*Y[0]+1.E2*Y[1]-10.0*Y[2]+Y[3];
DY[1] := -1.E3*Y[1]+10.0*Y[2]-10.0*Y[3];
DY[2] := -Y[2]+10.0*Y[3];
DY[3] := -0.1*Y[3];
end;
end.
Unit fjde90r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fjde90r(X :Real; var Y :Array of Real; var DF :Array of Real;
M :Integer);
implementation
procedure fjde90r(X :Real; var Y :Array of Real; var DF :Array of Real;
M :Integer);
begin
DF[0] := -1.E4;
DF[4] := 1.E2;
DF[8] := -10.0;
DF[12] := 1.0;
DF[1] := 0.0;
DF[5] := -1.E3;
DF[9] := 10.0;
DF[13] := -10.0;
DF[2] := 0.0;
DF[6] := 0.0;
DF[10] := -1.0;
DF[14] := 10.0;
DF[3] := 0.0;
DF[7] := 0.0;
DF[11] := 0.0;
DF[15] := -0.1;
end;
end.
Результаты:
после первого обращения к подпрограмме -
Y1 Y2
9.527772901474 - 01 9.951290901317 - 01
Y3 Y4
1.000043945169 + 00 9.999995117187 - 01
YX(1) YX(2)
9.527772877191 - 01 9.951290901299 - 01
YX(3) YX(4)
1.000043945183 + 00 9.999995117150 - 01
X H
4.882812500002 - 06 4.882812500002 - 06
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
2.270887004367 - 09 - 1.212659602364 - 13
RAB(3) RAB(4)
-3.637978807088 - 13 - 1.818989403544 - 13
после второго обращения к подпрограмме -
Y1 Y2
9.526829222541 - 01 9.951191388946 - 01
Y3 Y4
1.000044035131 + 00 9.999995107200 - 01
YX(1) YX(2)
9.526829198239 - 01 9.951191388900 - 01
YX(3) YX(4)
1.000044035176 + 00 9.999995107109 - 01
X H
4.892812499997 - 06 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.818989403544 - 13 - 1.212659602364 - 13
RAB(3) RAB(4)
- 2.425319204729 - 13 - 1.212659602364 - 13
после третьего обращения к подпрограмме -
Y1 Y2
9.524942147509 - 01 9.950992367185 - 01
Y3 Y4
1.000044215194 + 00 9.999995087192 - 01
YX(1) YX(2)
9.524942123180 - 01 9.950992367112 - 01
YX(3) YX(4)
1.000044215169 + 00 9.999995087064 - 01
X H
4.912812499994 - 06 4.000000000008 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.212659602364 - 13 - 1.212659602364 - 13
RAB(3) RAB(4)
- 1.212659602364 - 13 - 1.818989403544 - 13
после четвертого обращения к подпрограмме -
Y1 Y2
9.523055449436 - 01 9.950793349426 - 01
Y3 Y4
1.000044395201 + 00 9.999995067164 - 01
YX(1) YX(2)
9.523055425125 - 01 9.950793349353 - 01
YX(3) YX(4)
1.000044395165 + 00 9.999995067074 - 01
X H
4.932812499997 - 06 8.000000000017 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 6.063298011824 - 14 - 6.063298011824 - 14
RAB(3) RAB(4)
- 3.637978807088 - 13 - 6.063298011824 - 14
после пятого обращения подпрограмме -
Y1 Y2
9.527772901474 - 01 9.951290901317 - 01
Y3 Y4
1.000043945169 + 00 9.999995117187 - 01
YX(1) YX(2)
9.527772877100 - 01 9.951290901245 - 01
YX(3) YX(4)
1.000043945169 + 00 9.999995117078 - 01
X H
4.882812499996 - 06 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.818989403544 - 13 - 6.063298011824 - 14
RAB(3) RAB(4)
- 2.425319204729 - 13 - 1.212659602364 - 13
после шестого обращения к подпрограмме -
Y1 Y2
9.078026703492 - 01 9.902819081972 - 01
Y3 Y4
1.000087890148 + 00 9.999990234382 - 01
YX(1) YX(2)
9.078026657144 - 01 9.902819081881 - 01
YX(3) YX(4)
1.000087890114 + 00 9.999990234237 - 01
X H
9.765624999991 - 06 4.882812500002 - 06
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
2.162717767836 - 09 0.000000000000 + 00
RAB(3) RAB(4)
- 3.637978807088 - 13 - 6.063298011824 - 14
Пример 4.
Решается задача Коши из примера 3. Приводятся фрагмент вызывающей
программы для DE92R, подпрограмма F и результаты счета, полученные
после нескольких обращений к подпрограмме DE92R.
Unit tde92r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE92R_p, DE92R_p;
function tde92r: String;
implementation
function tde92r: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,Y4,T1,Y3,T2,AA,BB,C,Y2,T3,DD,FF,GG,QQ,Y1,XP :Real;
BUL :Boolean;
YX :Array [0..3] of Real;
YP :Array [0..3] of Real;
A :Array [0..15] of Real;
A1 :Array [0..15] of Real;
A2 :Array [0..15] of Real;
RAB :Array [0..15] of Real;
R1 :Array [0..3] of Real;
R2 :Array [0..3] of Real;
YD :Array [0..3] of Real;
label
_100,_101,_102,_103,_104,_105,_106;
begin
Result := ''; { результат функции }
M := 4;
X := 0.0;
YX[0] := 1.0;
YX[1] := 1.0;
YX[2] := 1.0;
YX[3] := 1.0;
HMIN := 1.E-10;
EPS := 1.E-4;
P := 1.E3;
JSTART := 0;
H := 0.01;
IH := 0;
_100:
IH := IH + 1;
DE92R (FDE92R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,A,A1,
A2,RAB,R1,R2,YD,IERR);
{ BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y4 := Exp(-0.1*X);
T1 := Exp(-X);
Y3 := -(9.1/0.9)*T1 + (10.0/0.9)*Y4;
T2 := Exp(-1.E3*X);
AA := -91.0/(0.9*999.0);
BB := 91.0/(0.9*999.9);
C := 1.0-AA-BB;
Y2 := C*T2+AA*T1+BB*Y4;
T3 := Exp(-1.E4*X);
DD := C/90.0;
FF := (100.0*AA+91.0/0.9)/9999.0;
GG := (100.0*BB-100.0/0.9+1.0)/9999.9;
QQ := 1.0-DD-FF-GG;
Y1 := QQ*T3+DD*T2+FF*T1+GG*Y4;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y1,Y2]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[Y3,Y4]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[YX[0],YX[1]]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[YX[2],YX[3],X,H]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
Result := Result + Format('%s',
[' PACПEЧATKA ПOГPEШHOCTИ ПPИБЛИЖEHHOГO PEШEHИЯ']) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A;
{ PACПEЧATKA ПOГPEШHOCTИ ПPИБЛИЖEHHOГO PEШEHИЯ }
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[RAB[0],RAB[1],RAB[2],RAB[3]]) + #$0D#$0A;
case IH of
1: goto _101;
2: goto _102;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _106;
end;
_101:
H := 1.E-8;
goto _100;
_102:
EPS := 1.E-2;
goto _100;
_103:
JSTART := -1;
EPS := 1.E-4;
goto _100;
_104:
JSTART := -1;
H := -1.E-8;
goto _100;
_105:
H := 0.01;
goto _100;
_106:
UtRes('tde92r',Result); { вывод результатов в файл tde92r.res }
exit;
end;
end.
Unit fde92r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde92r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
implementation
procedure fde92r(X :Real; var Y :Array of Real; var DY :Array of Real;
M :Integer);
begin
DY[0] := -1.E4*Y[0]+1.E2*Y[1]-10.0*Y[2]+Y[3];
DY[1] := -1.E3*Y[1]+10.0*Y[2]-10.0*Y[3];
DY[2] := -Y[2]+10.0*Y[3];
DY[3] := -0.1*Y[3];
end;
end.
Результаты:
после первого обращения к подпрограмме -
Y1 Y2
9.527772901474 - 01 9.951298901317 - 01
Y3 Y4
1.000043945169 + 00 9.999995117187 - 01
YX(1) YX(2)
9.527772906931 - 01 9.951290903537 - 01
YX(3) YX(4)
1.000043945183 + 00 9.999995117150 - 01
X H
4.882812500002 - 06 4.882812500002 - 06
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 7.571846557159 - 10 - 1.412748436753 - 11
RAB(3) RAB(4)
- 3.637978807088 - 13 - 1.818989403544 - 13
после второго обращения к подпрограмме -
Y1 Y2
9.526829222541 - 01 9.951191388946 - 01
Y3 Y4
1.000044035231 + 00 9.999995107200 - 01
YX(1) YX(2)
9.526829227989 - 01 9.951191391137 - 01
YX(3) YX(4)
1.000044035176 + 00 9.999995107109 - 01
X H
4.892812499997 - 06 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.212259602364 - 13 - 1.212659602364 - 13
RAB(3) RAB(4)
- 2.425319204729 - 13 - 1.212659602364 - 13
после третьего обращения к подпрограмме -
Y1 Y2
9.524942147509 - 01 9.950992367185 - 01
Y3 Y4
1.000044215194 + 00 9.999995087192 - 01
YX(1) YX(2)
9.524942152921 - 01 9.950992369349 - 01
YX(3) YX(4)
1.000044215169 + 00 9.999995087064 - 01
X H
4.912812499994 - 06 4.000000000008 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.212659602364 - 13 - 1.212659602364 - 13
RAB(3) RAB(4)
- 1.212659602364 - 13 - 1.818989403544 - 13
после четвертого обращения к подпрограмме -
Y1 Y2
9.523055449436 - 01 9.950793349426 - 01
Y3 Y4
1.000044395201 + 01 9.999995067164 - 01
YX(1) YX(2)
9.523055454865 - 01 9.950793351591 - 01
YX(3) YX(4)
1.000044395165 + 00 9.999995067074 - 01
X H
4.932812499997 - 06 8.000000000017 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 6.063298011814 - 13 0.000000000000 + 00
RAB(3) RAB(4)
- 3.637978807088 - 13 - 6.063298011824 - 14
после пятого обращения к подпрограмме -
Y1 Y2
9.527772901474 - 01 9.951290901317 - 01
Y3 Y4
1.000043945169 + 00 9.999995117187 - 01
YX(1) YX(2)
9.527772906868 - 01 9.951290903482 - 01
YX(3) YX(4)
1.000043945169 + 00 9.999995117078 - 01
X H
4.882812499996 - 06 - 2.000000000004 - 08
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
- 1.212659602364 - 13 - 6.063298011824 - 14
RAB(3) RAB(4)
- 2.425319204729 - 13 - 1.212659602364 - 13
после шестого обращения к подпрограмме -
Y1 Y2
9.078026703492 - 01 9.902819081972 - 01
Y3 Y4
1.000087890148 + 00 9.999990234382 - 01
YX(1) YX(2)
9.078026512007 - 01 9.902819086356 - 01
YX(3) YX(4)
1.000087890114 + 00 9.999990234237 - 01
X H
9.765624999991 - 06 4.882812500002 - 06
Распечатка контрольного члена погрешности:
RAB(1) RAB(2)
3.313774262397 - 09 - 1.552204291027 - 11
RAB(3) RAB(4)
- 3.637978807088 - 13 - 6.063298011824 - 14