ODEs - functional approach
Table of Contents
Solutions to IVPs as functions
Class [L]DTimeMap provides methods for computing solutions to IVPs that behave like regular functions. This mechanism is written in the spirit of functional programming. Using this method you will get a function, say solution, that represents solution to IVP and you will be able to evaluate this function at any intermediate time
DVector u = solution(t);
without additional integration of ODE.
This method is recommended if you
- do not know exact time
tat which you will need the solution or - you need solution values for many intermediate times
t
In order to obtain solution curve you have to
- define an instance of SolutionCurve - this is type defined in class [L]DTimeMap double initialTime = ....;DTimeMap::SolutionCurve solution(initialTime);
- specify initial condition and call your timeMap integrator After the last line the object solution becomes a function that you can evaluate at given time. We will describe details in the sequel.
The type DTimeMap::SolutionCurve provides (in particular) methods for
- obtaining the domain at which this curve is defined - this should be interval [initialTime,finalTime] double L = solution.getLeftDomain();double R = solution.getRightDomain();
- evaluating the function. For
![$ t\in[L,R]$](form_87.png)
you can compute which is an approximate solution to IVP at timeDVector u = solution(t);t.- Note
- An exception is thrown when the argument
tin call tosolution(t)is out of the domain
Complete example (from examples/odes/DTimeMapSolutionCurveExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace capd;
using namespace std;
int main(){
cout.precision(21);
// define an instance of class DMap that describes the vector field
LDMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
pendulum.setParameter("omega",1.);
int order=20;
LDOdeSolver solver(pendulum,order);
LDTimeMap timeMap(solver);
// specify initial condition
LDVector u(2);
u[0] = 1.;
u[1] = 2.;
long double initTime = 4;
long double finalTime = 10;
// define functional object
LDTimeMap::SolutionCurve solution(initTime);
// and integrate
timeMap(finalTime,u,solution);
// then you can ask for any intermediate point on the trajectory
cout << "domain = [" << solution.getLeftDomain() << "," << solution.getRightDomain() << "]\n";
cout << "solution(4) =" << solution(4) << endl;
cout << "solution(7) =" << solution(7) << endl;
cout << "solution(10)=" << solution(10) << endl;
return 0;
}
/* output:
domain = [4,10]
solution(4) ={1,2}
solution(7) ={-0.451828248033488868443,-1.67522823683216929268}
solution(10)={1.20483703334825757551,1.17041482661774720001}
*/
First order variational equations - functional approach
As in the case of solving IVPs (see Solutions to IVPs as functions) one can obtain solution to variational equation as a functional object that can be evaluated at any intermediate time. In principle there are two differences:
- We can specify initial condition for variational equations. This argument (initMatrix in the example below) is optional. If it is skipped then the Identity matrix will be chosen as an initial condition for variational equations.
- We have to send an additional argument to operator that receives monodromy matrix at the end of trajectory and indicates that we intend to integrate variational equations.
Complete example (from examples/odes/DTimeMapMonodromyMatrixCurveExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace capd;
using namespace std;
int main(){
// define an instance of class DMap that describes the vector field
LDMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
pendulum.setParameter("omega",1.);
int order=20;
LDOdeSolver solver(pendulum,order);
LDTimeMap timeMap(solver);
// specify initial condition
LDVector u(2);
u[0] = 1.;
u[1] = 2.;
long double initTime = 4;
long double finalTime = 10;
long double data[] = {1,1,3,4};
LDMatrix initMatrix(2,2,data); // {{1,1},{3,4}}
LDMatrix monodromyMatrix(2,2);
// define functional object
LDTimeMap::SolutionCurve solution(initTime);
// and integrate
timeMap(finalTime,u,initMatrix,monodromyMatrix,solution);
// then you can ask for any intermediate point on the trajectory
cout.precision(21);
cout << "domain = [" << solution.getLeftDomain() << "," << solution.getRightDomain() << "]\n";
cout << "solution(4) =" << solution(4) << endl;
cout << "solution(7) =" << solution(7) << endl;
cout << "solution(10)=" << solution(10) << endl;
cout << "monodromyMatrix(4)=" << solution.derivative(4) << endl;
cout << "monodromyMatrix(7)=" << solution.derivative(7) << endl;
cout << "monodromyMatrix(10)=" << solution.derivative(10) << endl;
cout << monodromyMatrix << endl;
return 0;
}
/* output:
domain = [4,10]
solution(4) ={1,2}
solution(7) ={-0.451828248033488868443,-1.67522823683216929268}
solution(10)={1.20483703334825757551,1.17041482661774720001}
monodromyMatrix(4)={
{1,1},
{3,4}
}
monodromyMatrix(7)={
{10.1969455920725050385,13.1893571903452595699},
{-3.06072164458348302597,-3.8608572219154975219}
}
monodromyMatrix(10)={
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
{
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
*/
Second order variational equations - functional approach
As in the case of solving IVPs (see Solutions to IVPs as functions) one can obtain solution to first and second order variational equation as a functional object that can be evaluated at any intermediate time. In principle there are two differences:
- We can specify initial condition for first and second order variational equations. These arguments (initMatrix and initHessian in the example below) are optional. If both of them are skipped then the Identity matrix will be chosen for first order variational equations and zero for the second order variational equations.
- We have to send two additional arguments to the operator that computes trajectory. These arguments receive monodromy matrix and normalized hassian at the end of trajectory, respectively. Moreover, they indicate that we intend to integrate second order variational equations.
Complete example (from examples/odes/DTimeMapHessianCurveExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace capd;
using namespace std;
cout << "hessian(" << t << ")={";
cout << " " << h(i,j,c);
cout << " }\n";
}
int main(){
// define an instance of class DMap that describes the vector field
LDMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);");
pendulum.setParameter("omega",1.);
int order=20;
LDC2OdeSolver solver(pendulum,order);
LDC2TimeMap timeMap(solver);
// specify initial condition
LDVector u(2);
u[0] = 1.;
u[1] = 2.;
long double initTime = 4;
long double finalTime = 10;
long double initDerData[] = {1,1,3,4};
long double initHessData[] = {1,3,3,0,2,1}; // {dxdx,dxdy,dydy} = {{1,3,3},{0,2,1}}
LDMatrix initMatrix(2,2,initDerData); // {{1,1},{3,4}}
LDHessian initHessian(2,2,initHessData);
LDMatrix D(2,2);
LDHessian H(2);
// define functional object
LDC2TimeMap::SolutionCurve solution(initTime);
// and integrate
timeMap(finalTime,u,initMatrix,initHessian,D,H,solution);
// then you can ask for any intermediate point on the trajectory
cout.precision(21);
cout << "domain = [" << solution.getLeftDomain() << "," << solution.getRightDomain() << "]\n";
cout << "solution(4) =" << solution(4) << endl;
cout << "solution(7) =" << solution(7) << endl;
cout << "solution(10)=" << solution(10) << endl;
cout << "monodromyMatrix(4)=" << solution.derivative(4) << endl;
cout << "monodromyMatrix(7)=" << solution.derivative(7) << endl;
cout << "monodromyMatrix(10)=" << solution.derivative(10) << endl;
cout << D << endl;
print(solution.hessian(4),4);
print(solution.hessian(7),7);
print(solution.hessian(10),10);
print(H,10);
return 0;
}
/* output:
domain = [4,10]
solution(4) ={1,2}
solution(7) ={-0.451828248033488868443,-1.67522823683216929268}
solution(10)={1.20483703334825757551,1.17041482661774720001}
monodromyMatrix(4)={
{1,1},
{3,4}
}
monodromyMatrix(7)={
{10.1969455920725050385,13.1893571903452595699},
{-3.06072164458348302597,-3.8608572219154975219}
}
monodromyMatrix(10)={
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
{
{-7.80565667723112838715,-10.0010040066793847611},
{-2.61890487701602830734,-3.48358623656441465785}
}
hessian(4)={ 1 3 3 0 2 1 }
hessian(7)={ 42.6352224242584004026 114.667151261346096119 73.2754844459227481224 17.2114064486085741219 42.4899350405156879182 26.8990700123135739949 }
hessian(10)={ -29.1661902479786202484 -79.9855736889665852526 -52.1961182105122759693 -14.8470974997925887675 -38.9769249314583959565 -24.2494060090330983866 }
hessian(10)={ -29.1661902479786202484 -79.9855736889665852526 -52.1961182105122759693 -14.8470974997925887675 -38.9769249314583959565 -24.2494060090330983866 }
*/
Higher order variational equations - functional approach
The class [L]DCnTimeMap can compute solutions to higher order variational equations as a functional objects that can be evaluated at any intermediate time. In principle there are two differences:
- We can specify initial condition for all variational equations. This argument (initJet in the example below) is optional. If it is skipped then the Identity matrix will be chosen for first order variational equations and zero for higher order variational equations.
- We have to send one additional arguments to the operator that computes trajectory. This argument receives jet of solution and indicates that we intend to integrate higher order variational equations.
Complete example (from examples/odes/DTimeMapJetCurveExample.cpp):
#include <iostream>
#include "capd/capdlib.h"
using namespace capd;
int main(){
// maximal degree of truncated Taylor series
const int degree = 4;
// define an instance of class DMap that describes the vector field
DMap pendulum("time:t;par:omega;var:x,dx;fun:dx,sin(omega*t)-sin(x);",degree);
pendulum.setParameter("omega",1.);
int order = 20;
DCnOdeSolver solver(pendulum,order);
DCnTimeMap timeMap(solver);
double initTime = 4;
double finalTime = 10;
// specify initial condition
DJet initJet(2,degree);
// set initial conditions
// main ODE
initJet(0) = 1; initJet(1) = 2;
// init matrix
initJet(0,0) = 1; initJet(0,1) = 1; initJet(1,0) = 3; initJet(1,1) = 4;
// init hessian
initJet(0,0,0) = 1; initJet(0,0,1) = 3; initJet(0,1,1) = 3; initJet(1,0,0) = 0; initJet(1,0,1) = 2; initJet(1,1,1) = 1;
// higher order coefficients remain zero
DJet J(2,degree);
// define functional object
DCnTimeMap::SolutionCurve solution(initTime);
// integrate
timeMap(finalTime,initJet,J,solution);
// print result
std::cout << "solution(7)=" << solution.jet(7).toString() << std::endl;
std::cout << "solution(10)="<< J.toString() << std::endl;
return 0;
}
/* output:
solution(7)=
value :
{0,0} : {-0.451828,-1.67523}
Taylor coefficients of order 1 :
{1,0} : {10.1969,-3.06072}
{0,1} : {13.1894,-3.86086}
Taylor coefficients of order 2 :
{2,0} : {42.6352,17.2114}
{1,1} : {114.667,42.4899}
{0,2} : {73.2755,26.8991}
Taylor coefficients of order 3 :
{3,0} : {145.438,187.641}
{2,1} : {605.182,743.839}
{1,2} : {812.539,969.994}
{0,3} : {353.441,416.38}
Taylor coefficients of order 4 :
{4,0} : {389.64,1142.97}
{3,1} : {2220.14,6172.14}
{2,2} : {4601.83,12302}
{1,3} : {4123.34,10734.3}
{0,4} : {1351.71,3462}
solution(10)=
value :
{0,0} : {1.20484,1.17041}
Taylor coefficients of order 1 :
{1,0} : {-7.80566,-2.6189}
{0,1} : {-10.001,-3.48359}
Taylor coefficients of order 2 :
{2,0} : {-29.1662,-14.8471}
{1,1} : {-79.9856,-38.9769}
{0,2} : {-52.1961,-24.2494}
Taylor coefficients of order 3 :
{3,0} : {-142.07,-97.0567}
{2,1} : {-575.954,-389.605}
{1,2} : {-757.59,-513.922}
{0,3} : {-323.487,-223.165}
Taylor coefficients of order 4 :
{4,0} : {-554.512,-714.868}
{3,1} : {-3064.1,-3836.76}
{2,2} : {-6216.23,-7630.76}
{1,3} : {-5495.57,-6668.73}
{0,4} : {-1790.96,-2161.71}
*/
