Integration and ODE#

Integral calculation tools#

To compute finite elements integrals by quadrature, XLiFE++ provides the following routines in order to compute standard quadrature points and weights on interval $[- 1, 1]$ :

the Gauss-Legendre formula (gaussLegendreRule)
the Gauss-Lobatto formula (gaussLobattoRule)
the Gauss-Jacobi formula (gaussJacobiRule)

The Gauss-Legendre and Gauss-Lobatto routines only return the $(n + 1) / 2$ first positive points in ascending order whereas Gauss-Jacobi routine returns $n$ points (negative and positive) in ascending order. Note that for low order, points and weights are tabulated (long double precision):

Reals abcs, weis;
Number n = 10;
gaussLegendreRule(n, abcs, weis);  // tabulated up to ordrer 16, then computed
gaussLobattoRule (n, abcs, weis);  // tabulated up to ordrer 16, then computed
gaussJacobiRule  (n, abcs, weis);  // tabulated up to 10 points

Uniform rectangle, trapeze and Simpson methods#

XLiFE++ provides some general tools to approximate 1D integrals:

\int_{a}^{b} f (t) d t .

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

rectangles method:: $\int_{a}^{b} f (t) d t \sim h \sum_{i = 0}^{n - 1} f (a + i h)$
trapezes method:: $\int_{a}^{b} f (t) d t \sim \frac{h}{2} (f (a) + 4 \sum_{i = 1}^{n / 2} f (a + (2 i - 1) h) + 2 \sum_{i = 1}^{n / 2 - 1} f (a + 2 i h) + f (b)) (n even)$
Simpson method:: $\int_{a}^{b} f (t) d t \sim \frac{h}{3} (f (a) + 2 \sum_{i = 1}^{n - 1} f (a + i h) + f (b))$

with $f$ a real or a 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) + 2 f ((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}$ by default):

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 π \frac{j k}{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 given vector 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.

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{array}{r} {\begin{cases} y^{'} (t) = f (t, y (t)) t \in] t_{0}; t_{f} [ \\ y (t_{0}) = y_{0} \end{cases} \end{array}

with $f :] t_{0}; t_{f} [\times V \to 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{array}{r} \begin{array}{lll} euler & : nb dt = 101 & \sup error = 0.00393671 \\ rk 4 & : nb dt = 101 & \sup error = 6.54499 e - 11 \\ ode 45 & : nb dt = 12 & \sup error = 5.25629 e - 08 \end{array} \end{array}

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);