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
31char 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
45namespace Lorene {
46
47Tbl 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
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic array class.
Definition tbl.h:161
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piece parabolae.
Lorene prototypes.
Definition app_hor.h:64