LORENE
integrate_1D.C
1
/*
2
* Integration of f(x) in the interval [xx(0), xx(n-1)], with non-equally spaced
3
* n-size xx grid.
4
*
5
* The function f is approximated by piecewise parabolae, The integral of f
6
* is set to 0 at xx(0).
7
*/
8
9
/*
10
* Copyright (c) 2015 Jerome Novak
11
*
12
* This file is part of LORENE.
13
*
14
* LORENE is free software; you can redistribute it and/or modify
15
* it under the terms of the GNU General Public License as published by
16
* the Free Software Foundation; either version 2 of the License, or
17
* (at your option) any later version.
18
*
19
* LORENE is distributed in the hope that it will be useful,
20
* but WITHOUT ANY WARRANTY; without even the implied warranty of
21
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22
* GNU General Public License for more details.
23
*
24
* You should have received a copy of the GNU General Public License
25
* along with LORENE; if not, write to the Free Software
26
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27
*
28
*/
29
30
31
char
integrate_1D_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/integrate_1D.C,v 1.1 2015/01/09 15:28:52 j_novak Exp $"
;
32
33
/*
34
* $Id: integrate_1D.C,v 1.1 2015/01/09 15:28:52 j_novak Exp $
35
* $Log: integrate_1D.C,v $
36
* Revision 1.1 2015/01/09 15:28:52 j_novak
37
* New integration function for general non-equally-spaced grids.
38
*
39
*
40
*/
41
42
// Headers Lorene
43
#include "tbl.h"
44
45
namespace
Lorene
{
46
47
Tbl
integ1D
(
const
Tbl
&
xx
,
const
Tbl
& ff) {
48
49
Tbl
resu
(ff) ;
50
if
(ff.get_etat() != ETATZERO) {
51
52
assert
(
xx
.get_etat() == ETATQCQ) ;
53
assert
(ff.get_etat() == ETATQCQ) ;
54
int
nx
=
xx
.get_taille() ;
55
assert
(
nx
> 2) ;
56
assert
(ff.get_taille() ==
nx
) ;
57
58
resu
.set(0) = 0. ;
59
double
x0
=
xx
(0) ;
60
double
x1
(0),
x2
(0),
x3
(0);
61
double
a1
(0),
a2
(0),
a3
(0);
62
double
b1
(0),
b2
(0),
b3
(0);
63
double
c1
(0),
c2
(0),
c3
(0) ;
64
65
for
(
int
i
=1;
i
<
nx
-1;
i
++) {
66
x1
=
xx
(
i
-1) ;
67
x2
=
xx
(
i
) ;
68
x3
=
xx
(
i
+1) ;
69
a1
= ff(
i
-1) / ( (
x1
-
x2
)*(
x1
-
x3
) ) ;
70
a2
= ff(
i
) / ( (
x2
-
x1
)*(
x2
-
x3
) ) ;
71
a3
= ff(
i
+1) / ( (
x3
-
x1
)*(
x3
-
x2
) ) ;
72
b1
=
a1
+
a2
+
a3
;
73
b2
= -(
x2
+
x3
)*
a1
- (
x1
+
x3
)*
a2
- (
x1
+
x2
)*
a3
;
74
b3
=
x2
*
x3
*
a1
+
x1
*
x3
*
a2
+
x1
*
x2
*
a3
;
75
if
(
i
==1) {
76
c1
=
b1
;
77
c2
=
b2
;
78
c3
=
b3
;
79
}
80
else
{
81
c1
= 0.5*(
b1
+
c1
) ;
82
c2
= 0.5*(
b2
+
c2
) ;
83
c3
= 0.5*(
b3
+
c3
) ;
84
}
85
resu
.set(
i
) =
resu
(
i
-1) +
c1
*(
x2
*
x2
*
x2
-
x0
*
x0
*
x0
)/3.
86
+ 0.5*
c2
*(
x2
*
x2
-
x0
*
x0
) +
c3
*(
x2
-
x0
) ;
87
c1
=
b1
;
88
c2
=
b2
;
89
c3
=
b3
;
90
x0
=
x2
;
91
}
92
93
x2
=
xx
(
nx
-1) ;
94
resu
.set(
nx
-1) =
resu
(
nx
-2) +
c1
*(
x2
*
x2
*
x2
-
x0
*
x0
*
x0
)/3.
95
+ 0.5*
c2
*(
x2
*
x2
-
x0
*
x0
) +
c3
*(
x2
-
x0
) ;
96
97
}
98
return
resu
;
99
}
100
101
}
// End of namespace Lorene
Lorene::Evolution_std
Time evolution with partial storage (*** under development ***).
Definition
evolution.h:371
Lorene::Tbl
Basic array class.
Definition
tbl.h:161
Lorene::integ1D
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piece parabolae.
Definition
integrate_1D.C:47
Lorene
Lorene prototypes.
Definition
app_hor.h:64
C++
Source
Non_class_members
Utilities
integrate_1D.C
Generated by
1.9.8