Home | Plots | FAQ | Documents | Simdem | Download |

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
- Quadrature of convolution integrals
- Models with swap over points or discontinuities

- Introduction
- Global storage and sub-models
- Special functions
- Operations with vectors
- Integer and logical functions
- 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

(See Bardsley and Prasad in Computers Chem. 21, 71-82, 1997)

Any serious method for supplying models as ASCII text files should use a readily understood language, it must accept all usual maths operations, including trig 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.

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.

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. k) 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)

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

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.

Contents 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 command quad The command 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

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.

Advice:

Adaptive quadrature cannot be expected to work correctly if the range of integration includes sharp spikes or long stretches of near-zero 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.

- 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 - i) Any sub-models can be used for i and j as long as i and j are less than or equal to 10
- ii) The lower limit for integration must be blim(1) whatever i or j are
- iii) The top limit for integration must be tlim(1) whatever i or j are
- iv) Extra parameters may need to be set to normalise the integral

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.

0 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.

0 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.

- 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.

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

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

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.

0.001 x 0.999 middle invn(x)to avoid nonzero IFAIL returns.

x 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.

10 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.

z 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).

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.

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.

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.

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).

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.

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 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.

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.

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.

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 three 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.

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

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

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.

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.

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

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

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

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

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)

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

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.

% Example: user supplied function of 1 variable ... a straight line ............. p(1) + p(2)*x ............. % 1 equation 1 variable 2 parameters % p(1) p(2) x multiply add f(1) % Here is the same model but with comments to help you understand 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: user supplied function of 1 variable ... single exponential .... f(x) = p(1)*exp[p(2)*x] .... .... % 1 equation 1 variable 2 parameters % p(2) x multiply exponential p(1) multiply f(1) % Here is the same model but with comments to help you understand the rules ========================================================================= % Example: user supplied function of 1 variable ... single exponential .... f(x) = p(1)*exp[p(2)*x] .... .... % 1 equation no. eqns. 1 variable no. variables (or differential equation) 2 parameters no. 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: 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 1 variable 3 parameters % x p(1) subtract p(2) divide normalcdf p(3) multiply f(1) % Here is the same model but with comments to help you understand 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: 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: 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 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) %

% 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)}

% 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: 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) %

Back to the Simfit home page