For detailed description see capd::map::Function and capd::map::Map
These classes provide methods for computing:
- values of functions/maps
- normalized derivatives (Taylor coefficients) of maps
- jet propagation through the map (in fact computation of a jet of composition of two maps)
Defined types and data structures
The main header file
#include "capd/capdlib.h"
defines the following types for computation in double (D), long double (LD) precision and interval (I) arithmetics.
- DFunction, LDFunction, IFunction - scalar valued functions

- DMap, LDMap, IMap - vector valued maps

- DHessian, LDHessian, IHessian - data structure to store second order normalized derivatives (Taylor coefficients) of maps
- DJet, LDJet, IJet - data structure to store multivariate polynomials, i.e. truncated Taylor series of maps
- Multiindex, Multipointer - data structures used to index jets.
Defining functions and maps
Maps can be parsed from a human readable string and/or C-routine. Performance of further computations does not depend on the way you define an object. This is for user convenience only - short functions can be defined as a string while large expressions are easier to encode as routines.
Parsing from a string. Syntax of the formula:
"[par:a,b,c...;][time:t;]var:x1,x2,...;fun:expression1,expression2,....;"
- var - names of subsequent arguments
- fun - expressions that define a map. You can use most elementary
- functions: sin, cos, exp, log, sqrt, sqr (square)
- operators: +,-,*,/,^ (power, integer or not. Exponent cannot depend on variables - we assume gradient of exponent is zero).
- constants: (0,1,2,-5,2.5,-0.25, etc)- we recommend usage of representable numbers only. If a constant is an interval or a floating point represented with high precision, one should use a parameter instead.
- par - parameters of the map. Derivatives of them with respect to main variables are assumed to be zero.
- time - a distinguished parameter with derivative dt/dt=1. Used to define time-dependent maps that stand for vector fields of nonautonomous ODEs.
- Attention
- The parser does not accept numerical constants in scientific notation. Use parameters, instead.
- Note
- Sections par and time are optional.
Example:
#include "capd/capdlib.h"
capd::IMap lorenz(
"par:s,r,q;var:x,y,z;fun:s*(y-x),x*(r-z)-y,x*y-q*z;");
lorenz.setParameter("s",10.);
lorenz.setParameter("r",28.);
Parsing from a C-routine. One has to write a global function (or a functor, lambda expression, can be class member function) that defines the map and then send it to the constructor of class Map.
The function that defines a map must have the following signature:
#include "capd/capdlib.h"
);
Here:
- t - is a distinguished time variable
- in - is a C-array of input (independent) variables
- dimIn - is number of input variables
- out - is a C-array of output (dependent) variables
- dimOut - is number of output variables
- params - is a C-array of parameters
- noParam - is number of parameters
See section Computing values and derivatives of maps for a complete example.
Computing values and derivatives of maps
Class Map provides methods for computation of values and derivatives. This is simply done by call to overloaded operator
() and/or method derivative
.
Complete example (from examples/maps/mapExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace std;
void pcr3bpVectorField(Node , Node in[], int , Node out[], int , Node params[], int )
{
Node mj = 1-params[0];
Node xMu = in[0] + params[0];
Node xMj = in[0] - mj;
Node xMuSquare = xMu^2;
Node xMjSquare = xMj^2;
Node ySquare = in[1]^2;
Node factor1 = mj*((xMuSquare+ySquare)^-1.5);
Node factor2 = params[0]*((xMjSquare+ySquare)^-1.5);
out[0] = in[2];
out[1] = in[3];
out[2] = in[0] - xMu*factor1 - xMj*factor2 + 2*in[3];
out[3] = in[1]*(1 - factor1 - factor2) - 2*in[2];
}
cout.precision(21);
int dimIn=4, dimOut=4, noParam=1;
LDMap f(pcr3bpVectorField,dimIn,dimOut,noParam);
f.setParameter(0, 0.125);
long double v[] = {2,3,4,5};
cout << "f(x)=" << y << endl;
cout << "Df(x)=" << Df << endl;
cout << "y-y2=" << y-y2 << endl;
cout << "Df-Df2=" << Df-Df2 << endl;
}
- Attention
- For rigorous computation you must use in the expression representable numbers only (both in the string parsed and routine parsed cases). If you want to use constants, like 1/3 then you can
- either set them as parameters of the system (strongly recommended) - do not forget to set the parameter value after creating an instance of Map.
- or wrap them into Node type (slower): String parser wraps automatically expressions like 1/3 to division of nodes.
Hessians of maps
Class Map provides algorithms for computation of normalized hessians of maps. One has to
- set maximal order of derivative to at least 2, when creating an instance of Map
DMap f(globalFunction,dimIn,dimOut,noParams,maxDerivative );
- define objects that will store computed derivative and hessian
- call operator
()
Now Hf
stores second order normalized derivatives. They can be accessed by
int fi =..., dxj = ..., dxk = ...;
double dfi_dxj_dxk = Hf(fi,dxj,dxk);
- Note
- Indexing of functions (
fi
) and derivatives (dxj,dxk
) starts from zero.
Complete example (from examples/maps/hessianExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace std;
void _f(Node , Node in[], int , Node out[], int , Node params[], int noParams)
{
out[0] = in[0]+in[1];
out[1] = -(in[2]+in[3]);
for(int i=0;i<noParams;++i){
Node temp = params[i]*
sqrt(
sqr(out[0] - in[0]) +
sqr(out[1] - in[3]));
out[1] = params[i]*
sqrt(
sqr(out[0] + in[1]) +
sqr(out[1] - in[2]));
out[0] = temp;
}
}
int dimIn=4, dimOut=2, noParam=5;
int maxDerivativeOrder = 2;
IMap f(_f,dimIn,dimOut,noParam,maxDerivativeOrder);
for(int i=0;i<noParam;++i)
cout.precision(17);
cout << "y=" << y << endl;
cout << "Df=" << Df << endl;
for(int fi=0;fi<dimOut;++fi)
for(int dx1=0;dx1<dimIn;++dx1)
for(int dx2=dx1;dx2<dimIn;++dx2)
cout << "Hf(" << fi << "," << dx1 << "," << dx2 << ")=" << Hf(fi,dx1,dx2) << endl;
}
Higher order Taylor coefficients of maps
Class Map provides algorithms for computation of normalized higher order derivatives of maps. One has to
- set maximal order of derivative to desired value, when creating an instance of Map
DMap f(globalFunction,dimIn,dimOut,maxRequestedDerivative );
- define objects that will store jets of maps. Class [L]DJet (and interval version IJet) represents truncated Taylor series of a map.
DJet fJet(dimOut,dimIn,truncationDegreeInclusive );
- call operator
()
Now fJet
stores value of f(x)
and all normalized derivatives of f
at x
up to order truncationDegreeInclusive
. For partial derivatives of order less or equal than three we provide direct access to coefficients by simple operator call
int fi =..., dxj = ..., dxk = ..., dxc = ...;
double dfi_dxj = fJet(fi,dxj);
double dfi_dxj_dxk = fJet(fi,dxj,dxk);
double dfi_dxj_dxk_dxc = fJet(fi,dxj,dxk,dxc);
- Note
- Indexing of functions (
fi
) and derivatives (dxj,dxk
) starts from zero.
- Attention
- For performance reasons indices in
fJet(fi,dxj,dxk,dxc)
must form a nondecreasing sequence! Otherwise behaviour is undefined.
Complete example (from examples/maps/jetExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace std;
void _f(Node , Node in[], int , Node out[], int , Node params[], int noParams)
{
out[0] = in[0]+in[1];
out[1] = -(in[1]+in[2]);
for(int i=0;i<noParams;++i){
Node temp = params[i]*
sqrt(
sqr(out[0] - in[1]) +
sqr(out[1] - in[2]));
out[1] = params[i]*
sqrt(
sqr(out[0] + in[2]) +
sqr(out[1] - in[0]));
out[0] = temp;
}
}
int dimIn=3, dimOut=2, noParam=5;
int maxDerivativeOrder = 3;
DMap f(_f,dimIn,dimOut,noParam,maxDerivativeOrder);
for(int i=0;i<noParam;++i)
f.setParameter(i, double(i+1)/9.);
double v[] = {2,3,4};
DJet jet(dimOut,dimIn,maxDerivativeOrder);
f(x,jet);
cout << "df1_dx0_dx0_dx2 = " << jet(1,0,0,2) << endl;
cout << jet.toString();
}
Indexing of higher order derivatives
Fourth and higher order derivatives are indexed by multipointers and multiindices.
Multiindex: is an element of
(sequence of integers) of the same length as the dimension of the domain of function, i.e.
.
i-th coefficient of a multiindex indicates the order of partial derivative with respect to i-th variable. For instance,
corresponds to normalized derivative 
Multipointer: is a nondecreasing sequence of integers of the length d, where d is the total order of partial derivative. Every element of multipointer is an index of variable with respect to which we take partial derivative. For instance, the following Multipointer
int data[] = {0,0,1,3,3};
is equivalent to Multiindex {2,1,0,2} and it corresponds to normalized derivative 
Given multipointer or multiindex one can access corresponding coefficient in the data structure Jet. One can also easily convert between multipointers and multiindices by constructor calls.
Complete example (from examples/maps/jetIndexingExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace std;
void _f(Node , Node in[], int , Node out[], int , Node params[], int noParams)
{
out[0] = in[0]+in[1];
out[1] = -(in[1]+in[2]);
for(int i=0;i<noParams;++i){
Node temp = params[i]*
sqrt(
sqr(out[0] - in[1]) +
sqr(out[1] - in[2]));
out[1] = params[i]*
sqrt(
sqr(out[0] + in[2]) +
sqr(out[1] - in[0]));
out[0] = temp;
}
}
int dimIn=3, dimOut=2, noParam=5;
int maxDerivativeOrder = 4;
DMap f(_f,dimIn,dimOut,noParam,maxDerivativeOrder);
for(int i=0;i<noParam;++i)
f.setParameter(i, double(i+1)/9.);
double v[] = {2,3,4};
DJet jet(dimOut,dimIn,maxDerivativeOrder);
f(x,jet);
int data[] = {0,2,2};
cout << "multiindex: " << m << endl;
cout << "vector of Taylor coefficients D^{m}f = " << jet(m) << endl;
cout << "selected Taylor coefficient D^{m}f_1 = " << jet(1,m) << endl << endl;
cout << "conversion to multipointer: " << mp << endl;
cout << "vector of Taylor coefficients D^(mp}f = " << jet(m) << endl;
cout << "selected Taylor coefficient D^(mp}f_1 = " << jet(1,m) << endl << endl;
do{
cout << p << ": " << jet(p) << endl;
}while(jet.hasNext(p));
return 0;
}
Jet transport - computing composition of jets
Class Map defines an operator for computating jet of a composition of the map and other jet sent as an argument. Given map f
and jet g at a point x
one can compute jet of f(g(x)) by simply call
Complete example (from examples/maps/jetTransportExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace std;
void _f(Node , Node in[], int , Node out[], int , Node [], int)
{
Node d =
sqr(in[0])+
sqr(in[1]);
out[0] = (
sqr(in[0])+in[1])/d;
out[1] = (
sqr(in[1])+in[0])/d;
}
int dimIn=2, dimOut=2, noParam=0;
int maxDerivativeOrder = 3;
DMap f(_f,dimIn,dimOut,noParam,maxDerivativeOrder);
double v[] = {2,3};
DJet fx(dimOut,dimIn,maxDerivativeOrder);
f(x,fx);
cout << "f(x):\n" << fx.toString() << endl;
return 0;
}