To solve a system of ordinary differential equations. More...

#include <ODESolver.h>

Public Member Functions
	ODESolver ()
	Default constructor.

	ODESolver (size_t nb_eq)
	Constructor providing the number of equations.

	ODESolver (TimeScheme s, real_t time_step=theTimeStep, real_t final_time=theFinalTime, size_t nb_eq=1)
	Constructor using time discretization data.

	~ODESolver ()
	Destructor.

void	set (TimeScheme s, real_t time_step=theTimeStep, real_t final_time=theFinalTime)
	Define data of the differential equation or system.

void	setNbEq (size_t n)
	Set the number of equations [Default: `1`].

void	setCoef (real_t a0, real_t a1, real_t a2, real_t f)
	Define coefficients in the case of a scalar differential equation.

void	setCoef (string a0, string a1, string a2, string f)
	Define coefficients in the case of a scalar differential equation.

void	setLinear ()
	Claim that ODE is linear.

void	setF (string f)
	Set time derivative, given as an algebraic expression, for a nonlinear ODE.

void	setF (string f, int i)
	Set time derivative, given as an algebraic expression, for a nonlinear ODE.

void	setDF (string df, int i, int j)
	Set time derivative of the function defining the ODE.

void	setdFdt (string df, int i)
	Set time derivative of the function defining the ODE.

void	setRK4RHS (real_t f)
	Set intermediate right-hand side vector for the Runge-Kutta method.

void	setRK4RHS (Vect< real_t > &f)
	Set intermediate right-hand side vector for the Runge-Kutta method.

void	setInitial (Vect< real_t > &u)
	Set initial condition for a first-oder system of differential equations.

void	setInitial (real_t u, int i)
	Set initial condition for a first-oder system of differential equations.

void	setInitial (Vect< real_t > &u, Vect< real_t > &v)
	Set initial condition for a second-order system of differential equations.

void	setInitialRHS (Vect< real_t > &f)
	Set initial RHS for a system of differential equations.

void	setInitial (real_t u, real_t v)
	Set initial condition for a second-order ordinary differential equation.

void	setInitial (real_t u)
	Set initial condition for a first-order ordinary differential equation.

void	setInitialRHS (real_t f)
	Set initial right-hand side for a single differential equation.

void	setMatrices (DMatrix< real_t > &A0, DMatrix< real_t > &A1)
	Define matrices for a system of first-order ODEs.

void	setMatrices (DMatrix< real_t > &A0, DMatrix< real_t > &A1, DMatrix< real_t > &A2)
	Define matrices for a system of second-order ODEs.

void	seODEVectors (Vect< real_t > &a0, Vect< real_t > &a1)
	Define matrices for an implicit nonlinear system of first-order ODEs.

void	seODEVectors (Vect< real_t > &a0, Vect< real_t > &a1, Vect< real_t > &a2)
	Define matrices for an implicit nonlinear system of second-order ODEs.

void	setRHS (Vect< real_t > &b)
	Set right-hand side vector for a system of ODE.

void	setRHS (real_t f)
	Set right-hand side for a linear ODE.

void	setRHS (string f)
	Set right-hand side value for a linear ODE.

void	setNewmarkParameters (real_t beta, real_t gamma)
	Define parameters for the Newmarxk scheme.

void	setConstantMatrix ()
	Say that matrix problem is constant.

void	setNonConstantMatrix ()
	Say that matrix problem is variable.

void	setLinearSolver (Iteration s=DIRECT_SOLVER, Preconditioner p=DIAG_PREC)
	Set linear solver data.

void	setMaxIter (int max_it)
	Set maximal number of iterations.

void	setTolerance (real_t toler)
	Set tolerance value for convergence.

real_t	runOneTimeStep ()
	Run one time step.

void	run (bool opt=false)
	Run the time stepping procedure.

size_t	getNbEq () const
	Return number of equations.

LinearSolver &	getLSolver ()
	Return LinearSolver instance.

real_t	getTimeDerivative (int i=1) const
	Get time derivative of solution.

void	getTimeDerivative (Vect< real_t > &y) const
	Get time derivative of solution (for a system)

real_t	get () const
	Return solution in the case of a scalar equation.

void	getPhase (Vect< real_t > &x, Vect< real_t > &v, size_t i=1)
	Get phase portrait vectors.

Detailed Description

To solve a system of ordinary differential equations.

The class ODESolver enables solving by a numerical scheme a system or ordinary differential equations taking one of the forms:

A linear system of differential equations of the first-order:
A₁(t)u'(t) + A₀(t)u(t) = f(t)
A linear system of differential equations of the second-order:
A₂(t)u''(t) + A₁(t)u'(t) + A₀(t)u(t) = f(t)
A system of ordinary differential equations of the form:
u'(t) = f(t,u(t))

The following time integration schemes can be used:

Forward Euler scheme (value: FORWARD_EULER) for first-order systems
Backward Euler scheme (value: BACKWARD_EULER) for first-order linear systems
Crank-Nicolson (value: CRANK_NICOLSON) for first-order linear systems
Heun (value: HEUN) for first-order systems
2nd Order Adams-Bashforth (value: AB2) for first-order systems
4-th order Runge-Kutta (value: RK4) for first-order systems
2nd order Backward Differentiation Formula (value: BDF2) for linear first-order systems
Newmark (value: NEWMARK) for linear second-order systems with constant matrices

Author: Rachid Touzani

Copyright: GNU Lesser Public License

Constructor & Destructor Documentation

◆ ODESolver()

ODESolver	(	TimeScheme	s,
		real_t	time_step = `theTimeStep`,
		real_t	final_time = `theFinalTime`,
		size_t	nb_eq = `1`
	)

Constructor using time discretization data.

Parameters

[in]	s	Choice of the scheme: To be chosen in the enumerated variable Scheme (see the presentation of the class)
[in]	time_step	Value of the time step. This value will be modified if an adaptive method is used. The default value for this parameter if the value given by the global variable `theTimeStep`
[in]	final_time	Value of the final time (time starts at 0). The default value for this parameter is the value given by the global variable `theFinalTime`
[in]	nb_eq	Number of differential equations (size of the system) [Default: `1`]

Member Function Documentation

◆ getPhase()

void getPhase	(	Vect< real_t > &	x,
		Vect< real_t > &	v,
		size_t	i = `1`
	)

Get phase portrait vectors.

This function gets vectors containing solution and its time derivative to determine phase portraite of the ode. The function is valid for a single ode only.

Parameters

x	Vector containing solution of the ode at each computed time step
v	Vector containing discrete time derivative of the ode at each computed time step
i	Component for which the phase is extracted [Dafault: `1`]

◆ getTimeDerivative() [1/2]

real_t getTimeDerivative ( int i = 1 ) const

Get time derivative of solution.

Return approximate time derivative of solution in the case of a single equation

Parameters

[in] i Index of component whose time derivative is sought

Returns: Time derivative of the i-th component of the solution

Remarks: If we are solving one equation, this parameter is not used.

◆ getTimeDerivative() [2/2]

void getTimeDerivative ( Vect< real_t > & y ) const

Get time derivative of solution (for a system)

Get approximate time derivative of solution in the case of an ODE system

Parameters

[out] y Vector containing time derivative of solution

◆ run()

void run ( bool opt = false )

Run the time stepping procedure.

Parameters

[in] opt Flag to say if problem matrix is constant while time stepping (true) or not (Default value is false)

Note: This argument is not used if the time stepping scheme is explicit

◆ runOneTimeStep()

real_t runOneTimeStep ( )

Run one time step.

Returns: Value of new time step if this one is updated

◆ seODEVectors() [1/2]

void seODEVectors	(	Vect< real_t > &	a0,
		Vect< real_t > &	a1
	)

Define matrices for an implicit nonlinear system of first-order ODEs.

The system has the nonlinear implicit form a1(u)' + a0(u) = 0 Vectors a0, a1 are given as references to class Vect.

Parameters

[in]	a0	Reference to vector in front of the 0-th order term (no time derivative)
[in]	a1	Reference to vector in front of the 1-st order term (first time derivative)

Remarks: This function has to be called at each time step

◆ seODEVectors() [2/2]

void seODEVectors	(	Vect< real_t > &	a0,
		Vect< real_t > &	a1,
		Vect< real_t > &	a2
	)

Define matrices for an implicit nonlinear system of second-order ODEs.

The system has the nonlinear implicit form a2(u)'' + a1(u)' + a0(u) = 0 Vectors a0, a1, a2 are given as references to class Vect.

Parameters

[in]	a0	Reference to vector in front of the 0-th order term (no time derivative)
[in]	a1	Reference to vector in front of the 1-st order term (first time derivative)
[in]	a2	Reference to vector in front of the 2-nd order term (second time derivative)

Remarks: This function has to be called at each time step

◆ set()

void set	(	TimeScheme	s,
		real_t	time_step = `theTimeStep`,
		real_t	final_time = `theFinalTime`
	)

Define data of the differential equation or system.

Parameters

[in]	s	Choice of the scheme: To be chosen in the enumerated variable Scheme (see the presentation of the class)
[in]	time_step	Value of the time step. This value will be modified if an adaptive method is used. The default value for this parameter if the value given by the global variable `theTimeStep`
[in]	final_time	Value of the final time (time starts at 0). The default value for this parameter is the value given by the global variable `theFinalTime`

◆ setCoef() [1/2]

void setCoef	(	real_t	a0,
		real_t	a1,
		real_t	a2,
		real_t	f
	)

Define coefficients in the case of a scalar differential equation.

This function enables giving coefficients of the differential equation as an algebraic expression of time t (see the function fparse)

Parameters

[in]	a0	Coefficient of the 0-th order term
[in]	a1	Coefficient of the 1-st order term
[in]	a2	Coefficient of the 2-nd order term
[in]	f	Value of the right-hand side

Note: Naturally, the equation is of the first order if a2=0

◆ setCoef() [2/2]

void setCoef	(	string	a0,
		string	a1,
		string	a2,
		string	f
	)

Define coefficients in the case of a scalar differential equation.

Parameters

[in]	a0	Coefficient of the 0-th order term
[in]	a1	Coefficient of the 1-st order term
[in]	a2	Coefficient of the 2-nd order term
[in]	f	Value of the right-hand side

Note: Naturally, the equation if of the first order if a2=0

◆ setConstantMatrix()

void setConstantMatrix ( )

Say that matrix problem is constant.

This is useful if the linear system is solved by a factorization method but has no effect otherwise

◆ setDF()

void setDF	(	string	df,
		int	i,
		int	j
	)

Set time derivative of the function defining the ODE.

This function enables prescribing the value of the 1-st derivative for a 1st order ODE or the 2nd one for a 2nd-order ODE. It is to be used for nonlinear ODEs of the form y'(t) = f(t,y(t)) or y''(t) = f(t,y(t),y'(t))
In the case of a system of ODEs, this function can be called once for each equation, given in the order of the unknowns

◆ setdFdt()

void setdFdt	(	string	df,
		int	i
	)

Set time derivative of the function defining the ODE.

This function enables prescribing the value of the 1-st derivative for a 1st order ODE or the 2nd one for a 2nd-order ODE. It is to be used for nonlinear ODEs of the form y'(t) = f(t,y(t)) or y''(t) = f(t,y(t),y'(t))
In the case of a system of ODEs, this function can be called once for each equation, given in the order of the unknowns

◆ setF() [1/2]

void setF ( string f )

Set time derivative, given as an algebraic expression, for a nonlinear ODE.

This function enables prescribing the value of the 1-st derivative for a 1st order ODE or the 2nd one for a 2nd-order ODE. It is to be used for nonlinear ODEs of the form y'(t) = f(t,y(t)) or y''(t) = f(t,y(t),y'(t))
In the case of a system of ODEs, this function can be called once for each equation, given in the order of the unknowns

Parameters

[in] f Expression of the function

◆ setF() [2/2]

void setF	(	string	f,
		int	i
	)

Set time derivative, given as an algebraic expression, for a nonlinear ODE.

This function enables prescribing the value of the 1-st derivative for a 1st order ODE or the 2nd one for a 2nd-order ODE. It is to be used for nonlinear ODEs of the form y'(t) = f(t,y(t)) or y''(t) = f(t,y(t),y'(t))
This function is to be used for the i-th equation of a system of ODEs

Parameters

[in]	f	Expression of the function
[in]	i	Index of equation. Must be not larger than the number of equations

◆ setInitial() [1/5]

void setInitial ( real_t u )

Set initial condition for a first-order ordinary differential equation.

Parameters

[in] u Initial condition (unknown) value

◆ setInitial() [2/5]

void setInitial	(	real_t	u,
		int	i
	)

Set initial condition for a first-oder system of differential equations.

Parameters

[in]	u	Initial condition for an unknown
[in]	i	Index of the unknown

◆ setInitial() [3/5]

void setInitial	(	real_t	u,
		real_t	v
	)

Set initial condition for a second-order ordinary differential equation.

Parameters

[in]	u	Initial condition (unknown) value
[in]	v	Initial condition (time derivative of the unknown) value

◆ setInitial() [4/5]

void setInitial ( Vect< real_t > & u )

Set initial condition for a first-oder system of differential equations.

Parameters

[in] u Vector containing initial condition for the unknown

◆ setInitial() [5/5]

void setInitial	(	Vect< real_t > &	u,
		Vect< real_t > &	v
	)

Set initial condition for a second-order system of differential equations.

Giving the right-hand side at initial time is somtimes required for high order methods like Runge-Kutta

Parameters

[in]	u	Vector containing initial condition for the unknown
[in]	v	Vector containing initial condition for the time derivative of the unknown

◆ setInitialRHS() [1/2]

void setInitialRHS ( real_t f )

Set initial right-hand side for a single differential equation.

Parameters

[in] f Value of right-hand side at initial time. This value is helpful for high order methods

◆ setInitialRHS() [2/2]

void setInitialRHS ( Vect< real_t > & f )

Set initial RHS for a system of differential equations.

Giving the right-hand side at initial time is somtimes required for high order methods like Runge-Kutta

Parameters

[in] f Vector containing right-hand side at initial time. This vector is helpful for high order methods

◆ setLinear()

void setLinear ( )

Claim that ODE is linear.

Claim that the defined ODE (or system of ODEs) is linear

◆ setLinearSolver()

void setLinearSolver	(	Iteration	s = `DIRECT_SOLVER`,
		Preconditioner	p = `DIAG_PREC`
	)

Set linear solver data.

Parameters

[in]	s	Solver identification parameter. To be chosen in the enumeration variable Iteration: DIRECT_SOLVER, CG_SOLVER, CGS_SOLVER, BICG_SOLVER, BICG_STAB_SOLVER, GMRES_SOLVER, QMR_SOLVER [Default: `DIRECT_SOLVER`]
[in]	p	Preconditioner identification parameter. To be chosen in the enumeration variable Preconditioner: IDENT_PREC, DIAG_PREC, ILU_PREC [Default: `DIAG_PREC`]

Note: The argument p has no effect if the solver is DIRECT_SOLVER

◆ setMatrices() [1/2]

void setMatrices	(	DMatrix< real_t > &	A0,
		DMatrix< real_t > &	A1
	)

Define matrices for a system of first-order ODEs.

Matrices are given as references to class DMatrix.

Parameters

[in]	A0	Reference to matrix in front of the 0-th order term (no time derivative)
[in]	A1	Reference to matrix in front of the 1-st order term (first time derivative)

Remarks: This function has to be called at each time step

◆ setMatrices() [2/2]

void setMatrices	(	DMatrix< real_t > &	A0,
		DMatrix< real_t > &	A1,
		DMatrix< real_t > &	A2
	)

Define matrices for a system of second-order ODEs.

Matrices are given as references to class DMatrix.

Parameters

[in]	A0	Reference to matrix in front of the 0-th order term (no time derivative)
[in]	A1	Reference to matrix in front of the 1-st order term (first time derivative)
[in]	A2	Reference to matrix in front of the 2-nd order term (second time derivative)

Remarks: This function has to be called at each time step

◆ setMaxIter()

void setMaxIter ( int max_it )

Set maximal number of iterations.

This function is useful for a non linear ODE (or system of ODEs) if an implicit scheme is used

Parameters

[in] max_it Maximal number of iterations [Default: 100]

◆ setNbEq()

void setNbEq ( size_t n )

Set the number of equations [Default: 1].

This function is to be used if the default constructor was used

◆ setNewmarkParameters()

void setNewmarkParameters	(	real_t	beta,
		real_t	gamma
	)

Define parameters for the Newmarxk scheme.

Parameters

[in]	beta	Parameter beta [Default: `0.25`]
[in]	gamma	Parameter gamma [Default: `0.5`]

◆ setNonConstantMatrix()

void setNonConstantMatrix ( )

Say that matrix problem is variable.

This is useful if the linear system is solved by a factorization method but has no effect otherwise

◆ setRHS() [1/2]

void setRHS ( real_t f )

Set right-hand side for a linear ODE.

Parameters

[in] f Value of the right-hand side for a linear ordinary differential equation

◆ setRHS() [2/2]

void setRHS ( Vect< real_t > & b )

Set right-hand side vector for a system of ODE.

Parameters

[in] b Vect instance containing right-hand side for a linear system of ordinary differential equations

◆ setRK4RHS() [1/2]

void setRK4RHS ( real_t f )

Set intermediate right-hand side vector for the Runge-Kutta method.

Parameters

[in] f Value of right-hand side

◆ setRK4RHS() [2/2]

void setRK4RHS ( Vect< real_t > & f )

Set intermediate right-hand side vector for the Runge-Kutta method.

Parameters

[in] f right-hand side vector

◆ setTolerance()

void setTolerance ( real_t toler )

Set tolerance value for convergence.

This function is useful for a non linear ODE (or system of ODEs) if an implicit scheme is used

Parameters

[in] toler Tolerance value [Default: 1.e-8]

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ ODESolver()

Member Function Documentation

◆ getPhase()

◆ getTimeDerivative() [1/2]

◆ getTimeDerivative() [2/2]

◆ run()

◆ runOneTimeStep()

◆ seODEVectors() [1/2]

◆ seODEVectors() [2/2]

◆ set()

◆ setCoef() [1/2]

◆ setCoef() [2/2]

◆ setConstantMatrix()

◆ setDF()

◆ setdFdt()

◆ setF() [1/2]

◆ setF() [2/2]

◆ setInitial() [1/5]

◆ setInitial() [2/5]

◆ setInitial() [3/5]

◆ setInitial() [4/5]

◆ setInitial() [5/5]

◆ setInitialRHS() [1/2]

◆ setInitialRHS() [2/2]

◆ setLinear()

◆ setLinearSolver()

◆ setMatrices() [1/2]

◆ setMatrices() [2/2]

◆ setMaxIter()

◆ setNbEq()

◆ setNewmarkParameters()

◆ setNonConstantMatrix()

◆ setRHS() [1/2]

◆ setRHS() [2/2]

◆ setRK4RHS() [1/2]

◆ setRK4RHS() [2/2]

◆ setTolerance()