Текст подпрограммы и версий qss2r_p.zip qss2e_p.zip |
Тексты тестовых примеров tqss2r_p.zip tqss2e_p.zip |
Вычисление определенного интеграла от векторной функции по обобщенной квадратурной формуле Симпсона.
Подпрограмма QSS2R вычисляет интеграл B ∫ F(x) dx от векторной функции F(x) = ( f1 (x), ... , fN (x) ) A с покомпонентной погрешностью B ЕРS ( 1 + | ∫ fi (x) dx | ) A
по обобщенной квадратурной формуле Симпсона с шагом h = (B - A)/M, где M - число частичных отрезков разбиения.
Я.М.Жилейкин, М.В.Соколовский. Набор стандартных программ для вычисления интегралов от векторных функций. Сб. "Методы и алгоритмы в численном анализе", Изд-во МГУ, 1981.
procedure QSS2R(var RINT :Array of Real; A :Real; B :Real; F :Proc_F3_A; N :Integer; EPS :Real; K :Integer; var IND :Array of Integer; var R :Array of Real; var IERR :Integer);
Параметры
RINT - | вещественный вектоp длины N, содержащий вычисленные значения интегралов; |
A, B - | заданные нижний и верхний пределы интегрирования (тип: вещественный); |
F - | имя подпрограммы, вычисляющей подинтегральные функции; |
N - | размерность вектоpа - функции F (x) (тип: целый); |
EPS - | заданная меpа погрешности вычисления интеграла (тип: вещественный); |
K - | целая переменная, задающая начальное число частичных отрезков разбиения; |
IND - | целый вектоp длины N, каждая компонента которого pавна числу частичных отрезков разбиения, при котоpом достигается заданная точность; |
R - | вещественный рабочий вектоp длины 3N; |
IERR - | целая переменная, служащая для диагностических сообщений: |
IERR = 65 - | когда заданная точность не может быть достигнута при максимально возможном числе (1048576) частичных отрезков разбиения. |
Версии
QSS2E - | вычисление с расширенной (Extended) точностью определенного интеграла от векторной функции по обобщенной квадратурной формуле Симпсона. |
Вызываемые подпрограммы
UTQS10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы QSS2R. |
UTQS12 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы QSS2E. |
Замечания по использованию
1. |
Первый оператор подпрограммы F должен иметь вид: | |
2. |
При EPS ≤ 0 происходит только один просчет при заданном значении K. | |
3. |
Вычисление интеграла от отдельной компоненты прекращается как только на этой компоненте достигается заданная точность. | |
4. |
Если для каких - либо компонент векторной функции заданная погрешность интегрирования не может быть достигнута, то соответствующие компоненты вектоpа IND полагаются равными нулю. | |
5. |
Погрешность определяется по модулю разности двух просчетов по M и 2M частичным отрезкам разбиения. Если заданная точность не достигнута, то значение M удваивается. | |
6. | В подпрограмме QSS2E параметры RINT, A, B, F, EPS, R имеют тип Extended. |
Unit TQSS2R_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc, UtRes_p, FQSS2R_p, QSS2R_p; function TQSS2R: String; implementation function TQSS2R: String; var _i,IERR :Integer; RINT :Array [0..4] of Real; R :Array [0..14] of Real; IND :Array [0..4] of Integer; begin Result := ''; { результат функции } QSS2R(RINT,0.0,1.0,FQSS2R,5,1.E-6,1,IND,R,IERR); RЕSULt := Result + Format('%s',[' RINT=']); Result := Result + #$0D#$0A; for _i:=0 to 4 do begin Result := Result + Format('%20.16f ',[RINT[_i]]); if ( ((_i+1) mod 4)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + Format('%s',[' IND=']); Result := Result + #$0D#$0A; for _i:=0 to 4 do begin Result := Result + Format('%10d ',[IND[_i]]); if ( ((_i+1) mod 4)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' IERR=']); Result := Result + Format('%5d ',[IERR]) + #$0D#$0A; UtRes('TQSS2R',Result); { вывод результатов в файл TQSS2R.res } exit; end; end. Unit FQSS2R_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc; procedure FQSS2R(X :Real; var Y :Array of Real; I :Array of Integer); implementation procedure FQSS2R(X :Real; var Y :Array of Real; I :Array of Integer); var J :Integer; label _1; begin Y[0] := X; for J:=1 to 4 do begin _1: Y[J+0] := Y[J-1]*X; end; end; end. Результаты: RINT IND 5.0000000000-01 1 3.3333333333-01 1 2.5000000000-01 1 2.0000000794-01 16 1.6666668653-01 16 IERR = 0