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.
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.
BurgersFPTest.cpp
KSPOTest.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.
SHTest.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.
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.
capd::jaco::DPDEDefinition
should be provided.
The example below presents what functions should be implemented on the example of the Burgers equationtemplate<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::FFT1DAll the appearing typesFFT1D;//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
Interval, Index1D, Norm1D, IntervalMatrix, ComplexScalar, ComplexDerivativePair
should be defined like in the source
files with examples.
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
DPDEInclusionCW diffIncl(m, dftPts, M, dftPts2, PI, nu, order, step, MaxNorm());with the following parameters
m
- 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.RealPolynomialBound
class. The values in the finite part are set by
using the operator[]
. 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 the setC
and setS
functions
RealPolynomialBound u_0(m, M, container); setC(u_0, 0); setS(u_0, 6);Then the set is constructed provided the initial condition
InclRect2Set set(u_0);
.
set.move(diffIncl);
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).
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).
In case of any questions feel free to contact me via .
There is a bitbucket.org private repository with all the source codes. If you need/would like to access email me, please.
The licence of the software package reads:
dPDEs - this program is an open research software performing rigorous integration in time of partial differential equations Copyright (C) 2010-2013 Jacek Cyranka This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.