Various utilities#

This chapter provides informations on tools not used by the FE machinery but may be of interest for end users.

Integral calculation tools#

Besides quadrature rules used by FE stuff, XLiFE++ provides some integration tools.

Uniform rectangle, trapeze and Simpson methods#

To approximate 1D integrals:

\[\int_a^b f(t) dt.\]

XLiFE++ handles the uniform rectangle, trapeze, Simpson method (\(h=\dfrac{b-a}{n}\)):

rectangles method:
\[\int_a^b f(t) dt \sim h\sum_{i=0}^{n-1} f(a+ih)\]
trapezes method:
\[\int_a^b f(t) dt \sim \dfrac{h}{2} \left(f(a)+4\sum_{i=1}^{n/2} f(a+(2i-1)h) +2\sum_{i=1}^{n/2-1} f(a+2ih)+f(b)\right)\ \ (n \mathrm{\;even})\]
Simpson method:
\[\int_a^b f(t) dt \sim \dfrac{h}{3} \left(f(a)+2\sum_{i=1}^{n-1} f(a+ih) + f(b)\right)\]

with \(f\) a real or complex function. They respectively integrate exactly P0 polynomials, P1 polynomials and P3 polynomials and approximate the integral with order 1, 2 and 3.

For each method, a model function (named rectangle(), trapz() or simpson()) is provided, enabling the user to supply a function with or without parameters, or the values of a function on a uniform set of points. Users have to specify the function, the interval bounds \(a,b\) and the number \(n\) of intervals. When a vector of function values is given, only the step \(h\) is required. For instance, using the trapeze method:

Real f(Real t) { return t*t; }
Complex fp(Real t, Parameters& pars)
{
  Int p=pars("p");
  return pow(t, p)*exp(i*t);
}
...
Number n=100;
Real intf = trapz(f, 0, 1, n);          // int_[0,1]f(t)dt
Parameters pa(2,"p");
Complex intfp = trapz(fp, pa, 0, 1, n); // using parameters
Vector<Real> vf(n+1);
Real h=1./n;
for (Number k=0; k<=n; k++) { vf[k]=f(k*h); }
Real intfv = trapz(vf, h); // using values instead of function

For advanced users, there exist also functions using iterators.

Adaptative trapeze#

XLiFE++ provides also an adaptive trapeze method using the improved trapeze method (\(0.25(b-a)(f(a)+2f((a+b)/2)+f(b))\)) to get an estimator. The user has to give its function (with parameters if required), the integral bounds and a tolerance factor (\(10^{-6}\) if not given):

Real intf = adaptiveTrapz(f, 0, 1);     // int_[0,1]f(t)dt with 1E-6 as tolerance
Parameters pa(2, "p");
Complex intfp = adaptiveTrapz(fp, pa, 0, 1, 1E-4); // with parameters and tolerance

Hint

Because at the first step, the adaptive method uses 5 points uniformly distributed, the result may be surprising if unfortunately the function takes the same value at these points!

Discrete fast Fourier transform (FFT)#

XLiFE++ provides the standard 1D discrete fast Fourier transform with \(2^n\) values:

\[g_k = \sum_{j=0}^{n-1} f_j e^{-2 i \pi \dfrac{jk}{n}}.\]

FFT and inverse FFT are computed using some functions addressing real or complex vectors of size \(2^n\):

Number ln=6, n=64;
Vector<Complex> f(n), g(n), f2(n);
for (Number i=0; i<n; i++) f[i]=std::sin(2*i*pi_/(n-1));
fft(f, g); ifft(g, f2); // f2 should be close to f
// or
Vector<Complex> fhat=fft(f), fihat=ifft(f);

Hint

If the vector you give is not of \(2^n\) length, fft/ifft computation will use only the first \(2^n\) values with \(2^n\) the greater value less than the vector size. Remember that even the input vector is a real vector the output is always a complex vector.

Advanced users can also use some functions fft/ifft addressing iterators that can handle any collection of real or complex values.

Filon’s method for oscillatory integrals#

Up to now, only one method is provided to compute 1D oscillatory integrals based on the Filon’s method. It is designed to compute an approximation of the following oscillatory integral:

\[I(x)=\int_0^T f(t) e^{-i x t}\]

where \(f\) is a slowly varying function with real or complex values. See FilonIM class for more details.

The FilonIM class provides mainly two constructors that pre-compute the values of \(f\) and \(f'\), a function that computes the integral at \(x\) (maybe a complex value). In constructors, ord_ corresponds to the interpolation used : 0 for Lagrange P0, 1 for Lagrange P1 and 2 for Hermite P3, f and df (if Hermite Filon is used) are some standard C++ functions of a real argument:

FilonIM(Number ord, T(f)(Real), Real tf, Number N);
FilonIM(Number ord, T(f)(Real), T(df)(Real), Real tf, Number N);
Complex compute(S& x) const;    // compute I(x)
Complex operator()(S& x) const; // compute I(x)
void print(ostream& os) const;  // print FilonIM on stream

This class is quite easy to use, for instance to compute the integral

\[I(x)=\int_0^{10} e^{-t} e^{-i x t}\]

by the P1-Filon method, define first a C++ function returning the value of \(f\) as a Complex, create a FilonIM object with your own parameters and finally call it passing an \(x\) value:

Complex f(Real x) { return Complex(exp(-x),0.); }
...
FilonIM fim(0, f, 10, 100);
Complex y = fim(0.5);

To invoke Hermite P3 Filon method, define also the derivative of \(f\):

Complex f(Real x) { return Complex(exp(-x),0.); }
Complex df(Real x) { return Complex(-exp(-x),0.); }
...
FilonIM fim(2, f, df, 10, 100);
Complex y = fim(0.5);

Hint

The P0 Filon method is of order 1, the P1 Filon is of order 2 while the Hermite P3 Filon is of order 4. The Hermite P3 Filon is less stable than the Lagrange Filon methods!

By default, FilonIM class is designed for complex value functions. In fact, this class is an alias of a template class that supports also real value functions. To explicitly use the real version of FilonIM class, you have to instantiate a FilonIMT<Real> object and pass to it some real value functions!

ODE solvers#

XLiFE++ provides some ODE solvers : Euler, Runge-Kutta 4 and Runge-Kutta 45 that is an adaptive step time solver based on the Dortmand-Prince method (see http://en.wikipedia.org/wiki/Dormand-Prince_method). All solvers are templated by the type of the state \(y\in V\) involved in the first order differential equation:

\[\begin{split}\left\{ \begin{array}{l} y'(t) = f(t,y(t)) \;\;\; t \in ]t_0; t_f[\\ y(t_0)=y_0\\ \end{array} \right.\end{split}\]

with \(f:]t_0; t_f[\times V \rightarrow V\) a non stiff function.

The user function f may be any C++ function of the form T& f(Real t, const T& y, T& fty) (T any scalar/vector/matrix type):

Real f(Real t, const Real& y, Real& fty){...; return fty;} // real scalar ODE
Complex& f(Real t, const Complex& y, Complex& fty)         // complex scalar ODE
Reals& f(Real t, const Reals& y, Reals& fty)               // real vector ODE
Complexes& f(Real t, const Complexes& y, Complexes& fty)   // complex vector ODE

ODE solvers are objects of template classes EulerT, RK4T and Ode45T but they can also be invoked by simple functions

template <typename T>
T& (*fun)(Real, const Real& y, Real& fty); // c++ function alias
Vector<T> euler(fun f, Real a, Real b, Real dt, const T& y0);
Vector<T> rk4(fun f, Real a, Real b, Real dt, const T& y0);
pair<Vector<Real>,Vector<T>> ode45(fun f, Real a, Real b, Real dt, const T& y0, [Real prec]);

In adaptive RK45 method, the time step may grow to becomes very large. To still get a meaningful number of times step, the time step is limited by the initial guess dt. Even if the initial guess dt is very large, the method should be convergent if the problem is not stiff.

The following example deals with the linear pendulum (ODE of order 2):

Real omg=1., th0=pi_/4, thp0=0.;
Reals& f(Real t, const Reals& y, Reals& fyt)
{
  fyt.resize(2);
  fyt[0]=y[1]; fyt[1]=-omg*omg*y[0];
  return fyt;
}
Reals yex(Real t)
{ //exact solution
  Reals yy(2,0.);
  yex[0] = thp0*sin(omg*t)/omg+ th0*cos(omg*t);
  yex[1] = thp0*cos(omg*t)- th0*omg*sin(omg*t);
  return yex;
}
...
Real t0=0,tf=1,dt=0.01, t=t0, errsup=0;
Reals y0(2,0.);y0[0]=th0; y0[1]=thp0;
Vector<Reals> sol = euler(f,t0,tf,dt,y0);
for(Number i=0;i<sol.size(); i++, t+=dt)
  errsup=max(errsup,norm(sol[i]-yex(t)));
theCout<<"euler : nb dt="<<sol.size()<<" error sup = "<<errsup<<eol;

t=t0; errsup=0;
sol = rk4(f,t0,tf,dt,y0);
for(Number i=0;i<sol.size(); i++, t+=dt)
errsup=max(errsup,norm(sol[i]-yex(t)));
theCout<<"rk4 : nb dt="<<sol.size()<<" error sup = "<<errsup<<eol;

t=t0; errsup=0;
pair<Reals,Vector<Reals>> solp = ode45(f,t0,tf,0.1,y0,1E-6);
Number n=solp.first.size();
for(Number i=0;i<n; i++)
errsup=std::max(errsup,norm(solp.second[i]-yex(solp.first[i])));
theCout<<"ode45 : nb dt="<<n<<" error sup = "<<errsup<<eol;

Because, ode45 is an adaptive method, the time step is not constant. So the ode45 function returns times and state values as a pair of vectors. The previous code gives

\[\begin{split}\begin{array}{lll} \mathrm{euler} &\ : \mathrm{nb\ dt} =101 &\mathrm{sup\ error} = 0.00393671 \\ \mathrm{rk4} &\ : \mathrm{nb\ dt} =101 &\mathrm{sup\ error} = 6.54499e-11 \\ \mathrm{ode45} &\ : \mathrm{nb\ dt} =12 &\mathrm{sup\ error} = 5.25629e-08 \end{array}\end{split}\]

showing the advantage to use Ode45!

Tip

Contrary to Euler and RK4 method, Ode45 handles other parameters :

  • nbtry_ the maximum number of tries (12) when adapting the time step,

  • minscale_ (0.125) and maxscale_ (4) to limit the variation of the time step at each iteration.

To modify these internal parameters, use the extended function

Number nbt = 20;
Real mins = 0.05, maxs = 2;
pair<Reals,Vector<Reals>> solp = ode45(f,t0,tf,0.1,y0,nbt,mins,maxs,1E-6);
On this page