|
Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://num-anal.srcc.msu.su/lib_na/cat/q_htm_p/qsgce_p.htm
Дата изменения: Tue Nov 17 12:59:15 2015 Дата индексирования: Sun Apr 10 01:13:06 2016 Кодировка: Windows-1251 |
|
Текст подпрограммы и версий qsgce_p.zip |
Тексты тестовых примеров tqsgce_p.zip tqsgce1_p.zip |
Быстрое вычисление узлов и весов квадратурной формулы Гаусса с расширенной (Extended) точностью.
Вычисляются узлы X ( k ) и веса C ( k ), k = 1, 2, ..., N квадратуры Гаусса на отрезке [- 1, 1] при любом N (2 ≤ N < ∞) за O ( N ) операций с расширенной (Extended) точностью.
Л.Г. Васильева, Я.М. Жилейкин, О быстром вычислении узлов и весов квадратуры Гаусса, ЖВМ и МФ, 2003.
procedure qsgce(N :Integer; N1 :Integer; var X :Array of Extended;
var C :Array of Extended);
Параметры
| N - | число узлов квадратуры Гаусса (тип: целый); |
| N1 - |
N1 = N / 2, если N - четное, N1 = ( N + 1) / 2, если N - нечетное, (тип: целый); |
| X - | вектор длины N1 для вычисленных узлов квадратуры Гаусса (тип: вещественный); |
| C - | вектор длины N1 для вычисленных весов квадратуры Гаусса (тип: вещественный). |
Вызываемые подпрограммы: нет
Замечания по использованию
| 1. |
Требуется, чтобы N ≥ 2. | |
| 2. |
Узлы квадратур Гаусса располагаются симметрично на отрезке [- 1, 1] относительно начала координат: X ( k ) = - X ( N - k + 1), а веса в симметричных точках равны C ( k ) = C ( N - k + 1), k = 1, 2, ..., N1. Поэтому подпрограмма QSGCD вычисляет N1 узлов и весов Гаусса только для отрезка [0, 1]. | |
| 3. |
Узлы и веса X ( k ) и C ( k ) вычисляются с точностью 10- 13. |
1.
Unit tqsgce_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, QSGCE_p;
function tqsgce: String;
implementation
function tqsgce: String;
var
N,N1,I :Integer;
X :Array [0..2] of Extended;
C :Array [0..2] of Extended;
begin
Result := ''; { результат функции }
{ ЕСЛИ N - ЧЕТНОЕ, ТО N1=N/2; ЕСЛИ N - НЕЧЕТНОЕ, ТО N1=(N+1)/2 }
N := 5;
N1 := 3;
QSGCE(5,3,X,C);
Result := Result + Format('%s',[' EPS=1.D-15' + #$0D#$0A]) + #$0D#$0A;
Result := Result + Format('%s',[' N=']);
Result := Result + Format('%5d ',[N]) + #$0D#$0A;
Result := Result + Format('%s',
[' МАССИВ УЗЛОВ Х(N1)' + #$0D#$0A]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to N1 do
begin
Result := Result + Format('%20.16f ',[X[I-1]]) + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',
[' МАССИВ ВЕСОВ С(N1)' + #$0D#$0A]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to N1 do
begin
Result := Result + Format('%20.16f ',[C[I-1]]) + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
UtRes('tqsgce',Result); { вывод результатов в файл tqsgce.res }
end;
end.
Результаты:
X(1) = .906179845938664 C(1) = .236926885056189
X(2) = .538469310105683 C(2) = .478628670499366
X(3) = .000000000000000 C(3) = .568888888888889
2.
Unit tqsgce1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, QSGCE_p;
function tqsgce1: String;
implementation
function tqsgce1: String;
var
N,N1,I :Integer;
X :Array [0..255] of Extended;
C :Array [0..255] of Extended;
begin
Result := ''; { результат функции }
{ ЕСЛИ N - ЧЕТНОЕ, ТО N1=N/2; ЕСЛИ N - НЕЧЕТНОЕ, ТО N1=(N+1)/2 }
N := 512;
N1 := 256;
QSGCE(N,N1,X,C);
Result := Result + Format('%s',[' EPS=1.D-15' + #$0D#$0A]) + #$0D#$0A;
Result := Result + Format('%s',[' N=']);
Result := Result + Format('%5d ',[N]) + #$0D#$0A;
Result := Result + Format('%s',
[' МАССИВ УЗЛОВ Х(N1)' + #$0D#$0A]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to N1 do
begin
Result := Result + Format('%20.16f ',[X[I-1]]) + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',
[' МАССИВ ВЕСОВ С(N1)' + #$0D#$0A]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to N1 do
begin
Result := Result + Format('%20.16f ',[C[I-1]]) + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
UtRes('tqsgce1',Result); { вывод результатов в файл tqsgce1.res }
exit;
end;
end.
Результаты:
X(1) = .999988990984382 C(1) = .000028252637374
X(2) = .999941994606846 C(2) = .000065765731660
......................................................................................................
X(255) = .009194771386433 C(255) = .006129674838036
X(256) = .003064962185159 C(256) = .006129905175406