1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
{ ******************************************************************
Gauss-Legendre integration
****************************************************************** }
unit ugausleg;
interface
uses
utypes;
function GausLeg(Func : TFunc; A, B : Float) : Float;
{ Integral from A to B }
function GausLeg0(Func : TFunc; B : Float) : Float;
{ Integral from 0 to B }
function Convol(Func1, Func2 : TFunc; T : Float) : Float;
{ Convolution product at time T }
implementation
const Npts = 8; { Number of points / 2 }
const Root : array[1..Npts] of Float =
(0.0950125098376370440185,
0.281603550778258913230,
0.458016777657227386342,
0.617876244402643748447,
0.755404408355003033895,
0.865631202387831743880,
0.944575023073232576078,
0.989400934991649932596);
const Weight : array[1..Npts] of Float =
(0.189450610455068496285,
0.182603415044923588867,
0.169156519395002538189,
0.149595988816576732081,
0.124628971255533872052,
0.095158511682492784810,
0.062253523938647892863,
0.027152459411754094852);
function GausLeg(Func : TFunc; A, B : Float) : Float;
{ ------------------------------------------------------------------
Computes the integral of function Func from A to B
by the Gauss-Legendre method
------------------------------------------------------------------ }
var
P, Q, Sum, X, Y : Float;
I : Integer;
begin
P := 0.5 * (B + A);
Q := 0.5 * (B - A);
Sum := 0.0;
for I := 1 to Npts do
begin
X := Q * Root[I];
Y := Func(P + X) + Func(P - X);
Sum := Sum + Weight[I] * Y;
end;
GausLeg := Q * Sum;
end;
function GausLeg0(Func : TFunc; B : Float) : Float;
{ ------------------------------------------------------------------
Computes the integral of function Func from 0 to B
by the Gauss-Legendre method
------------------------------------------------------------------ }
var
P, Sum, X, Y : Float;
I : Integer;
begin
P := 0.5 * B;
Sum := 0;
for I := 1 to Npts do
begin
X := P * Root[I];
Y := Func(P + X) + Func(P - X);
Sum := Sum + Weight[I] * Y;
end;
GausLeg0 := P * Sum;
end;
function Convol(Func1, Func2 : TFunc; T : Float) : Float;
{ ------------------------------------------------------------------
Computes the convolution product of two functions Func1 and Func2
at time T by the Gauss-Legendre method
------------------------------------------------------------------ }
var
P, PpX, PmX, Sum, X, Y : Float;
I : Integer;
begin
P := 0.5 * T;
Sum := 0.0;
for I := 1 to Npts do
begin
X := P * Root[I];
PpX := P + X;
PmX := P - X;
Y := Func1(T - PpX) * Func2(PpX) + Func1(T - PmX) * Func2(PmX);
Sum := Sum + Weight[I] * Y;
end;
Convol := P * Sum;
end;
end.
|