Commit 83b8b60440beb64b30ed08642d2c5ddd1c1c82ff

  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Committer)
  • Tue May 17 18:45:45 EEST 2011
  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Author)
  • Tue May 17 18:45:45 EEST 2011
Class to handle vectors written in Chebyshev basis
Makefile
(2 / 1)
  
22LDFLAGS = `gsl-config --libs` -lm -L /linuxfs-home/student/hejajama/lib/lib
33
44SOURCES = src/main.cpp src/amplitude.cpp src/tools.cpp src/datafile.cpp \
5 src/solver_force.cpp src/solver_chebyshev.cpp src/chebyshev_amplitude.cpp
5 src/solver_force.cpp src/solver_chebyshev.cpp src/chebyshev_amplitude.cpp \
6 src/chebyshev.cpp
67FFTSOURCES = src/fft4g.c
78LIBBCISOURCES = #libbci-1.1.0/bci.c libbci-1.1.0/tdspl.c libbci-1.1.0/tools.c
89OBJECTS=$(SOURCES:.cpp=.o)
src/chebyshev.cpp
(108 / 0)
  
1/*
2 * BK equation solver
3 * Heikki Mäntysaari <heikki.mantysaari@jyu.fi>, 2011
4 */
5
6#include "tools.hpp"
7#include "chebyshev.hpp"
8#include <gsl/gsl_chebyshev.h>
9#include <gsl/gsl_integration.h>
10#include <cmath>
11
12const unsigned int MAXITER_INNERPROD = 1000;
13const REAL INNERPROD_ACCURACY = 0.0001;
14
15/*
16 * Compute the inner product between two Chebyshevs
17 * E.g. integrate the product function over [-1,1] with weight
18 * (1-x^2)^{-1/2}
19 */
20struct Inthelper_innerprod
21{
22 ChebyshevBasis* a;
23 ChebyshevBasis* b;
24
25
26};
27
28REAL Inthelperf_innerprod(REAL x, void* p)
29{
30 Inthelper_innerprod *par = (Inthelper_innerprod*) p;
31 return 1.0/std::sqrt(1.0-SQR(x))*par->a->Evaluate(x)*par->b->Evaluate(x);
32}
33
34REAL ChebyshevBasis::InnerProduct(ChebyshevBasis &vec)
35{
36 Inthelper_innerprod helper;
37 helper.a=this;
38 helper.b=&vec;
39
40 gsl_function int_helper;
41 int_helper.function=&Inthelperf_innerprod;
42 int_helper.params=&helper;
43 REAL result, abserr;
44 /*
45 gsl_integration_workspace *workspace
46 = gsl_integration_workspace_alloc(MAXITER_INNERPROD);
47 int status=gsl_integration_qag(&int_helper, -1, 1, 0,
48 INNERPROD_ACCURACY, MAXITER_INNERPROD,
49 GSL_INTEG_GAUSS51, workspace, &result, &abserr);
50 gsl_integration_workspace_free(workspace);
51 */
52 gsl_integration_cquad_workspace * workspace =gsl_integration_cquad_workspace_alloc(200);
53 int status =gsl_integration_cquad(&int_helper, -1, 1, 0, INNERPROD_ACCURACY,
54 workspace, &result, &abserr, NULL);
55 gsl_integration_cquad_workspace_free(workspace);
56
57 return result;
58
59
60}
61
62REAL ChebyshevBasis::Evaluate(REAL x)
63{
64 if (std::abs(x)>1)
65 {
66 cerr << "Asked the value of Chebyshev polynomial at x=" << x << " "
67 << LINEINFO << endl;
68 return 0;
69 }
70
71 return gsl_cheb_eval(cheb, x);
72
73}
74
75ChebyshevBasis::ChebyshevBasis(unsigned int d)
76{
77 degree = d;
78 cheb = gsl_cheb_alloc (degree);
79 cheb->a=-1.0;
80 cheb->b=1.0;
81 cheb->order=degree;
82}
83
84REAL ChebyshevBasis::Component(unsigned int n)
85{
86 if (n>degree)
87 {
88 cerr<< "Asked chebyshev-vector component " << n << ", but "
89 << "dimension is " << degree << " at " << LINEINFO << endl;
90 return 0.0;
91 }
92
93 return cheb->c[n];
94}
95
96void ChebyshevBasis::SetComponent(unsigned int c, REAL val)
97{
98 if (c>degree)
99 {
100 cerr<< "Tried to set component " << c << " to chebyshev vector but the "
101 << "dimension is " << degree << " at " << LINEINFO << endl;
102 return;
103 }
104
105 if (c==0) val*=2.0;
106 cheb->c[c]=val;
107
108}
src/chebyshev.hpp
(37 / 0)
  
1/*
2 * BK equation solver
3 * Heikki Mäntysaari <heikki.mantysaari@jyu.fi>, 2011
4 */
5
6#ifndef _CHEBYSHEV_H
7#define _CHEBYSHEV_H
8
9/*
10 * Chebyshev polynomial based basis
11 * This class represents a vector for which the components are coefficients of
12 * the Chebyshev polynomials
13 * In that space we define an inner product
14 * a \cdot b = \int_{-1}^1 (1-x)^{-1/2} a(x)b(x)
15 */
16
17#include "config.hpp"
18#include <gsl/gsl_chebyshev.h>
19#include <vector>
20
21class ChebyshevBasis
22{
23 public:
24 ChebyshevBasis(unsigned int d);
25 REAL Component(unsigned int n);
26 REAL InnerProduct(ChebyshevBasis &vec);
27 REAL Evaluate(REAL x);
28 void SetComponent(unsigned int c, REAL val);
29
30
31 private:
32 unsigned int degree;
33 gsl_cheb_series *cheb;
34
35};
36
37#endif
src/main.cpp
(2 / 1)
  
1010#include "solver_chebyshev.hpp"
1111#include "chebyshev_amplitude.hpp"
1212#include "tools.hpp"
13#include "chebyshev.hpp"
1314#include <gsl/gsl_errno.h>
1415#include <cmath>
1516#include <fstream>
3737};
3838
3939int main(int argc, char* argv[])
40{
40{
4141 ChebyshevAmplitudeSolver N;
4242 //BruteForceSolver N;
4343 REAL minktsqr=DEFAULT_MINKTSQR;