Manual of scientific software for validated (rigorous, interval) numerics for partial
differential equations.
Software constructing reachability bounds for the solution set of parabolic PDEs.
Public repository with the most recent version
Introduction
The goal of this document is to provide a manual on the usage of the software package for validated integration of partial differential equations.
The installation guide can be found here.
The purpose of the software is to provide the research community with a tool (solver) for the rigorous integration of dPDEs (dissipative partial differential equations, a subclass of PDEs) with periodic boundary conditions.
By rigorous integration I mean a process in which a compact and connected set containing the unique solution is propagated along the time interval. The rigorous integration process guarantees to include the exact solution in the output set. Whereas in any standard (non-rigorous) integration process a point, represented by floating point numbers, is being propagated without indication how close is the approximated solution to the real one.
Fig. Difference between a regular non-rigorous method with a rigorous method.
Observe that in the rigorous method setting, the interval boxes, which are being propagated in time, can be of any finite dimension. In case of partial differential equations those boxes are in fact infinite dimensional. Before an algorithm dealing with infinite dimensional sets has been worked out, a suitable theory had been established. Refer to the work in References for more details.
One of the purposes of the software package is to permit performing computer assisted proofs of existence of certain solutions, including periodic orbits, connections between equilibria, existence of attractors (including globally attracting ones).
To the best of my knowledge this is the first and the only software package for validated numerical integration of partial differential equations available at the moment.
There exists several packages capable of performing the described task, but for ordinary differential equations (ODEs) only. The case of PDEs is by a great deal different from ODEs, as PDEs can be viewed as infinite dimensional ODEs. Let me mention the CAPD library, and the VNODE library. Here can be found the list of example results , computer assisted proofs that have been obtained by using the CAPD rigorous ODE solver.
Current development
Currently I am working on new features, including adapting the software for the bifurcation analysis of 2D Navier-Stokes related systems. Once this is done I am going to introduce the computer assisted proofs package for 2D Navier-Stokes related systems.Boundary conditions and type of problems that could be considered.
Currently, the software allow to deal with scalar one-dimensional equations with periodic boundary condition. Additionally it is required that in the given PDE the elliptic operator dominates the nonlinear operator in some sense.
Eigenfunctions of the Dirichlet ($\varphi_k(x)=\sin{kx},\ k\in\mathbb{N}$) and Neumann problem ($\varphi_k(x)=\cos{kx},\ k\in\mathbb{N}$) for the Laplacian on an interval can also be used.
``Solve'' partial differential equation included in package
In the package found here there are included a few examples, in the source files a partial differential equation is defined, and then is rigorously integrated until a condition is fulfilled. The examples are as followsBurgers equation
\[ \begin{aligned} &u_t+u\cdot u_x-\nu u_{xx}=0,\\ &u(x,0)=u_0(x),\quad x\in\mathbb{R},\\ &u(x,t)=u(x+2k\pi, t),\quad x\in\mathbb{R},\ t\in[0, \mathcal{T}),\ k\in\mathbb{Z}. \end{aligned} \] An example how to integrate an initial condition for this equation is found in the source fileBurgersFPTest.cpp
Kuramoto-Sivashinsky equation
\[ \begin{aligned} &u_t=-\nu u_{xxxx}-u_{xx}+(u^2)_x,\\ &u(x,0)=u_0(x),\quad x\in\mathbb{R},\\ &u(x,t)=-u(-x,t),\quad u(x,t)=u(x+2k\pi, t),\quad x\in\mathbb{R},\ t\in[0, \mathcal{T}),\ k\in\mathbb{Z}. \end{aligned} \] An example how to integrate an initial condition for this equation is found in the source filesKSPOTest.cpp,
KSUnstPOTest.cpp.
In the provided examples a box is propagated along some periodic orbits taken from the paper P. Zgliczyński, Rigorous Numerics for Dissipative PDEs III. An effective algorithm for rigorous integration of dissipative PDEs, Topological Methods in Nonlinear Analysis, 36 (2010) 197--262. Initial conditions are some boxes on one stable periodic orbit and one unstable periodic orbit.
Swift-Hohenberg equation
\[ \begin{aligned} &u_t=\left(\nu-\left(1+\frac{\partial^2}{\partial x^2}\right)^2\right)u-u^3,\label{eq:SH}\\ &u(x,0)=u_0(x),\quad x\in\mathbb{R},\\ &u(x,t)=u(-x,t),\quad u(x,t)=u(x+2k\pi, t),\quad x\in\mathbb{R},\ t\in[0, \mathcal{T}),\ k\in\mathbb{Z}. \end{aligned} \] An example how to integrate an initial condition for this equation is found in the source fileSHTest.cpp
In the provided example a box close to the unstable fixed point is being propagated until it is attracted by an attracting fixed point.
Fig. Setting for Swift-Hohenberg equation.
Cahn-Hillard equation and Diblock Copolymer model
\[ \begin{aligned} &u_t=-(u_{xx}+\lambda f(u))_{xx}-\lambda\sigma(u-\mu),\label{eq:CH}\\ &u(x,0)=u_0(x),\quad x\in\mathbb{R},\\ &\mu=\int_0^1{u(x)}\,dx\\ &u_x=u_{xxx}=0,\quad x=0,\ 2\pi. \end{aligned} \]``Solve'' a custom partial differential equation
As for now, the software can be easily adapted to any equation which has the nonlinearity given by the following polynomial \[ \begin{aligned} \frac{\partial^i u^p}{\partial x^i}, \end{aligned} \] i.e. $i$-th derivative of $u^p$. , for instance \[ \begin{aligned} &(u^2),\\ &(u^2)_x=2u\cdot u_x,\\ &(u^3),\\ &(u^3)_{x}=3u^2 u_x,\\ &(u^3)_{xx}=3(2u\cdot u_x^2+u^2\cdot u_{xx}),\\ &\text{etc...} \end{aligned} \]
Here I present a case study illustrating how to define a ``custom'' partial differential equation and how to integrate this equation in time.
- First a class defining the given equation inheriting from
capd::jaco::DPDEDefinitionshould be provided. The example below presents what functions should be implemented on the example of the Burgers equation - Typedefs. Then, in the main source .cpp file the types should be defined regarding the equation and the boundary conditions that are going to be used. The example shows the case of the periodic-odd bd. cond.. The important definitions are presented in the following piece of code
-
As ``scalar'' values for the solver we use the so-called first order jets in order to obtain the partial derivatives with respect to the initial condition. Those partial derivatives are needed to obtain a suitable coordinate frame for the Lohner representation of a set. First order jets are first degree homogeneous multivariate polynomials alike. An algebra of first order jets can be defined -- the higher order terms resulting from multiplication are simply being truncated.
By using this approach the variational equations does not have to be provided separately at all. The same solver is working for the any type of periodic boundary conditions (cosine, sine or both), moreover automatic optimizations are performed to avoid superfluous computations, some information can be found in my PhD dissertation. And the vector field that is going to be integrated does not have to be redefined whenever the type of periodic boundary conditions are changed.
Those first order jets has to be initialized with the desired boundary conditions and the dimension. This is done in the following piece of code.
//begin FOJ initialization for FFT integrator capd::jaco::DPDEContainer container; //2.set here the subspace of the initial condition e.g. setToRealValuedOdd means that the initial condition is real valued odd container.setToRealValuedEven(); FOJ1D::initialize(ModesContainer1D::modes2arraySizeStatic(m), container); //end FOJ initialization -
Then the actual integrator is constructed
DPDEInclusionCW diffIncl(m, dftPts, M, dftPts2, PI, nu, order, step, MaxNorm());with the following parametersm- the Galerkin projection dimensiondftPts- the number of discrete points used by FFT (for the projection of size $m$), should satisfy $> 2m+1$, there are some subtle issues related to the aliasing,M- the Tail dimension (the modes out of a Galerkin projection are denoted by ``Tail''),dftPts2- the number of discrete points used by FFT (for the projection of size $M$), should satisfy $> 2M+1$, there are some subtle issues related to the aliasing,PI- interval enclosing the value of $\pi$,nu- the parameter appearing in the elliptic operator (for some equations viscosity),order- order of the Taylor method used,step- step size (constant),MaxNorm()- the norm used in calculating the logarithmic norm.
-
The initial condition should be provided in the
RealPolynomialBoundclass. The values in the finite part are set by using theoperator[]. The polynomial bounds for the infinite part $C>0$ and $s>0$ such that $|a_k|\leq\frac{C}{|k|^s}$ are being set by thesetCandsetSfunctionsRealPolynomialBound u_0(m, M, container); setC(u_0, 0); setS(u_0, 6);Then the set is constructed provided the initial conditionInclRect2Set set(u_0);. -
The set is propagated by calling (repeatedly)
set.move(diffIncl);
template<class PolyBdT> class Burgers : public PolyBdT::SubspaceType, public capd::jaco::DPDEDefinition{ public: RealType nu; int m_p; int m_d; int m_r; int m_sufficientlyLarge; ComplexType m_N_coeff; static const double S_DISSIPATIVE = -1; Burgers(RealType nu_) : nu(nu_){ m_p=2.; //partial derivative order in the elliptic operator m_d=1.; //spatial dimension (currently only 1 is supported) m_r=1.; //derivative order in the nonlinear operator (will work only for the second degree polynomial in the nonlinear part) m_sufficientlyLarge=m_d+m_p+1; //minimal regularity of solutions for which we prove existence and uniqueness, leave at it is m_N_coeff = -0.5; //constant coefficient in front of the polynomial in the nonlinear part } RealType ni(int k) const{ return nu; } RealType ni(const IndexType& index) const{ return nu; } const ComplexType& Ncoeff() const{ return m_N_coeff; } //leave this function as it is, the dissipative directions threshold should be adjusted bool isDissipative(int k){ if(lambda(k) < S_DISSIPATIVE){ return true; } return false; } // i is position in an array of a mode. Returns real part of the eigenvalue \lambda RealType lambda(int i) const{ IndexType j=array2modeIndex(i); return -j.squareEuclNorm() * ni(j); } // k is the index of mode, lambda_k from the paper. RealType lambda_k(int k) const{ return -k * k * ni(k); } RealType lambda_k(const IndexType& k) const{ return -k.squareEuclNorm() * ni(k);} // function returning V(K)=\{ \inf{v(|k|} | |k|>=K} \}, see definition of dissipative PDE in the paper. // Eigenvalues of a dPDE satisfies \lambda_k=-\nu(|k|)|k|^p . RealType V(int K) const{ return leftBound(ni(K)); } RealType V(const IndexType& k) const{ return leftBound(ni(k)); } // returns smallest integer K larger than k, such that f(i) is monotonously non increasing function for i>K. // Where $f(i) = e^{h\lambda_k} k^r$ int maximumPoint(const RealType& h, int r, int k) const{ return ceil(sqrt(RealType(r/(2.*nu*h)).rightBound())); } };
typedef capd::jaco::KS<RealPolynomialBound> KS; //the class defining the equation
typedef capd::jaco::Odd<Interval, Index1D, Norm1D, IntervalMatrix> OddSubspace; //the class defining the boundary conditions, here there can be used two other possible classes: Even, Real
typedef capd::jaco::RealPolynomialBound<ComplexScalar, Index1D, 0, OddSubspace> RealPolynomialBound; //the class defining the polynomial bounds - finite representation of infinite dimensional sets
typedef capd::jaco::FFT1D FFT1D;//the class implementing the Fast Fourier Transform algorithm
//{
typedef capd::jaco::DPDE2 DPDE; //when the nonlinear operator is a second degree polynomial use this definition,
typedef capd::jaco::DPDE3 DPDE; //whereas when the nonlinear operator is a third degree polynomial use this definition
//}
typedef capd::jaco::FirstOrderJet FOJ1D; //definition of complex first order jets, this is extended ``scalar'' type allowing to obtain the Jacobian matrix
typedef capd::jaco::ComplexPolyBdJetOptimized ModesContainer;
typedef capd::jaco::FFT1D JetFFT1D;
//{
typedef capd::jaco::DPDE2 JetDPDE; //when the nonlinear operator is a second degree polynomial use this definition,
typedef capd::jaco::DPDE3 JetDPDE; //whereas when the nonlinear operator is a third degree polynomial use this definition
//}
typedef capd::jaco::FFTTaylorDynSys FFTDynSys;//the integrator for integration of the finite part of the equations (a Galerkin projection)
typedef FFTDynSys::ModesContainerType ModesContainer1D;
typedef capd::jaco::DPDEInclusionCW DPDEInclusionCW;//differential inclusion, solve infinite dimensional system of equations
typedef capd::diffIncl2::InclRect2Set InclRect2Set;//this is the Lohner representation of sets for the differential inclusion
All the appearing types Interval, Index1D, Norm1D, IntervalMatrix, ComplexScalar, ComplexDerivativePair should be defined like in the source
files with examples.
Included FFT library
As a byproduct of this research I have developed a nice template FFT (Fast Fourier Transform) library in C++. This is included in the software package. The FFT portion of the library can be used with intervals or any class with overloaded basic algebraic operations (addition and multiplication).
Implementation, dependence on external libraries, and possible use of tools (Poincare map, covering relations) from CAPD library
The package has been written in C++ programming language. It makes use of the C++ code templates to achieve generality. The only external dependencies of the library are the FILIB interval arithmetic library and the CAPD library the modules for the Lohner representation of sets, and solving differential inclusions are used. All required source codes are included in the software package.
The code is CAPD library compatible. Meaning that the software package can be seen as a module to the CAPD library. For instance the solver from the package inherits from the CAPD solver class to allow use of the other CAPD tools like Poincare map or covering relations (but this feature has not been tested yet).
Concluding remarks
In case of any questions feel free to contact me via jcyranka at g mail
Github repository https://github.com/dzako/dpdes