This file implements the functions of the semi-implicit integration method and polynomial extrapolation in C. Based on the work of Bader and Deuflhard. Specific functions of the algorithm. More...

#include "bader.h"
#include "mbs_message.h"
#include "mbs_bader.h"
#include "mbs_errors_names.h"
#include <math.h>

Functions
int	bader (double y[], double dydx[], int nv, double xx, double htry, double eps, double yscal[], double hdid, double hnext, double h_max, int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn dd)
	Semi - implicit extrapolation step for integrating stiffo.d.e. More...

int	simpr (double y[], double dydx[], int nvar, double xs, double htot, int nstep, double yout[], int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn *dd)
	Performs one step of semi-implicit midpoint rule. More...

int	mmid (double y[], double dydx[], int nvar, double xs, double htot, int nstep, double yout[], int(derivs)(double, double[], double[], MbsData , MbsDirdyn ), MbsData s, MbsDirdyn *dd)
	Modified midpoint step. More...

void	pzextr (int iest, double xest, double yest[], double yz[], double dy[], int nv, MbsData s, MbsDirdyn dd)
	Use polynomial extrapolation to evaluate nv functions at x = 0 by fitting a polynomial to a sequence of estimates with progressively smaller values x = xest, and corresponding function vectors yest[1..nv].This call is number iest in the sequence of calls.Extrapolated function values are output as yz[1..nv], and their estimated error is output as dy[1..nv]. More...

void	rzextr (int iest, double xest, double yest[], double yz[], double dy[], int nv, MbsData s, MbsDirdyn dd)
	Exact substitute for pzextr, but uses diagonal rational function extrapolation instead of polynomial extrapolation. More...

Variables
static double	dmaxarg1

static double	dmaxarg2

static double	maxarg1

static double	maxarg2

static double	sqrarg

static double	dminarg1

static double	dminarg2

static double	minarg1

static double	minarg2

Detailed Description

This file implements the functions of the semi-implicit integration method and polynomial extrapolation in C. Based on the work of Bader and Deuflhard. Specific functions of the algorithm.

Creation date: April 2018

Author: Sebastien Timmermans

\source Bader, G., and Deuflhard, A semi-implicit mid-point rule for stiff systems of ordinary differential equations, P. 1983, Numerische Mathematik, vol. 41, pp. 373398.

(c) Universite catholique de Louvain

Function Documentation

◆ bader()

int bader	(	double	y[],
		double	dydx[],
		int	nv,
		double *	xx,
		double	htry,
		double	eps,
		double	yscal[],
		double *	hdid,
		double *	hnext,
		double	h_max,
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Semi - implicit extrapolation step for integrating stiffo.d.e.

s, with monitoring of local truncation error to adjust stepsize.Input are the dependent variable vector y[1..nv] and its derivative dydx[1..nv] at the starting value of the independent variable x.Also input are the stepsize to be attempted htry, the required accuracy eps, and the vector yscal[1..nv] against which the error is scaled. On output, y and x are replaced by their new values, hdid is the stepsize that was actually accomplished, and hnext is the estimated next stepsize.d erivs is a user - supplied routine that computes the derivatives of the right - hand side with respect to x, while jacobn(a fixed name) is a user - supplied routine that computes the Jacobi matrix of derivatives of the right - hand side with respect to the components of y.Be sure to set htry on successive steps to the value of hnext returned from the previous step, as is the case if the routine is called by odeint.

Parameters

y	Initial y values, index starting at 0
dydx	Initial y derivative values, index starting at 0
nv	Dimension of y
xx	Value of time at this time step
htry	The total step to be taken
eps	The required accuracy
yscal	The vector [0..n-1] against which the error is scaled, index starting at 0
hdid	The stepsize that was actually accomplished
hnext	The estimated next stepsize
h_max	The maximum stepsize allowed
derivs	The function computing f'
s	The computed MBS structure
dd	The associated MbsDirdyn structure

Returns: error_code Return <0 if something went wrong, 0 otherwise

if mmid is called, compute Jacobian (dfdx and dfdy)

◆ mmid()

int mmid	(	double	y[],
		double	dydx[],
		int	nvar,
		double	xs,
		double	htot,
		int	nstep,
		double	yout[],
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Modified midpoint step.

At xs, input the depend variable vector y [1 ... nvar] and its derivative vector dydx[1...nvar]. Also input is htot, the total step to be made, and nstep, the number of substeps to be used. The output is returned as yout[1...nvar], which need not be a distinct array from y; if it is distinct, however, then y and dydx are returned undamaged.

Parameters

y	Initial y values, index starting at 0
dydx	Initial y derivative values, index starting at 0
nvar	Dimension of y
xs	Value of time at this time step
htot	The total step to be taken
nstep	The number of substeps to be used
yout	solution of the incremented variables
derivs	The function computing f'
s	The computed MBS structure
dd	The associated MbsDirdyn structure

Returns: error_code Return the error, <0 if problem, 0 if everything went well

◆ pzextr()

void pzextr	(	int	iest,
		double	xest,
		double	yest[],
		double	yz[],
		double	dy[],
		int	nv,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Use polynomial extrapolation to evaluate nv functions at x = 0 by fitting a polynomial to a sequence of estimates with progressively smaller values x = xest, and corresponding function vectors yest[1..nv].This call is number iest in the sequence of calls.Extrapolated function values are output as yz[1..nv], and their estimated error is output as dy[1..nv].

Parameters

iest	Number of this call in the sequence of calls
xest	Value of x at which to evaluate the function
yest	Corresponding function vectors
yz	Extrapolated function values output vector
dy	Estimated error output vector
nv	Dimension of y, number of functions
s	The computed MBS structure
dd	The associated MbsDirdyn structure

◆ rzextr()

void rzextr	(	int	iest,
		double	xest,
		double	yest[],
		double	yz[],
		double	dy[],
		int	nv,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Exact substitute for pzextr, but uses diagonal rational function extrapolation instead of polynomial extrapolation.

◆ simpr()

int simpr	(	double	y[],
		double	dydx[],
		int	nvar,
		double	xs,
		double	htot,
		int	nstep,
		double	yout[],
		int()(double, double[], double[], MbsData , MbsDirdyn *)	derivs,
		MbsData *	s,
		MbsDirdyn *	dd
	)

Performs one step of semi-implicit midpoint rule.

Input are the dependent variable y[1 ... nvar], its derivative dydx[1...nvar], the derivative of the right-hand side with respect to x, dfdx[1..nvar] and the jacobian dfdy[1...nvar][1...nvar] at xs. Also input are htot, the total step to be taken, and nstep, the number of substeps to be used. The output is returned as yout[1...nvar], derivs is the user-supplied routine that calculate dydx.

Parameters

y	Initial y values, index starting at 0
dydx	Initial y derivative values, index starting at 0
nvar	Dimension of y
xs	Value of time at this time step
htot	The total step to be taken
nstep	The number of substeps to be used
yout	solution of the incremented variables
derivs	The function computing f'
s	The computed MBS structure
dd	The associated MbsDirdyn structure

Returns: error_code Return the error, <0 if problem, 0 if everything went well

calls an external function

Variable Documentation

◆ dmaxarg1

double dmaxarg1

static

◆ dmaxarg2

double dmaxarg2

static

◆ dminarg1

double dminarg1

static

◆ dminarg2

double dminarg2

static

◆ maxarg1

double maxarg1

static

◆ maxarg2

double maxarg2

static

◆ minarg1

double minarg1

static

◆ minarg2

double minarg2

static

◆ sqrarg

double sqrarg

static

Functions

Variables

Detailed Description

Function Documentation

◆ bader()

◆ mmid()

◆ pzextr()

◆ rzextr()

◆ simpr()

Variable Documentation

◆ dmaxarg1

◆ dmaxarg2

◆ dminarg1

◆ dminarg2

◆ maxarg1

◆ maxarg2

◆ minarg1

◆ minarg2

◆ sqrarg