User defined models
The SIMFIT package has a comprehensive library of models for simulating or fitting, which will cover most situations. However, you may have to construct your own user-defined models in cases involving procedures like the following ones.
- Systems of nonlinear differential equations and Jacobians
- Sets of linked or independent nonlinear equations
- Evaluation of special functions
- Probability integrals and inverses
- Integrating one or more functions of one or several variables
- Locating roots of nonlinear equations in one or several variables
- Constrained nonlinear optimization using partial derivatives
- Evaluating or fitting convolution integrals
- Models with swap over points or discontinuities
- Introduction
- Global storage and sub-models
- Special functions
- Operations with vectors
- Integer and logical functions
- Using standard expressions
- Summary of commands
- Examples
- Example 1: straight line
- Example 2: single exponential
- Example 3: normal integral
- Example 4: capillary diffusion
- Example 5: damped simple harmonic motion
- Example 6: Lotka-Voterra predator-prey differential equations
- Example 7: convolution integral
- Example 8: Rosenbrock's function
- Example 9: helix in parametric form
1. Introduction
Any serious method for supplying models as ASCII text files should use a readily understood language, it must accept all usual maths operations, including trigonometric and hyperbolic functions, and it should allow the use of special functions and calls to other library (DLL) routines. The SIMFIT curve fitting and simulation programs can use such models, as long as the ASCII files are formatted correctly and the syntax has been checked using program usermod.Using standard expressions (New feature at Version 7.1.6)
Standard mathematical expressions can be used at any time in the model defining section of a user-supplied file as long as they are contained in a begin{expression} ... end{expression} structure. For example the code
begin{expression}
f(1) = p(1) + p(2)x + p(3)x^2 + p(4)x^3
end{expression}
would define model number 1, i.e. f(1), as a cubic using four variable parameters
p(1), p(2), p(3), and p(4), and an independent variable x. When a model file is
analysed by SIMFIT it accepts the input to create a virtual stack in reverse Polish,
that is, last in first out, or post fix notation. However, when a begin{expression}
token is encountered it transforms the equations encountered from standard mathematical
equations into reverse Polish until an end{expression} token closes the environment.
So, for most users, there would never be any need to learn reverse Polish and models
can simply be created in usual mathematical notation using program USERMOD.
Note that models involving more advanced operations such as quadrature, root finding, Chebyshev expansions, convolution integrals, vector algebra, or conditional execution, etc. may require combining standard notation with explicit reverse Polish.
Using reverse Polish (i.e., postfix, or calculator) notation
When developing your own model, you must understand that a reverse Polish (last-in-first-out, or post-fix) procedure works by making a stack such that, at the end of the stack operations, the stack has only one element, which is the value of the model. The test files have comments (to the right of the commands, or following the % characters) to help you understand the way models of any degree of complexity can be constructed.The functions that can be used (see usermod1.tf1, etc.)
In addition to all the standard mathematical functions you can use functions such as erfc, the gamma function, normal integrals, etc. The list has now been extended to include calls to numerical routines, loops, conditionals, iterations, etc.
The SIMFIT method is unique in that it allows the user to optimise the execution stack by using duplicate, pop and so on, and by allowing many functions to be defined at the same time. However, since few people program calculators or write in PostScript nowadays, it may seem hard to realise that anything can be programmed without using parentheses if you use post-fix. Many students with no programming experience have been able to write files for systems of differential equations after a few hours practise so it can't be that hard. Start with usermod1.tf1, which is a simple line y = mx + c. Read the file then input it into program usermod to evaluate, plot, find zeros, estimate areas, etc. and it will all suddenly become clear. Just do it.
Rules for supplying a user defined function for simulation or fitting
Please observe the use of the special symbol % in this file. The symbol % starting a line is an escape sequence to indicate a change in the meaning of the input stream, e.g. from text to parameters, from parameters to model instructions, from model to Jacobian, etc. Characters occurring more than eight places after the first non-blank character are interpreted as comments and text here is ignored when the model is parsed. The % symbol MUST be used to indicate:-
i) start of the file
ii) start of the model parameters
iii) start of the model equations
iv) end of model equations (start of Jacobian with diff. eqns.)
The file you supply must have EXACTLY the format now described.
- The file must start with a % symbol indicating where text starts The next lines must be the name/details you choose for the model. This would normally be at least 4 not greater than 24 lines. This text is only to identify the model and is not used by simfit. The end of this section is marked by a % symbol. The next three lines define the type of model.
- The first of these lines must indicate the number of equations in the model, e.g. 1 equation, 2 equations, 3 equations, etc.
- The next must indicate the number of independent variables as in:- 1 variable, 2 variables, 3 variables, etc. or else it could be differential equation to indicate that the model is one or a set of ordinary differential equations with one independent variable.
- The next line must define the number of parameters in the model.
- With differential equations, the last parameters are reserved to set the values for the integration constants y0(i), which can be either estimated or fixed as required. For example, if there are n equations and m parameters are declared in the file, only m-n can be actually used in the model, since y0(i) = p(m-n+i) for i = 1, 2, ..., n.
- Lines are broken up into tokens by spaces.
- Only the first token in each line matters after the model starts.
- Comments begin with % and are added just to explain what's going on.
- Usually the comments beginning with a % can be omitted.
- Critical lines starting with % must be present as explained above. The model operations then follow, one per line until the next line starting with a % character indicates the end of the model.
- Numbers can be in any format, e.g. 2, 1.234, 1.234E-6, 1.234E6
- The symbol f(i) indicates that model equation i is evaluted at this point.
- Differential equations can define the Jacobian after defining
the model. If there are n differential equations of the form
dy(i)/dx = f(i)(x, y(1), y(2), ..., y(n))then the symbol y(i) is used to put y(i) on the stack and there must be a n by n matrix defined in the following way. The element J(a,b) is indicated by putting j(n*(b-1) + a) on the stack. That is the columns are filled up first. For instance with 3 equations you would have a Jacobian J(i,j) = df(i)/dy(j) defined by the sequence:J(1,1) = j(1), J(1,2) = j(4), J(3,1) = j(7) J(2,1) = j(2), J(2,2) = j(5), J(3,2) = j(8) J(3,1) = j(3), J(3,2) = j(6), J(3,3) = j(9)
Permitted operations and the effects produced
x : stack -> stack, x y : stack -> stack, y z : stack -> stack, z add : stack, a, b -> stack, (a + b) subtract : stack, a, b -> stack, (a - b) multiply : stack, a, b -> stack, (a*b) divide : stack, a, b -> stack, (a/b) p(i) : stack -> stack, p(i) ... i can be 1, 2, 3, etc f(i) : stack, a -> stack ...evaluate model since now f(i) = a power : stack, a, b -> stack, (a^b) squareroot : stack, a -> stack, sqrt(a) exponential : stack, a -> stack, exp(a) tentothepower : stack, a -> stack, 10^a ln (or log) : stack, a -> stack, ln(a) log10 : stack, a -> stack, log(a) (to base ten) pi : stack -> stack, 3.1415927 sine : stack, a -> stack, sin(a) ... radians not degrees cosine : stack, a -> stack, cos(a) ... radians not degrees tangent : stack, a -> stack, tan(a) ... radians not degrees arcsine : stack, a -> stack, arcsin(a) ... radians not degrees arccosine : stack, a -> stack, arccos(a) ... radians not degrees arctangent : stack, a -> stack, arctan(a) ... radians not degrees sinh : stack, a -> stack, sinh(a) cosh : stack, a -> stack, cosh(a) tanh : stack, a -> stack, tanh(a) exchange : stack, a, b -> stack, b, a duplicate : stack, a -> stack, a, a pop : stack, a, b -> stack, a absolutevalue : stack, a -> stack, abs(a) negative : stack, a -> stack , -a minimum : stack, a, b -> stack, min(a,b) maximum : stack, a, b -> stack, max(a,b) gammafunction : stack, a -> stack, gamma(a) lgamma : stack, a -> stack, ln(gamma(a)) normalcdf : stack, a -> stack, phi(a) integral from -infinity to a erfc : stack, a -> stack, erfc(a) y(i) : stack -> stack, y(i) Only diff. eqns. j(i) : stack, a -> stack J(i-(i/n), (i/n)+1) Only diff. eqns. *** : stack -> stack, *** ... *** can be any number
Error handling
As the stack is evaluated, action is taken to avoid underflow, overflow and forbidden operations, like 1/x as x tends to zero or taking the log or square root of a negative number etc.2. Global storage and sub-models
This document describes enhancements to the Simfit user-supplied model options at Version 5.4 release 4.037. It is a supplement to w_readme.f5 (which describes the basic techniques) and will be of interest to users who want to define subsidiary models for function evaluation, composition of functions, root finding, adaptive quadrature, constrained optimisation, conditional branching, etc. to be called as independent procedures from within a main model. Note that Simfit creates a separate virtual stack for each model and procedure, and so the only communication between the stacks is the passing of arguments, retrieval of results and sharing of storage and parameter values if requested.
a) The command put b) The command get c) The command get3 d) The command epsabs e) The command epsrel f) The command blim g) The command tlim h) The command putpar i) The command value j) The commands quad and convolute k) The command root l) The command value3 m) The command order n) The command middle o) Syntax for subsidiary models p) Rules for subsidiary models q) Nesting subsidiary models r) IFAIL values for D01AJF, D01EAF and C05AZF s) Examples of subsidiary models t) Synonyms and abbreviations
a) The command put(i)
This command is used to store the results of a sub-calculation. The top stack element is transferred into storage location i and is available globally to all sub-models called by the main model. The stack length is decreased by one, i.e. the element transferred by put(i), and any existing storage element in location i is overwritten. The main model is checked for use of put(k) with no corresponding get(k) unless the main model calls sub-models. Subsidiary models are not checked, of course.b) The command get(j)
This command is used to get the results from a sub-calculation. The storage element previously placed in location j is copied onto the top of the stack. The stack length is increased by one, i.e. the element copied by get(j), and storage location j is unchanged. The main model is checked for get(k) with no corresponding put(k). Subsidiary models are not checked, of course.c) The command get3(i,j,k)
This command is practically equivalent to if ... elseif ... else. If the top stack element is negative get(i) is invoked, if it is zero (to machine precision) get(j) is invoked and, if it is positive get(k) is invoked. The stack length is unchanged, as the top element is replaced by storage location i, j or k. Storage locations i, j and k are unchanged. The main model is checked for corresponding puts. Subsidiary models are not checked, of course. The command order is provided to increase the versatility of the command get3(.,.,.).d) The command epsabs
This command sets the absolute error for subsequent calculations. The default value of espabs (1.0e-3, i.e. absolute error tolerance) is replaced by the top stack element. The stack length is reduced by one, i.e. the value used to re-define epsabs. If epsabs is used inside a subsidiary model, the effects are confined locally to the submodel.e) The command epsrel
This command sets the relative error for subsequent calculations. The default value of esprel (1.0e-6, i.e. relative error tolerance) is replaced by the top stack element. The stack length is reduced by one, i.e. the value used to re-define epsrel. If epsrel is used inside a subsidiary model, the effects are confined locally to the submodel.f) The command blim(i)
This command sets the lower-limit vector for subsequent calculations. The top stack element is used for blim(i), the bottom limit for root finding or lower limit for quadrature. The stack length is reduced by one, i.e. the value used to define blim(i). If blim is used inside a subsidiary model, the effects are confined locally to the submodel.g) The command tlim(i)
This command sets the upper-limit vector for subsequent calculations. The top stack element is used for tlim(i), the top limit for root finding or upper limit for quadrature. The stack length is reduced by one, i.e. the value used to define tlim(i). If tlim is used inside a subsidiary model, the effects are confined locally to the submodel.h) The command putpar
This communicates parameters p(i) from the main model to a sub-model. It must be used to transfer the current parameter values into global storage locations if it is wished to use them in a subsidiary model. Unless the command putpar is used in a main model, the sub-models have no access to parameter values to enable the command p(i) to add parameter i to the sub-model stack. The stack length is unchanged. Note that the storage locations for putpar are initialised to zero so, if you do not use putpar at the start of the main model, calls to p(i) in subsequent sub-models will not lead to a crash, they will simply use p(i) = 0. The command putpar cannot be used in a subsidiary model, of course. Note that putpar should be only used when absolutely necessary as it slows down model evaluation.i) The command value(i)
This command is used to evaluate a subsidiary model. It uses the current values for independent variables to evaluate subsidiary model i. The stack length is increased by one, as the value returned by function evaluation is added to the top of the stack. The command putpar must be used before value(i) if it wished to use the main parameters in subsidiary model number i.Advice:
It is important to make sure that a subsidiary model is correct, by testing it as a main model, if possible, before using it as a subsidiary model. You must be careful that the independent variables passed to the submodel for evaluation are the ones intended. Of course, value(i) can call sub-models which themselves can call root, and/or quad. At whatever level value(i) is called, the independent variables will be the those at that level of nesting, i.e. top-level from a main model and dummy variables in a sub-model. Expert users will immediately realise that there is nothing to prevent them designing a sub-model that is evaluated without using any parameters or arguments, or perhaps just using gets and/or puts, which gives great scope for ingenuity in defining special functions.
j) The command quad(i)
This command is used to estimate an integral by adaptive quadrature. It uses the epsabs, epsrel, blim and tlim values to integrate model i and place the return value on the stack. The values assigned to the blim and tlim arrays are the limits for the integration. If the model i requires j independent variables, then j blim and tlim values must be set before quad(i) is used. The length of the stack increases by one, the return value placed on the top of the stack. The command putpar must be used before quad(i) if it is wished to use the main parameters in the subsidiary model number i.Advice:
Adaptive quadrature cannot be expected to work correctly if the range values, e.g. the extreme upper tail of a decaying exponential. The routines used (D01AJF and D01EAF) cannot really handle infinite ranges, but excellent results can be obtained using commonsense extreme limits, e.g. several relaxation times for a decaying exponential, rather than an arbitrarily large number. With difficult problems it will probably be necessary to increase epsrel and epsabs.
The special case of convolution integrals = f*g
These are done using the command convolute(i,j) which has the following special features.- i) convolute(i,j) calculates the integral from A = blim(1) to B = tlim(1)
for submodels i (i.e. f(.)) and j (i.e. g(.)) as in
f*g = integral from A to B of f(u)g(B - u)du - ii) Any sub-models can be used for i and j as long as i and j are less than or equal to 10
- iii) The lower limit for integration must be blim(1) whatever i or j are
- iv) The top limit for integration must be tlim(1) whatever i or j are
- v) Extra parameters may need to be set to normalise the integral
k) The command root(i)
This command is used to estimate a zero of a sub-model iteratively. It uses the epsabs, epsrel, blim and tlim values to find a root for model i and places the return value on the stack. The values assigned to blim(1) and tlim(1) are the limits for root location. The length of the stack increases by one, the root value placed on the top of the stack. The command putpar must be used before root(i) if it wished to use the main parameters in the subsidiary model i.Advice:
The limits A = blim(1) and B = tlim(1) are used as starting estimates to bracket the root. If f(A)*f(B) > 0 then the range (A,B) is expanded by up to ten orders of magnitude (without changing blim(1) or tlim(1)) until f(A)*f(B) < 0. If this or any other failure occurs, the root is returned as zero. Note that A and B will not change sign, so you can search for, say, just positive roots. If this is too restrictive, make sure blim(1)*tlim(1) < 0. C05AZF is used and, with difficult problems, it will probably be necessary to increase epsrel.l) The command value3(i,j,k)
This is a very powerful command which is capable of many applications of the form: if ... elseif ... else. If the top stack value is negative value(i) is called, if it is zero (to machine precision), value(j) is called, and if it is positive value(k) is called. It relies on the presence of correctly formatted sub-models i, j and k of course, but the values returned by sub-models i, k and k are arbitrary, as almost any code can be employed in models i, j and k. The top value of the stack is replaced by the value returned by the appropriate sub-model. The command order is provided to increase the scope for using the command value3(.,.,.).m) The command order
This is a mechanism for putting -1, 0 or 1 on the stack to control the effects of the commands get3(.,.,.) and value3(.,.,.). What happens is that, before the command order is used, there must be three numbers on the stack in the order ...,a,w,b where, typically, a and b would be limits, such that a < b, and w would either be below these limits, between these limits, or above these limits. If w =< a then the stack ...,a,w,b will become ...,-1, but if a < w =< b then ...,a,w,b will become ...,0, while, if w > b, then ...,a,w,b will become ...,1. For example, the code0 x 4 order value3(1,2,3) f(1)is used in the test file updownup.mod to define a model that has two swap-over points where the model definition changes, i.e. sub-model 1 for x =< 0, sub-model 2 for 0 < x =< 4, and sub-model 3 for x > 4. Note that the command order takes three arguments off the stack and replaces them with one argument, so the stack length is decreased by two.
n) The command middle
There must be three numbers on the stack ...,a,w,b as for the command order, with a < b. If w < a the stack will become ...,a, if w > b then the stack will become ...,b, otherwise, for a < w =< b the stack will become ...,w. This command is used to prevent underflow, overflow or calculations that would lead to singularities. For example, the code0 x 1 middlewould leave a number w on the stack, where 0 =< w =< 1, and w = x only if 0 =< x =< 1. Note that the command order takes three arguments off the stack and replaces them with one argument, so the stack length is decreased by two.
o) Syntax for subsidiary models
The format for defining sub-models is very strict and must be used exactly as now described. Suppose you want to use n independent equations. Then n user files are developed and, when they have been tested, they are pasted in order directly after the end of the main model, each surrounded by a begin{model(i)} and end{model(i)} command. So, if you want to use a particular model as a sub-model, you first of all develop it using program usermod then, when it is satisfactory, just add it to the main model. However, note that sub-models are subject to several restrictions.p) Rules for subsidiary models
- Sub-model files must be placed in numerically increasing order at the end of the main model file. Model parsing is abandoned if a sub-model occurs out of sequence.
- There must be no spaces or non-model lines between the main model and the subsidiary models, or between any subsequent sub-models.
- Sub-models cannot be differential equations.
- Sub-models of the type being described must define just one equation.
- Sub-models are not tested for consistent put and get commands, since puts might be defined in the main model, etc.
- Sub-models cannot use putpar, since putpar only has meaning in a main model.
- Sub-models can use the commands value(i), root(j) and quad(k), but it is up to users to make sure that all calls are consistent.
- When the command value(i) is used, the arguments passed to the sub-model for evaluation are the independent variables at the level at which the command is used. For instance if the main model uses value(i) then value(i) will be evaluated at the x, y, z, etc. of the main model, but with model(i) being used for evaluation. Note that y(k) must be used for functions with more than three independent variables, i.e. when x, y and z no longer suffice. It is clear that if a model uses value(i), then the number of independent variables in that model must be equal to or greater than the number of independent variables in sub-model(i).
- When the commands root(i) and quad(j) are used, the independent variables in the sub-model numbers i and j are always dummy variables.
- When developing models and subsidiary models independently you may get error messages about x not being used, or a get without a corresponding put. Often these can be suppressed by using a pop until the model is developed. For instance x followed by pop will silence any messages about x not being used, etc.
q) Nesting subsidiary models
The subsidiary models can be considered to be independent except when there is a clash that would lead to recursion. For instance, value(1) can call model(1) which uses root(2) to find a root of model(2), which calls quad(3) to integrate model(3). However, at no stage can there be simultaneous use of the same model as value(k), and/or quad(k), and/or root(k). The same subsidiary model cannot be used by more than any one instance of value/quad/root at the same time. Just commonsense really, virtual stack k for model k can only be used for one function evaluation at a time. Obviously there can be at most one instance each of value, root and quad active simultaneously.r) IFAIL values for D01AJF, D01AEF and C05AZF
Since these iterative techniques may be used inside optimisation or numerical integration procedures, the soft fail option IFAIL = 1 is used. If the Simfit version of these routines is used, a silent exit will occur and failure will not be communicated to users. So it is up to users to be very careful that all calls to quadrature and root finding are consistent and certain to succeed. Default function values of zero are returned on failure.s) Examples
The test file updown.mod illustrates how to cause a model to swap from one definition to another at critical values of the independent variable, updownup.mod illustrates how to construct a model with two swap-over points, while the series of test files usermodx.tf? give examples of various techniques using sub-models for evaluation, root finding and adaptive quadrature.t) Synonyms and abbreviations
The following sets of commands are equivalent:-sub, minus, subtract mul, multiply div, divide sqrt, squarero, squareroot exp, exponent, exponential ten, tentothe, tentothepower ln, log sin, sine cos, cosine tan, tangent asin, arcsin, arcsine acos, arccos, arccosine atan, arctan, arctangent dup, duplicate exch, exchange, swap del, delete, pop abs, absolute neg, negative min, minimum max, maximum phi, normal, normalcd, normalcdf abserr, epsabs relerr, epsrel
3. Special functions
This document describes how to call the special functions commonly used in numerical analysis, statistics, mathematical simulation and data fitting, by one-line commands in Simfit user-supplied models. The options provided at Version 5.4 release 3.039 are described, and a knowledge of the Simfit user-defined model syntax, as described in the files w_readme.f5 and w_readme.f6, is presumed.
Command NAG Description
======= === ===========
arctanh(x) S11AAF Inverse hyperbolic tangent
arcsinh(x) S11AAF Inverse hyperbolic sine
arccosh(x) S11AAF Inverse hyperbolic cosine
ai(x) S17AGF Airy function Ai(x)
dai(x) S17AJF Derivative of Ai(x)
bi(x) S17AHF Airy function Bi(x)
dbi(x) S17AKF Derivative of Bi(x)
besj0(x) S17AEF Bessel function J0
besj1(x) S17AFF Bessel function J1
besy0(x) S17ACF Bessel function Y0
besy1(x) S17ADF Bessel function Y1
besi0(x) S18ADF Bessel function I0
besi1(x) S18AFF Bessel function I1
besk0(x) S18ACF Bessel function K0
besk1(x) S18ADF Bessel function K1
phi(x) S15ABF Normal cdf
phic(x) S15ACF Normal cdf complement
erf(x) S15AEF Error function
erfc(x) S15ADF Error function complement
dawson(x) S15AFF Dawson integral
ci(x) S13ACF Cosine integral Ci(x)
si(x) S13ADF Sine integral Si(x)
e1(x) S13AAF Exponential integral E1(x)
ei(x) ...... Exponential integral Ei(x)
rc(x,y) S21BAF Elliptic integral RC
rf(x,y,z) S21BBF Elliptic integral RF
rd(x,y,z) S21BCF Elliptic integral RD
rj(x,y,z,r) S21BDF Elliptic integral RJ
sn(x,m) S21CAF Jacobi elliptic function SN
cn(x,m) S21CAF Jacobi elliptic function CN
dn(x,m) S21CAF Jacobi elliptic function DN
ln(1+x) S01BAF ln(1 + x) for x near zero
mchoosen(m,n) ...... Binomial coefficient
gamma(x) S13AAF Gamma function
lngamma(x) S14ABF log Gamma function
psi(x) S14ADF Digamma function, (d/dx)log(Gamma(x))
dpsi(x) S14ADF Trigamma function, (d^2/dx^2)log(Gamma(x))
igamma(x,a) S14BAF Incomplete Gamma function
igammac(x,a) S14BAF Complement of Incomplete Gamma function
fresnelc(x) S20ADF Fresnel C function
fresnels(x) S20ACF Fresnel S function
bei(x) S19ABF Kelvin bei function
ber(x) S19AAF Kelvin ber function
kei(x) S19ADF Kelvin kei function
ker(x) S19ACF Kelvin ker function
cdft(x,m) G01EBF cdf for t distribution
cdfc(x,m) G01ECF cdf for chi-square distribution
cdff(x,m,n) G01EDF cdf for F distribution (m = num, n = denom)
cdfb(x,a,b) G01EEF cdf for beta distribution
cdfg(x,a,b) G01EFF cdf for gamma distribution
invn(x) G01FAF inverse normal
invt(x,m) G01FBF inverse t
invc(x,m) G01FCF inverse chi-square
invb(x,a,b) G01FEF inverse beta
invg(x,a,b) G01FFF inverse gamma
spence(x) ...... Spence integral: 0 to x of -(1/y)log|(1-y)|
clausen(x) ...... Clausen integral: 0 to x of -log(2*sin(t/2))
struveh(x,m) ...... Struve H function order m (m = 0, 1)
struvel(x,m) ...... Struve L function order m (m = 0, 1)
kummerm(x,a,b)...... Confluent hypergeometric function M(a,b,x)
kummeru(x,a,b)...... U(a,b,x), b = 1 + n, the logarithmic solution
lpol(x,m,n) ...... Legendre polynomial of the 1st kind, P_n^m(x),
-1 =< x =< 1, 0 =< m =< n
abram(x,m) ...... Abramovitz function order m (m = 0, 1, 2), x > 0,
integral: 0 to infinity of t^m exp( - t^2 - x/t)
debye(x,m) ...... Debye function of order m (m = 1, 2, 3, 4)
(m/x^m)[integral: 0 to x of t^m/(exp(t) - 1)]
fermi(x,a) ...... Fermi-Dirac integral (1/Gamma(1 + a))[integral:
0 to infinity t^a/(1 + exp(t - x))]
heaviside(x,a)...... Heaviside unit impulse function h(x - a)
delta(i,j) ...... Kronecker delta function
impulse(x,a,b)...... Unit impulse function (small b for Dirac delta)
spike(x,a,b) ...... Triangular spike unit impulse function
gauss(x,a,b) ...... Gauss unit impulse function (normal pdf)
sqwave(x,a) ...... Square wave, amplitude 1, period 2a
rtwave(x,a) ...... Rectified triangle, amplitude 1, period 2a
mdwave(x,a) ...... Morse dot wave, amplitude 1, period 2a
stwave(x,a) ...... Sawtooth wave, amplitude 1, period a
rswave(x,a) ...... Rectified sine wave, amplitude 1, period pi/a
shwave(x,a) ...... Rectified half sine, amplitude 1, period 2*pi/a
uiwave(x,a,b) ...... Unit impulse, area 1, period a, width b
Summary
Note that, prior to Version 5.4 release 3.010, only the first eight characters were read of each line, but that has now been increased. Also, to allow users to document their models, from now on all lines starting with a ! within the main model will be ignored and treated as comment lines.
Any of the above commands included as a line in a Simfit model or sub-model simply takes the top stack element as argument and replaces it by the function value. The NAG routines indicated can be consulted for details, as equivalent routines, agreeing very closely with the NAG specifications, are used. The soft fail (IFAIL = 1) options have been used so the simulation will not terminate on error condition, a default value will be returned. Obviously it is up to users to make sure that sensible arguments are supplied, for instance positive degrees of freedom, F or chi-square arguments, etc. To help avoid problems that could arise with arguments outside limits, e.g. probabilities less than zero or greater than one, the command middle (synonym mid) is provided.
Using the command middle to avoid singularities
This command is very useful for avoiding undeflow and overflow, etc. and takes three arguments a, w, b, where a < b. If w =< a then a is placed on the stack, if w > b then b is placed on the stack, otherwise w is placed on the stack. For instance, to obtain an inverse for the normal distribution which requires 0 =< x =< 1, we could use the code0.001 x 0.999 middle invn(x)to avoid nonzero IFAIL returns.
Functions with one argument
The top stack element will be popped and used as an argument, so the routines can be used in several ways. For instance the following codex phi(x) f(1)would simulate model 1 as a normal cdf, while the code
get(4) phi(x) f(3)would return model three as the normal cdf for whatever was stored in storage location 4.
Functions with two arguments
The top stack element is popped and used as x, while the second is popped and used as a, m, or y, as the case may be. For instance the code10 x cdft(x,m)would place the t distribution cdf with 10 degrees of freedom on the stack, while the code
5 0.975 invt(x,m)would place the critical value for a two-tailed t test with 5 degrees of freedom at a confidence level of 95% on the stack.
Functions with three or more arguments
The procedure is a simple extension of that described for functions of two arguments. First the stack is prepared as ...u,v,w,z,y,x but, after the function call, it would be ...u,v,w,f(x,y,z). For example, the codez y x rf(x,y,z) f(11)would return model 11 as the elliptic function RF, since f(11) would have been defined as a function of at least three variables. However, the code
get(3) get(2) get(1) rd(x,y,z) 1 add f(7)would define model(7) as one plus the elliptic function RD evaluated at whatever was stored in locations 3 (i.e. z), 2 (i.e. y) and 1 (i.e. x).
Test files
Three new test files have been supplied to illustrate these commands:usermods.tf1: special functions with one argument usermods.tf2: special functions with two arguments usermods.tf2: special functions with three argumentsThese should be used in program usermod by repeatedly, editing, reading in the edited files, simulating, etc. to explore the options. Users can choose which of the options provided is to be used, by simply uncommenting the desired option and leaving all the others commented. Note that these are all set up for f(1) as a function of one variable and that, by commenting and removing comments so that only one command is active at any one time, the models can be plotted as continuous functions. Alternatively singly calculated values can be compared to tabulated values, which should be indistiguishable if your editing is correct.
4. Operations with vectors
This document describes enhancements to the user-defined model syntax, at version 5.4 release 4.011, designed to increase the scope for model simulation, fitting and plotting where vector arithmetic is involved.
There are techniques to allow initialisation of vectors, i.e. arrays of numerical constants during first pass, commands to accelerate model evaluation, methods to ease the development of models that require vector norms or dot products, and options to simplify repetitive operations such as iterative calculations (i.e. loops).
In particular, 1-line commands to evaluate poynomials or approximate functions by Chebyshev series are described. A knowledge of the Simfit user-defined model syntax described in w_readme.f5, w_readme.f6, and w_readme.f7, is assumed.
The command store(j)
This command is similar to the put(j) command, but there is an important difference; the command put(j) is executed every time the model is evaluated, but the command store(j) is only executed when the model file is parsed for the first time. So store(j) is really equivalent to a data initialisation statement at compile time. For instance, the code
3 store(14)
would initialise store(14) to the value 3. If no further put(14) is used, then storage location 14 would preserve the value 3 for all subsequent calculations in the main model or any sub-model, so that storage location 14 could be regarded as a global constant. Of course any put(14) in the main model or any sub-model would overwrite storage location 14. The main use for the store command is to define special values that are to be used repeatedly for model evaluation, e.g. coefficients for a Chebyshev expansion. For this reason there is another very important difference between put(j) and store(j); store(j) must be preceeded by a literal constant, e.g. 3.2e-6, and cannot be assigned as the end result of a calculation, because storage initialisation is done before calculations.
To summarise: store(j) must be preceded by a numerical value, when it pops this top stack element after copying it into storage location j. So the stack length is decreased by one, to initialise storage location j, but only on the first pass through the model, i.e. when parsing.
The command storef(file)
Since it is tedious to define more than a few storage locations using the command store(j), the command storef(*.*) provides a mechanism for initialising an arbitrary number of storage locations at first pass using contiguous locations. The file specified by the storef command is read and, if it is a Simfit vector file, all the successive components are copied into corresponding storage locations. An example of this is the test model file cheby.mod (and the related data file cheby.dat) which should be run using program usermod to see how a global vector is set up for a Chebshev approximation. Other uses could be when a model involves calculations with a set of fixed observations, such as a time series.
To summarise: the command storef(mydatya) will read the components of any n-dimensional Simfit vector file, mydata, into n successive storage locations starting at position 1, but only when the model file is parsed at first pass. Subsequent use of put(j) or store(j) for j in the range (1,n) will overwrite the previous effect of storef(mydata).
The command poly(x,m,n)
This evaluates m terms of a polynomial by Horner's method of nested multiplication, with terms starting at store(n) and proceeding as far as store(m + n - 1). The polynomial will be of degree m - 1 and it will be evaluated in ascending order. For example, the code
1 store(10) 0 store(11) -1 store(12) 10 3 x poly(x,m,n)
will place the value of f(x) = 1 - x^2 on the stack, where x is the local argument. Of course, the contents of the storage locations can also be set by put(j) commands which would then overwrite the previous store(j) command. For instance, the following code
5 put(12) 10 3 2 poly(x,m,n)
used after the previous code, would now place the value 21 on the stack, since f(t) = 1 + 5t^2 = 21, and the argument is now t = 2.
To summarise: poly(x,m,n) evaluates a polynomial of degree m - 1, using successive storage locations n, n + 1, n + 2, ..., n + m - 1, i.e. the constant term is storage location n, and the coefficient of degree m - 1 is storage location m + n - 1. The argument is whatever value is on the top of the stack when poly(x,m,n) is invoked. This command takes three arguments x, m, n off the stack and replaces them by one value, so the stack is decreased by two elements. If there is an error in m or n, e.g. m or n negative, there is no error message, and the value f(x) = 0 is returned.
The command cheby(x,m,n)
The Chebyshev coefficients are first stored in locations n to n + m - 1, then the command cheby(x,m,n) will evaluate a Chebyshev expansion using the Broucke method with m terms. Note that the first term must be twice the constant term since, as usual, only half the constant is used in the expansion. This code, for instance, will return the Chebyshev approximation to exp(x).2.532132 store(20) 1.130318 store(21) 0.271495 store(22) 0.044337 store(23) 0.005474 store(24) 20 5 x cheby(x,m,n)Note that, if the numerical values are placed on the stack sequentially, then they obviously must be peeled off in reverse order, as follows.
2.532132 1.130318 0.271495 0.044337 0.005474 store(24) store(23) store(22) store(21) store(20) 20 5 x cheby(x,m,n)To summarise: cheby(x,m,n) evaluates a Chebyshev approximation with m terms, using successive storage locations n, n + 1, n + 2, ..., n + m - 1, i.e. twice the constant term is in storage location n, and the coefficient of T(m - 1) is in storage location m + n - 1. The argument is whatever value is on the top of the stack when cheby(x,m,n) is invoked. This command takes three arguments x, m, n off the stack and replaces them by one value, so the stack is decreased by two elements. If there is an error in x, m or n, e.g. x not in (-1,1), or m or n negative, there is no error message, and the value f(x) = 0 is returned. Use test model cheby.mod with program usermod to appreciate this command.
The commands l1norm(m,n), l2norm(m,n) and linorm(m,n)
The Lp norms are calculated for a vector with m terms, starting at storage location n, i.e. l1norm calculates the sum of the absolute values, l2norm calculates the Euclidean norm, while linorm calculates the infinity norm (that is, the largest absolute value in the vector).
It should be emphasised that l2norm(m,n) puts the Euclidean norm on the stack, that is the length of the vector (the square root of the sum of squares of the elements) and not the square of the distance. For example, the code
2 put(5) -4 put(6) 3 put(7) 4 put(8) 1 put(9) l1norm(3,5)
would place 9 on the stack, while the command l2norm(5,5) would put 6.78233 on the stack, and the command linorm(5,5) would return 4.
To summarise: these commands take two arguments off the stack and calculate either the sum of the absolute values, the square root of the sum of squares, or the largest absolute value in m successive storage locations starting at location n. The stack length is decreased by one since m and n are popped and replaced by the norm. There are no error messages and, if an error is encountered, a zero value is returned.
The commands sum(m,n) and ssq(m,n)
As there are occasions when it is useful to be able to add up the signed values or the squares of values in storage, these commands are provided. For instance, the code
1 2 3 4 put(103) put(102) put(101) put(100) 100 4 sum(m,n) f(1) 101 3 ssq(m,n) f(2)
would assign 10 to function 1 and 29 to function 2.
To summarise: these commands take two arguments off the stack and then replace them with either the sum of m storage locations starting at position n, or the sum of squares of m storage locations starting at position n, so decreasing the stack length by 1.
The command dotprod(l,m,n)
This calculates the scalar product of two vectors of length l which are stored in successive locations starting at positions m and n.
To summarise: The stack length is decreased by 2, as three arguments are consumed, and the top stack element is then set equal to the dot product, unless an error is encountered when zero is returned.
Commands to use mathematical constants
The following commands are provided to faciltate model building.
pi = 3.141592653589793e+00 ... pi
piby2 = 1.570796326794897e+00 ... pi divided by two
piby3 = 1.047197551196598e+00 ... pi divided by three
piby4 = 7.853981633974483e-01 ... pi divided by four
twopi = 6.283185307179586e+00 ... pi multiplied by two
root2pi = 2.506628274631000e+00 ... square root of two pi
deg2rad = 1.745329251994330e-02 ... degrees to radians
rad2deg = 5.729577951308232e+01 ... radians to degrees
root2 = 1.414213562373095e+00 ... square root of two
root3 = 1.732050807568877e+00 ... square root of three
eulerg = 5.772156649015329e-01 ... Euler's gamma
lneulerg =-5.495393129816448e-01 ... log (Euler's gamma)
To summarise: these constants are merely added passively to the stack
and do not affect any existing stack elements. To use the constants,
the necessary further instructions would be required. So, for instance,
to transform degrees into radial measure, the code
94.25 deg2rad multiplywould replace the 94.25 degrees by the equivalent radians.
5. Integers and logical functions
Enhancements to user defined model commands at release 4.016, involving integer and logical functions, and dynamically selected user-defined model functions with arbitrary arguments.Integer functions
Sometimes integers are needed in models, for instance, as exponents, as summation indices, as logical flags, as limits in do loops, or as pointers in case constructs, etc.
So there are special integer functions that take the top argument off the stack whatever number it is (say x) then replace it by an appropriate integer as follows.
int(x): replace x by the integer part of x nint(x): replace x by the nearest integer to x sign(x): replace x by -1 if x < 0, by 0 if x = 0, or by 1 if x > 0
When using integers with simfit models it must be observed that only double precision floating point numbers are stored, and all calculations are done with such numbers, so that 0 actually means 0.0 to machine precision. So, for instance, when using these integer functions with real arguments to create logicals or indices for summation, etc. the numbers on the stack that are to be used as logicals or integers are actually transformed dynamically into integers when required at run time, using the equivalent of nint(x) to generate the appropriate integers. Because of this, you should note that code such as
... 11.3 19.7 1.2 int(x) 2.9 nint(x) divide
would result in 1.0/3.0 being added to the stack (i.e. 0.3333...) and not 1/3 (i.e 0) as it would for true integer division, leading to the stack
..., 11.3, 19.7, 0.3333333
Logical functions
Logical variables are stored in the global storage vector as either 1.0 (so that nint(x) = 1 = true) or as 0.0 (so that nint(x) = 0 = false). The logical functions either generate logical variables by testing the magnitude of the arbitrary stack value (say x) with respect to zero (to machine precision) or they accept only logical arguments (say m or or n) and return an error message if the stack values are not pre-set to 0.0 or 1.0. Note that logical variables (i.e. Booleans) can be stored using put(i) and retrieved using get(i), so that logical tests of any order of complexity can be constructed.lt0(x) replace x by 1 if x < 0, otherwise by 0 le0(x) replace x by 1 if x =< 0, otherwise by 0 eq0(x) replace x by 1 if x = 0, otherwise by 0 ge0(x) replace x by 1 if x >= 0, otherwise by 0 gt0(x) replace x by 1 if x > 0, otherwise by 0 not(m) replace m by NOT(m), error if m not 0 or 1 and(m,n) replace m and n by AND(m,n), error if m or n not 0 or 1 or(m,n) replace m and n by OR(m,n), error if m or n not 0 or 1 xor(m,n) replace m and n by XOR(m,n), error if m or n not 0 or 1
Conditional execution
Using these integer and logical functions in an appropriate sequence interspersed by put(.) and get(.) commands, any storage location (say j) can be set up to test whether any logical condition is true or false. So, the commands if(.) and ifnot(.) are provided to select model features depending on logical variables. The idea is to calculate the logical variables using the integer and logical functions, then load them into storage using put(.) commands. The if(.) and ifnot(.) commands inspect the designated storage locations and return 1 if the storage location has the value 1.0 (to machine precision), or 0 otherwise, even if the location is not 0.0 (to machine precision). The logical values returned model code is executed whereas, if a 0 is returned, the next line is missed out.if(j) execute the next line only if storage(j) = 1.0 ifnot(j) execute the next line only if storage(j) = 0.0Note that very extensive logical tests and blocks of code for conditional executions, do loops, while and case constructs can be generated by using these logical functions sequentially but, because not all the lines of code will be active, the parsing routines will indicate the number of if(.) and ifnot(.) commands and the resulting potentially unused lines of code. This information is not important for correctly formatted models, but it can be used to check or debug code if required.
Consult the test file if.mod to see how to use logical functions.
Arbitrary functions with arbitrary arguments
The submodels described so far for evaluation, integration, root finding, supplied by the simfit calls to the model evaluation subroutines. However, sometimes it is useful to be able to evaluate a submodel with arbitrary arguments added to the stack, or arguments that are functions of the main arguments. Again, it is useful to be able to evaluate an arbitrary function chosen dynamically from a set of submodels indexed by an integer parameter calculated at run time, rather than read in at compile time when the model is first parsed. So, to extend the user-defined model syntax, the command user1(x,m) is provided. The way this works involves three steps:1) an integer (m) is put on the stack to denote the required model, 2) calculations are performed to put the argument (x) on the stack, then 3) the user defined model is called and the result placed on the stack.For instance the code
... 14.7 3 11.3 user1(x,m)would result in
..., 14.7, 12.5if the value of submodel number 3 is 12.5 at an argument of 11.3.
Similar syntax is used for functions of two and three variables, i.e.
user1(x,m) user2(x,y,m) user3(x,y,z,m)Clearly the integer m can be literal, calculated or retrieved from storage, but it must correspond to a submodel that has been defined in the sequence of submodels, and the calculated arbitrary arguments x, y, z must be sensible values. For instance the commands
2 x user1(x,m)and
value(2)are equivalent. However the first form is more versatile, as the model number (m, 2 in this case) and argument (x, the dummy value in this case) can be altered dynamically as the result of stack calculations, while the second form will always invoke the case with m = 2 and x = the subroutine argument to the model.
The model file user1.mod illustrates how to use the user1(x,m) command.
6. Using standard expressions
Standard expressions using mathematical notation can be used at any point inside the model definition section of the user defined model file. For example the code
begin{expression}
f(1) = [p(1) + p(2)x]/[1.0 + p(3)x + p(4)x^2]
end{expression}
would define model number 1 as a 2:2 rational function such as are
used in steady state enzyme kinetics.
Another example could be the lines
%
begin{expression}
f(1) = p(1)y(1) - p(2)y(1)y(2)
f(2) = -p(3)y(2) + p(4)y(1)y(2)
end{expression}
%
begin{expression}
j(1) = p(1) - p(2)y(2)
j(2) = p(4)y(2)
j(3) = -p(2)y(1)
j(4) = -p(3) + p(4)y(1)
end{expression}
%
defining the Lotka-Volterra predator-prey differential equations with Jacobian.
7. Summary of commands
The syntax is defined by Bardsley and Prasad (Computers Chem. 21, 71-82, 1997) and in the Simfit reference manual and the Simfit readme files. Program usermod can test codes and check models for errors, and there are example test model files that can be run to illustrate the various techniques.Commands described in w_readme.f5 (The basic commands)
Command Stack effect ------- ------------ x stack -> stack, x y stack -> stack, y z stack -> stack, z add stack, a, b -> stack, (a + b) subtract stack, a, b -> stack, (a - b) multiply stack, a, b -> stack, (a*b) divide stack, a, b -> stack, (a/b) p(i) stack -> stack, p(i) ... i can be 1, 2, 3, etc f(i) stack, a -> stack ...evaluate model since now f(i) = a power stack, a, b -> stack, (a^b) squareroot stack, a -> stack, sqrt(a) exponential stack, a -> stack, exp(a) tentothepower stack, a -> stack, 10^a ln (or log) stack, a -> stack, ln(a) log10 stack, a -> stack, log(a) (to base ten) pi stack -> stack, 3.1415927 sine stack, a -> stack, sin(a) ... radians not degrees cosine stack, a -> stack, cos(a) ... radians not degrees tangent stack, a -> stack, tan(a) ... radians not degrees arcsine stack, a -> stack, arcsin(a) ... radians not degrees arccosine stack, a -> stack, arccos(a) ... radians not degrees arctangent stack, a -> stack, arctan(a) ... radians not degrees sinh stack, a -> stack, sinh(a) cosh stack, a -> stack, cosh(a) tanh stack, a -> stack, tanh(a) exchange stack, a, b -> stack, b, a duplicate stack, a -> stack, a, a pop stack, a, b -> stack, a absolutevalue stack, a -> stack, abs(a) negative stack, a -> stack , -a minimum stack, a, b -> stack, min(a,b) maximum stack, a, b -> stack, max(a,b) gammafunction stack, a -> stack, gamma(a) lgamma stack, a -> stack, ln(gamma(a)) normalcdf stack, a -> stack, phi(a) integral from -infinity to a erfc stack, a -> stack, erfc(a) y(i) stack -> stack, y(i) Only diff. eqns. j(i) stack, a -> stack J(i-(i/n), (i/n)+1) Only diff. eqns. *** stack -> stack, *** ... *** can be any number
Commands described in w_readme.f6 (global storage and sub-models)
Command Effect
------- ------
begin{model(i)} Start of definition of sub-model i
end{model(i)} End of definition of sub-model i
put(i) Cut top stack element into storage location i
get(i) Copy storage location i to top of stack
get3(i,j,k) Get one of (i,j,k) following use of command order
epsabs Re-define absolute error tolerance
epsrel Re-define relative error tolerance
blim(i) Re-define bottom limit i
tlim(i) Re-define top limit i
putpar Copy parameters into global storage for sub-models
value(i) Evaluate model i
quad(i) Integrate model i
convolute(i,j) Convolute models i and j
root(i) Estimate a root of model i
value3(i,j,k) Evaluate one of (i,j,k) following use of command order
order -1,0,1 depending on values (use before get3 and value3)
middle Reflect a value back between limits if necessary
Commands described in w_readme.f7 (Special functions)
Command NAG Description
------- ---- -----------
arctanh(x) S11AAF Inverse hyperbolic tangent
arcsinh(x) S11AAF Inverse hyperbolic sine
arccosh(x) S11AAF Inverse hyperbolic cosine
ai(x) S17AGF Airy function Ai(x)
dai(x) S17AJF Derivative of Ai(x)
bi(x) S17AHF Airy function Bi(x)
dbi(x) S17AKF Derivative of Bi(x)
besj0(x) S17AEF Bessel function J0
besj1(x) S17AFF Bessel function J1
besy0(x) S17ACF Bessel function Y0
besy1(x) S17ADF Bessel function Y1
besi0(x) S18ADF Bessel function I0
besi1(x) S18AFF Bessel function I1
besk0(x) S18ACF Bessel function K0
besk1(x) S18ADF Bessel function K1
phi(x) S15ABF Normal cdf
phic(x) S15ACF Normal cdf complement
erf(x) S15AEF Error function
erfc(x) S15ADF Error function complement
dawson(x) S15AFF Dawson integral
ci(x) S13ACF Cosine integral Ci(x)
si(x) S13ADF Sine integral Si(x)
e1(x) S13AAF Exponential integral E1(x)
ei(x) ...... Exponential integral Ei(x)
rc(x,y) S21BAF Elliptic integral RC
rf(x,y,z) S21BBF Elliptic integral RF
rd(x,y,z) S21BCF Elliptic integral RD
rj(x,y,z,r) S21BDF Elliptic integral RJ
sn(x,m) S21CAF Jacobi elliptic function SN
cn(x,m) S21CAF Jacobi elliptic function CN
dn(x,m) S21CAF Jacobi elliptic function DN
ln(1+x) S01BAF ln(1 + x) for x near zero
mchoosen(m,n) ...... Binomial coefficient
gamma(x) S13AAF Gamma function
lngamma(x) S14ABF log Gamma function
psi(x) S14ADF Digamma function, (d/dx)log(Gamma(x))
dpsi(x) S14ADF Trigamma function, (d^2/dx^2)log(Gamma(x))
igamma(x,a) S14BAF Incomplete Gamma function
igammac(x,a) S14BAF Complement of Incomplete Gamma function
fresnelc(x) S20ADF Fresnel C function
fresnels(x) S20ACF Fresnel S function
bei(x) S19ABF Kelvin bei function
ber(x) S19AAF Kelvin ber function
kei(x) S19ADF Kelvin kei function
ker(x) S19ACF Kelvin ker function
cdft(x,m) G01EBF cdf for t distribution
cdfc(x,m) G01ECF cdf for chi-square distribution
cdff(x,m,n) G01EDF cdf for F distribution (m = num, n = denom)
cdfb(x,a,b) G01EEF cdf for beta distribution
cdfg(x,a,b) G01EFF cdf for gamma distribution
invn(x) G01FAF inverse normal
invt(x,m) G01FBF inverse t
invc(x,m) G01FCF inverse chi-square
invb(x,a,b) G01FEF inverse beta
invg(x,a,b) G01FFF inverse gamma
spence(x) ...... Spence integral: 0 to x of -(1/y)log|(1-y)|
clausen(x) ...... Clausen integral: 0 to x of -log(2*sin(t/2))
struveh(x,m) ...... Struve H function order m (m = 0, 1)
struvel(x,m) ...... Struve L function order m (m = 0, 1)
kummerm(x,a,b)...... Confluent hypergeometric function M(a,b,x)
kummeru(x,a,b)...... U(a,b,x), b = 1 + n, the logarithmic solution
lpol(x,m,n) ...... Legendre polynomial of the 1st kind, P_n^m(x),
-1 =< x =< 1, 0 =< m =< n
abram(x,m) ...... Abramovitz function order m (m = 0, 1, 2), x > 0,
integral: 0 to infinity of t^m exp( - t^2 - x/t)
debye(x,m) ...... Debye function of order m (m = 1, 2, 3, 4)
(m/x^m)[integral: 0 to x of t^m/(exp(t) - 1)]
fermi(x,a) ...... Fermi-Dirac integral (1/Gamma(1 + a))[integral:
0 to infinity t^a/(1 + exp(t - x))]
heaviside(x,a)...... Heaviside unit impulse function h(x - a)
delta(m,n) ...... Kronecker delta impulse function
impulse(x,a,b)...... Unit impulse function (small b for Dirac delta)
spike(x,a,b) ...... Unit triangular impulse function
gauss(x,a,b) ...... Unit Gauss pdf impulse function
sqwave(x,a) ...... Square wave, amplitude 1, period 2a
rtwave(x,a) ...... Rectified triangle, amplitude 1, period 2a
mdwave(x,a) ...... Morse dot wave, amplitude 1, period 2a
stwave(x,a) ...... Sawtooth wave, amplitude 1, period a
rswave(x,a) ...... Rectified sine wave, amplitude 1, period pi/a
shwave(x,a) ...... Sine half-wave, amplitude 1, period 2*pi/2
uiwave(x,a,b) ...... Unit impulse, area 1, period a, width b
Commands described in w_readme.f8 (Operations with vectors)
Command Effect ------- ------ store(j) Initialise global store(j) storef(file) Initialise global store(1) to global store(n) from a file poly(x,m,n) Evaluate a polynomial cheby(x,m,n) Evaluate a Chebyshev expansion l1norm(m,n) Evaluate a L_1 norm l2norm(m,n) Evaluate a L_2 norm linorm(m,m) Evaluate a L_infinity norm sum(m,n) Evaluate the sum of vector components ssq(m,n) Evaluate the sum of squares of vector components dotprod(l,m,n) Evaluate the scalar product of two vectors
Commands to use mathematical constants
Command Value Comment
------- ----- -------
pi 3.141592653589793e+00 pi
piby2 1.570796326794897e+00 pi divided by two
piby3 1.047197551196598e+00 pi divided by three
piby4 7.853981633974483e-01 pi divided by four
twopi 6.283185307179586e+00 pi multiplied by two
root2pi 2.506628274631000e+00 square root of two pi
deg2rad 1.745329251994330e-02 degrees to radians
rad2deg 5.729577951308232e+01 radians to degrees
root2 1.414213562373095e+00 square root of two
root3 1.732050807568877e+00 square root of three
eulerg 5.772156649015329e-01 Euler's gamma
lneulerg -5.495393129816448e-01 log (Euler's gamma)
Commands described in w_readme.f9 (integers and arbitrary functions)
Note: x, y, z, m and n are the stack values but j is a literal constant, while logicals 1 = true and 0 = false are stored as floating point numbers, to be converted into integers using nint(.) internally when required.Command Effect ------- ------ int(x) x replaced by the integer part of x nint(x) x replaced by the nearest integer to x sign(x) x replaced by -1 if x < 0, by 0 if x = 0 or by 1 if x > 0 lt0(x) x replaced by 1 if x < 0, otherwise by 0 le0(x) x replaced by 1 if x =< 0, otherwise by 0 eq0(x) x replaced by 1 if x = 0, otherwise by 0 ge0(x) x replaced by 1 if x >= 0, otherwise by 0 gt0(x) x replaced by 1 if x > 0, otherwise by 0 not(m) m replaced by not(m), error if m not 0 or 1 and(m,n) m and n replaced by and(m,n), error if m and n not 0 or 1 or(m,n) m and n replaced by or(m,n), error if m and n not 0 or 1 xor(m,n) m and n replaced by xor(m,n), error if m and n not 0 or 1 if(j) execute the next line iff storage(j) = 1.0 ifnot(j) execute the next line iff storage(j) = 0.0 user1(x,m) model m evaluated at x user2(x,y,m) model m evaluated at x, y user3(x,y,z,m) model m evaluated at x, y, z
8. Examples
To appreciate the SIMFIT conventions, examine the following test files:-
deqmod?.tf? ...systems of differential eqns. (require deqpar?.tf?)
usermod1.tf?...functions of 1 independent variable
usermod2.tf?...functions of 2 independent variables
usermod3.tf?...functions of 3 independent variables
usermodd.tf?...one differential eqn.
The files usermod1.tf? are functions of 1 variable, usermod2.tf? are
functions of two variables, usermod3.tf? are functions of 3 variables,
while usermodd.tf? are differential equations. Models deqmod?.tf? are
models for several functions at the same time, e.g. sets of differential
equations as used by deqsol or sqpfit.
For differential equations, the Jacobian must be supplied if you are using program qnfit, but program deqsol will calculate the Jacobian by numerical approximation if you end the model with % then follow with a line starting with a %, instead of the definitions for j(i).
Program usermod can be used to develop, then test a model before trying to use it. The idea is first to examine the test files supplied, then read them into usermod for practise before finally developing your own models. Any file you supply to deqsol, makdat, qnfit, etc. must be formatted exactly like the test files provided.
Model 1 to 6 are first displayed using expressions and then using the less familar reverse Polish notation.
Example 1: straight line
Model as an expression
%
Example: p(1) + p(2)*x
%
1 equation
1 variable
2 parameters
%
begin{expression}
f(1) = p(1) + p(2)x
end{expression}
%
Same model in reverse Polish with comments to explain the rules
% start of text defining model indicated by % Example: user supplied function of 1 variable ... a straight line ............. p(1) + p(2)*x ............. % end of text, start of parameters indicated by % 1 equation number of equations to define the model 1 variable number of variables (or differential equation) 2 parameters number of parameters in the model % end of parameters, start of model indicated by % p(1) put p(1) on the stack: stack = p(1) p(2) put p(2) on the stack: stack = p(1), p(2) x put an x on the stack: stack = p(1), p(2), x multiply multiply top elements: stack = p(1), p(2)*x add add the top elements: stack = p(1) + p(2)*x f(1) evaluate the model f(1) = p(1) + p(2)*x % end of the model definitions indicated by %
Example 2: single exponential
Model as an expression
%
Example: f(x) = p(1)*exp[p(2)*x]
%
1 equation
1 variable
2 parameters
%
begin{expression}
f(1) = p(1)exp[p(2)x]
end{expression}
%
Same model in reverse Polish with comments to explain the rules
%
Example: user supplied function of 1 variable ... single exponential
....
f(x) = p(1)*exp[p(2)*x] ....
....
%
1 equation number of equations
1 variable number of variables (or differential equation)
2 parameters number of parameters in this model
%
p(2) put p(2) on the stack: stack=p(2)
x put an x on the stack: stack=p(2), x
multiply multiply top elements: stack=p(2)*x
exponential exponential operator : stack=exp[p(2)*x]
p(1) put p(1) on the stack: stack=exp[p(2)*x], p(1)
multiply multiply top elements: stack=p(1)*exp[p(2)*x]
f(1) f(1) = p(1)*exp[p(2)*x]
%
Example 3: normal integral
Model as an expression
%
Example: normal integral
%
1 equation
1 variable
3 parameters
%
begin{expression}
f(1) = p(3)normalcdf[(x - p(1))/p(2)]
end{expression}
%
Same model in reverse Polish with comments to explain the rules
%
Example: user supplied function of 1 variable ... normal integral
i.e., integral from -infinity to x of
p(3)/[p(2)*sqrt(2*pi)]*exp{(-1/2)*[(u - p(1))/p(2)]^2}
%
1 equation number of equations
1 variable number of variables (or differential equation)
3 parameters number of parameters in this model
%
x stack=x
p(1) stack=x, p(1)
subtract stack=x - p(1)
p(2) stack=x - p(1), p(2)
divide stack=[x - p(1)]/p(2)
normalcdf stack=integral
p(3) stack=integral, p(3)
multiply stack=p(3)*integral
f(1) evaluate model
%
Example 4: capillary diffusion
Model as an expression
%
Example: f(x) = p(1)*erfc[x/(2*sqrt(p(2))]
%
1 equation
1 variable
2 parameters
%
begin{expression}
f(1) = p(1)erfc[x/(2sqrt(p(2))]
end{expression}
%
Same model in reverse Polish
% Example: user supplied function of 1 variable ... capillary diffusion f(x) = p(1)*erfc[x/(2*sqrt(p(2))] % 1 equation 1 variable 2 parameters % x p(2) squareroot 2 multiply divide erfc p(1) multiply f(1) %
Example 5: damped simple harmonic motion
Model as an expression
%
Example: f(x) = p(4)*exp[-p(3)*x]*cos[p(1)*x - p(2)]
%
1 equation
1 variable
4 parameters
%
begin{expression}
f(1) = p(4)exp[-p(3)x]cos[p(1)x - p(2)]
end{expression}
%
Same model in reverse Polish.
%
Example: user supplied function of 1 variable ... damped SHM
Damped simple harmonic motion in the form
f(x) = p(4)*exp[-p(3)*x]*cos[p(1)*x - p(2)]
where p(i) >= 0
%
1 equation
1 variable
4 parameters
%
p(1)
x
multiply
p(2)
subtract
cosine
p(3)
x
multiply
negative
exponential
multiply
p(4)
multiply
f(1)
%
Example 6: Lotka-Voterra predator-prey differential equations
Model using expressions
%
Example: Lotka-Volterra predator-prey equations
differential equations: f(1) = dy(1)/dx
= p(1)*y(1) - p(2)*y(1)*y(2)
f(2) = dy(2)/dx
= -p(3)*y(2) + p(4)*y(1)*y(2)
jacobian: j(1) = df(1)/dy(1)
= p(1) - p(2)*y(2)
j(2) = df(2)/dy(1)
= p(4)*y(2)
j(3) = df(1)/dy(2)
= -p(2)*y(1)
j(4) = df(2)/dy(2)
= -p(3) + p(4)*y(1)
initial condition: y0(1) = p(5), y0(2) = p(6)
%
2 equations
differential equation
6 parameters
%
begin{expression}
f(1) = p(1)y(1) - p(2)y(1)y(2)
f(2) = -p(3)y(2) + p(4)y(1)y(2)
end{expression}
%
begin{expression}
j(1) = p(1) - p(2)y(2)
j(2) = p(4)y(2)
j(3) = -p(2)y(1)
j(4) = -p(3) + p(4)y(1)
end{expression}
%
The coding for the model is now finished but, optionally, parameter starting
values, curve fitting limits and the range of integration can be appended.
If these are not supplied program DEQSOL will request an initialisation file.
begin{limits}
0 1.0 3
0 1.0 3
0 1.0 3
0 1.0 3
0 1.0 3
0 0.5 3
end{limits}
begin{range}
121
0
10
end{range}
Same model in reverse Polish.
%
Example of a user supplied pair of differential equations
file: deqmod2.tf2 (typical parameter file deqpar2.tf2)
model: Lotka-Volterra predator-prey equations
differential equations: f(1) = dy(1)/dx
= p(1)*y(1) - p(2)*y(1)*y(2)
f(2) = dy(2)/dx
= -p(3)*y(2) + p(4)*y(1)*y(2)
jacobian: j(1) = df(1)/dy(1)
= p(1) - p(2)*y(2)
j(2) = df(2)/dy(1)
= p(4)*y(2)
j(3) = df(1)/dy(2)
= -p(2)*y(1)
j(4) = df(2)/dy(2)
= -p(3) + p(4)*y(1)
initial condition: y0(1) = p(5), y0(2) = p(6)
Note: the last parameters must be y0(i) in differential equations
%
2 equations
differential equation
6 parameters
%
y(1)
y(2)
multiply
duplicate
p(2)
multiply
negative
p(1)
y(1)
multiply
add
f(1)
p(4)
multiply
p(3)
y(2)
multiply
subtract
f(2)
%
p(1)
p(2)
y(2)
multiply
subtract
j(1)
p(4)
y(2)
multiply
j(2)
p(2)
y(1)
multiply
negative
j(3)
p(4)
y(1)
multiply
p(3)
subtract
j(4)
%
Example 7: convolution integral
Model using expressions
%
convolution integral: from 0 to x of f(u)*g(x - u) du, where
f1(t) = f(t) = exp(-p(1)*t)
f2(t) = g(t) = [p(2)^2]*t*exp(-p(2)*t)
f3(t) = f*g = f1*f2
--------------------------------------------------------------
This demonstrates how to define 2 equations as sub-models, using
expressions A, B, C to communicate parameters to the sub-models,
and the command convolute(1,2) to integrate sub-models 1 and 2
(by adaptive quadrature) from blim(1) = 0 to t = tlim(1) = x.
Precision of D01AJF quadrature is controlled by epsabs and epsrel
and blim(1) and tlim(1) must be used for the convolution limits
which, in this case are 0 to x, where x > 0 by assumption.
Note that usually extra parameters must be supplied if it wished
to normalise so that the integral of f or g or f*g is specified
(e.g.,equals 1) over the total range of possible integration.
This must often be done, e.g., if g(.) is a density function.
The gamma distribution normalising factor p(2)**2 is stored as
C in this example to avoid unnecessary re-calculation.
%
3 equations
1 variable
2 parameters
%
begin{expression}
A = p(1)
B = p(2)
C = p(2)*p(2)
end{expression}
1
x
user1(x,m)
f(1)
2
x
user1(x,m)
f(2)
0.0001
epsabs
0.001
epsrel
0
blim(1)
x
tlim(1)
convolute(1,2)
f(3)
%
begin{model(1)}
%
Example: exponential decay, exp(-p(1)*x)
%
1 equation
1 variable
0 parameter
%
begin{expression}
f(1) = exp(-A*x)
end{expression}
%
end{model(1)}
begin{model(2)}
%
Example: gamma density of order 2
%
1 equation
1 variable
0 parameters
%
begin{expression}
f(1) = C*x*exp(-B*x)
end{expression}
%
end{model(2)}
Advice: consult w_readme.f6 for details of how to define sub models
Same model using reverse Polish
%
An example using two sub-models for a convolution integral f*g
=================================================================
This demonstrates how to define 2 equations as sub-models, using
the command putpar to communicate parameters to the sub-models,
and the command convolute(1,2) to integrate sub-models 1 and 2
(by adaptive quadrature) from blim(1) = 0 to t = tlim(1) = x.
Precision of D01AJF quadrature is controlled by epsabs and epsrel
and blim(1) and tlim(1) must be used for the convolution limits
which, in this case are 0 to x, where x > 0 by assumption.
.................................................................
integral: from 0 to x of f(u)*g(x - u) du, where
f(t) = exp(-p(1)*t)
g(t) = [p(2)^2]*t*exp(-p(2)*t)
.................................................................
Note that usually extra parameters must be supplied if it wished
to normalise so that the integral of f or g or f*g is specified
(e.g.,equals 1) over the total range of possible integration.
This must often be done, e.g., if g(.) is a density function.
The gamma distribution normalising factor p(2)**2 is stored in
this example to avoid unnecessary re-calculation.
%
1 equation
1 variable
2 parameters
%
putpar
p(2)
p(2)
multiply
put(1)
0.0001
epsabs
0.001
epsrel
0
blim(1)
x
tlim(1)
convolute(1,2)
f(1)
%
begin{model(1)}
%
Example: exponential decay, exp(-p(1)*x)
%
1 equation
1 variable
1 parameter
%
p(1)
x
multiply
negative
exponential
f(1)
%
end{model(1)}
begin{model(2)}
%
Example: gamma density of order 2
%
1 equation
1 variable
2 parameters
%
p(2)
x
multiply
negative
exponential
x
multiply
get(1)
multiply
f(1)
%
end{model(2)}
Example 8: Rosenbrock's function
Model using expressions
%
Example: Rosenbrock's 2-dimensional function for optimisation
f(1) = 100(y - x^2)^2 + (1 - x)^2
f(2) = g(1)
= d(f(1))/dx
= -400x(y - x^2) - 2(1 - x)
f(3) = g(2)
= d(f(1))/dy
= 200(y - x^2)
%
3 equations
2 variables
0 parameters
%
begin{expression}
A = y - x^2
B = 1 - x
f(1) = 100*A^2 + B^2
f(2) = -400A*x - 2B
f(3) = 200A
end{expression}
%
Same model in reverse Polish
%
Rosenbrock's 2-dimensional function for optimisation
f(1) = 100(y - x^2)^2 + (1 - x)^2
f(2) = g(1)
= d(f(1))/dx
= -400x(y - x^2) - 2(1 - x)
f(3) = g(2)
= d(f(1))/dy
= 200(y - x^2)
%
3 equations
2 variables
0 parameters
%
y
x
x
multiply
subtract
put(1)
1.0
x
subtract
put(2)
get(1)
get(1)
multiply
100.0
multiply
get(2)
get(2)
multiply
add
f(1)
get(1)
x
multiply
-400.0
multiply
get(2)
-2.0
multiply
add
f(2)
get(1)
200.0
multiply
f(3)
%
Example 9: Helix in parameteric form
Model using expressions
%
Example: X = A*cos(t), Y = B*sin(t), Z = C*t
where: t = x, A = p(1), B = p(2), C = p(3)
and X(t) = f(1), Y(t) = f(2), Z(t) = f(3)
%
3 equations
1 variable
3 parameters
%
begin{expression}
f(1) = p(1)cos(x)
f(2) = p(2)sin(x)
f(3) = p(3)x
end{expression}
%
Same model in reverse Polish
%
Example: the helix
X = A*cos(t), Y = B*sin(t), Z = C*t
where: t = x, A = p(1), B = p(2), C = p(3)
and X(t) = f(1), Y(t) = f(2), Z(t) = f(3)
%
3 equations
1 variable
3 parameters
%
x
cos
p(1)
multiply
f(1)
x
sin
p(2)
multiply
f(2)
x
p(3)
multiply
f(3)
%
Simfit home page