This header defines specific integrators functions in C. More...

#include "mbs_data.h"
#include "mbs_aux.h"

Enumerations
enum	{ RK4, Dopri5, Rosenbrock, EulerEx, Eulaire, EulerIm, Bader, WMethods, AlphaM, Custom }
	Enumeration of the possible integrator for dirdyn. More...

Functions
void	set_integrator (MbsDirdyn *dd)
	Set the function pointer in MbsDirdyn. More...

void	print_warnings_integrator (MbsData mbs_data, MbsDirdyn mbs_dd, int type_of_integrator)
	Check the options set by user to warn him when he modified an unused option for the integrator. More...

void	print_warnings_constant_step_integrator (MbsDirdyn mbs_dd, char integrator_name)
	Print the warning message, called by print_warnings_integrator for constant step size dt integrators. More...

void	print_warnings_explicit_integrator (MbsDirdyn mbs_dd, char integrator_name)
	Print the warning message, called by print_warnings_integrator for explicit integrators. More...

int	rk4 (double y[], double dydx[], int n, double x, double h, double yout[], int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn *dd)
	Runge Kutta 4 integrator implementation Given values for the variables y[1..n] and their derivatives dydx[1..n] known at x, use the fourth-order Runge-Kutta method to advance the solution over an interval h and return the incremented variables as yout[1..n], which need not be a distinct array from y. More...

int	rosenbrock (int n, int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), double x, double y[], double eps, double hmax, double htry, long nmax, double dydx[], double yscal[], double hnext, MbsData s, MbsDirdyn dd, double hdid)
	Fourth-order Rosenbrock step for integrating stiff problems, with monitoring of local truncation error to adjust stepsize. More...

int	ThetaSC (double y[], double dydx[], int n, double x, double h, double yout[], int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn *dd)
	Unknown integrator. More...

int	mbs_estim_jac_acc (double x, double htry, double y[], double dydx[], int compute_dfdx, double dfdx[], double *dfdy, int n, int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn dd)
	Evaluate the Jacobian of a function accelerations using finite difference. More...

int	mbs_freeze_jac (int freeze_index, int next_freeze_index, double x, double h, double y[], double dydx[], int compute_dfdx, double dfdx[], double dfdy, double dydx_freeze[], double dfdx_freeze[], double dfdy_freeze, int n, int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn mbs_dd)
	Freeze the Jacobian of an integrator structure until mbs_dd->options->n_freeze. More...

int	euler_implicit (double y[], double dydx[], int n, double x, double h, double yout[], int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn mbs_dd)
	Euler Implicit integrator implementation. More...

int	w_methods (double y[], double dydx[], int n, double x, double h, double yout[], int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn *mbs_dd)
	W Methods integrator implementation. More...

Detailed Description

This header defines specific integrators functions in C.

Creation date: May 2018

Author: Sebastien Timmermans, Olivier Lantsoght

(c) Universite catholique de Louvain

Enumeration Type Documentation

◆ anonymous enum

anonymous enum

Enumeration of the possible integrator for dirdyn.

// Set the integrator in main : mbs_dirdyn->options->integrator = Dopri5; mbs_dirdyn->options->resfilename = "dirdyn_Dopri5";

Adding a custom integrator can be done using the "Custom" keyword and modifying the files mbs_custom.c / .h

Enumerator
RK4
Dopri5
Rosenbrock
EulerEx
Eulaire
EulerIm
Bader
WMethods
AlphaM
Custom

Function Documentation

◆ euler_implicit()

int euler_implicit	(	double	y[],
		double	dydx[],
		int	n,
		double *	x,
		double	h,
		double	yout[],
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	mbs_dd
	)

Euler Implicit integrator implementation.

@source (adapted from) Arnold M. et al, Linearly implicit time integration methods in real-time applications: DAEs and stiff ODEs Multibody System Dynamics, 2007, 17:99â117

Parameters

y	state vector of variable
dydx	derivative of y by x
n	number of variables
x	current time value, at which the function needs to be computed
h	time step
yout	solution of the incremented variables
derivs	the function computing f'
s	the MbsData structure of the model on which dirdyn analysis is computed.
mbs_dd	the MbsDirdyn structure related to the integration.

Returns: error_code Return a negative number if something went wrong, 0 otherwise

calls an external function

A = (- h* Jv - h*h * Jp )

A = ( I - * Jv - h * h * Jp)

add 1/h in diagonal terms (identity matrix)

Compute the product : Jp * Xv

B = qdd + h * Jp * Xv

delta p = Xv + h * delta v

with delta v = B

Xv (n+1) with delta v = B

Xp (n+1)

◆ mbs_estim_jac_acc()

int mbs_estim_jac_acc	(	double	x,
		double	htry,
		double	y[],
		double	dydx[],
		int	compute_dfdx,
		double	dfdx[],
		double **	dfdy,
		int	n,
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Evaluate the Jacobian of a function accelerations using finite difference.

Modify and do not restore dydx. WARNING: All vector/matrix indexes start at 0! Update the value of dfdx and dfdy.

Parameters

x	current time value, at which the function needs to be computed
htry	time step
y	state vector of variable
compute_dfdx	flag to ON to compute the dfdx terms, default = 1 ; if compute_dfdx = 0 => dfdx is filled with zeroes
dydx	derivative of y by x
dfdx	derivative of f by x
dfdy	derivative of f by y
n	number of variables
derivs	the function computing f'
s	the MbsData structure of the model on which dirdyn analysis is computed.
dd	the MbsDirdyn structure related to the integration.

Returns: Error status, <0 in case of failure.

◆ mbs_freeze_jac()

int mbs_freeze_jac	(	int	freeze_index,
		int *	next_freeze_index,
		double	x,
		double	h,
		double	y[],
		double	dydx[],
		int	compute_dfdx,
		double	dfdx[],
		double **	dfdy,
		double	dydx_freeze[],
		double	dfdx_freeze[],
		double **	dfdy_freeze,
		int	n,
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	mbs_dd
	)

Freeze the Jacobian of an integrator structure until mbs_dd->options->n_freeze.

WARNING: All vector/matrix indexes start at 0!

Parameters

freeze_index	current count of refreezing steps (if 0, recompute Jac)
next_freeze_index	[out] pointor to next refreezing step
x	current time value, at which the function needs to be computed
h	time step
y	state vector of variable
dydx	derivative of y by x
compute_dfdx	flag to ON to compute the dfdx terms, default = 1 ; if compute_dfdx = 0 => dfdx is filled with zeroes
dfdx	derivative of f by x
dfdy	derivative of f by y
dydx_freeze	pointor to freezed derivative of y by x
dfdx_freeze	pointor to freezed derivative of f by x
dfdy_freeze	pointor to freezed derivative of f by y
n	number of variables
derivs	the function computing f'
s	the MbsData structure of the model on which dirdyn analysis is computed.
mbs_dd	the MbsDirdyn structure related to the integration.

Returns: Error status, <0 in case of failure.

◆ print_warnings_constant_step_integrator()

void print_warnings_constant_step_integrator	(	MbsDirdyn *	mbs_dd,
		char *	integrator_name
	)

Print the warning message, called by print_warnings_integrator for constant step size dt integrators.

Parameters

mbs_dd	MbsDirdyn structure
integrator_name	string to print in the warning message, corresponding to the integrator (RK4, EulerEx, ...)

◆ print_warnings_explicit_integrator()

void print_warnings_explicit_integrator	(	MbsDirdyn *	mbs_dd,
		char *	integrator_name
	)

Print the warning message, called by print_warnings_integrator for explicit integrators.

Parameters

mbs_dd	MbsDirdyn structure
integrator_name	string to print in the warning message, corresponding to the integrator (RK4, EulerEx, ...)

◆ print_warnings_integrator()

void print_warnings_integrator	(	MbsData *	mbs_data,
		MbsDirdyn *	mbs_dd,
		int	type_of_integrator
	)

Check the options set by user to warn him when he modified an unused option for the integrator.

Parameters

mbs_data	MbsData structure
mbs_dd	MbsDirdyn structure
type_of_integrator	ID if the used integrator (RK4, Dopri5...)

◆ rk4()

int rk4	(	double	y[],
		double	dydx[],
		int	n,
		double	x,
		double	h,
		double	yout[],
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Runge Kutta 4 integrator implementation Given values for the variables y[1..n] and their derivatives dydx[1..n] known at x, use the fourth-order Runge-Kutta method to advance the solution over an interval h and return the incremented variables as yout[1..n], which need not be a distinct array from y.

The user supplies the routine derivs(x,y,dydx), which returns derivatives dydx at x.

Parameters

y	state vector of variable
dydx	derivative of y by x
n	number of variables
x	current time value, at which the function needs to be computed
h	time step
yout	solution of the incremented variables
derivs	The function computing f'
s	the MbsData structure of the model on which dirdyn analysis is computed.
dd	the MbsDirdyn structure related to the integration.

Returns: Error status, <0 in case of failure.

◆ rosenbrock()

int rosenbrock	(	int	n,
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		double *	x,
		double	y[],
		double	eps,
		double	hmax,
		double	htry,
		long	nmax,
		double	dydx[],
		double	yscal[],
		double *	hnext,
		MbsData *	s,
		MbsDirdyn *	dd,
		double *	hdid
	)

Fourth-order Rosenbrock step for integrating stiff problems, with monitoring of local truncation error to adjust stepsize.

Input are the dependent variable vector y[1..n] and its derivative dydx[1..n] at the starting value of the independent variable x. WARNING: this feature has not been heavily tested !

@source: H. H. Rosenbrock, "Some general implicit processes for the numerical solution of differential equations", The Computer Journal (1963) 5(4): 329-330 Shampine, L.F. 1982, ACM Transactions on Mathematical Software, vol. 8, pp. 93â113

Parameters

n	Dimension of y
derivs	The function computing f'
x	Value of x
y	Initial y values, index starting at 0
eps	The required accuracy
hmax	The maximum stepsize allowed
htry	The stepsize to be attempted
nmax	The maximal number of allowed steps to reach tolerances
dydx	Initial y derivative values, index starting at 0
yscal	The vector [0..n-1] against which the error is scaled, index starting at 0
hnext	The estimated next stepsize
s	The MbsData structure of the model on which dirdyn analysis is computed.
dd	The MbsDirdyn structure related to the integration.
hdid	The stepsize that was actually accomplished

Returns: error_code Return a negative number if something went wrong, 0 otherwise

calls an external function

◆ set_integrator()

void set_integrator ( MbsDirdyn * dd )

Set the function pointer in MbsDirdyn.

Parameters

dd	MbsDirdyn structure to adapt

◆ ThetaSC()

int ThetaSC	(	double	y[],
		double	dydx[],
		int	n,
		double	x,
		double	h,
		double	yout[],
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Unknown integrator.

Returns: Error status, <0 in case of failure.

◆ w_methods()

int w_methods	(	double	y[],
		double	dydx[],
		int	n,
		double	x,
		double	h,
		double	yout[],
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	mbs_dd
	)

W Methods integrator implementation.

@source (adapted from) Arnold M. et al, Linearly implicit time integration methods in real-time applications: DAEs and stiff ODEs Multibody System Dynamics, 2007, 17:99â117

Parameters

y	state vector of variable
dydx	derivative of y by x
n	number of variables
x	current time value, at which the function needs to be computed
h	time step
yout	solution of the incremented variables
derivs	the function computing f'
s	the MbsData structure of the model on which dirdyn analysis is computed.
mbs_dd	the MbsDirdyn structure related to the integration.

Returns: error_code Return a negative number if something went wrong, 0 otherwise

calls an external function

A = (- h* Jv - h*h * Jp )

A = ( I - * Jv - h * h * Jp)

add 1/h in diagonal terms (identity matrix)

partial

useless if s<2

sum alpha + gamma

Jp * Xv

Complete right-hand term B

B = f + h * Jp * Xv gamma +

delta p = Xv + h * delta v

Enumerations

Functions

Detailed Description

Enumeration Type Documentation

◆ anonymous enum

Function Documentation

◆ euler_implicit()

◆ mbs_estim_jac_acc()

◆ mbs_freeze_jac()

◆ print_warnings_constant_step_integrator()

◆ print_warnings_explicit_integrator()

◆ print_warnings_integrator()

◆ rk4()

◆ rosenbrock()

◆ set_integrator()

◆ ThetaSC()

◆ w_methods()