SOLVE: linear and nonlinear equations and data fit (H)

NUL

required

• chi() is an

  inline expression

  or a

  callback function

  that SOLVE tries to minimize. This can be used to determine
  • the root(s) of an expression, e.g.:
    • NUL=x^2 -1
  • the minimum of a function, e.g.:
    • NUL=x^2 + 1
  • the parameters of a function to least-square fit some data like
    • NUL = ( theory(nr) - data(nr) ) / error(nr)
  The function can include

  tolerance weights

  for each data point (equation) to account for relative weights or errors, like error(nr) in the example.

Unknown

• scalar argument

  required

  for cases with a single root or a single least-square fit parameter
• x is the initial guess (default 0) on input, and is the root on output
• the Callback function in NUL=... can be entered directly in this case (e.g. NUL=x^2-1)
• AXIS(MiN=..., MaX=...) and LINE(..., Points=...) come handy to visualize the equation
• Example: find a root and a function minimum of a

  quintic equation

  (1 root, 1 equation)
  • chi2 = SOLVE( Unknown=x, NUL=x^5 - x + 1 )
  results usually depend on initial x, e.g.:
  • initial x > -0.6: x=0.6687397766, chi2=0.2162322132 (

    function minimum

    )
  • initial x ≤ -0.6: x=-1.16729231, chi2=9.34876755E-9 (

    root of function

    )

• a vector argument is

  required

  for cases with m>1 roots. The solution vector must be dimensioned with m elements:
  • REAL :: roots(m)
  • SOLVE( Unknown=roots, NUL=null(), DataIdx=nr )
• example of a least-square fit to the Gaussian y=h*EXP(-((x-x0)/w)^2). Measured data are in the vectors x and y
  • parms = ALIAS(h, x0, w) ! for better legibility form a vector from the 3 parameters
  • parms = (8, 5, 2) ! guess some initial parameters for h, x0, w
  • chi2 = SOLVE(Unknown=parms, DATA=m, DataIdx=nr, NUL=1-h*EXP( -((x(nr)-x0)/w)^2)/y(nr) )
• note the normalization of the NUL expression to y(nr) to force equal weights for all data points in the example.
• common practice is also the normalization to the measuring errors.

DataIdx

required

• DataIdx is needed to evaluate the correct NUL=... equation in the callback function
• Example:
  •    XY = ALIAS(x, y) ! same as REAL XY(2) with named elements: x == XY(1) and y == XY(2)
  •    XY = (2, 7) ! define initial values for x and y
  •    chi2 = SOLVE( U=XY, NUL=null(), DataIdx=nr ) ! callback function is null()
  •    ! return values are x=7.497601303, y=2.768650594, chi2=3.7395492E-7
  • END ! of "main"
  • FUNCTION null() ! called by the solver of the SOLVE function. All variables are global
  •    IF(nr == 1) null = x - SIN(x)*COSH(y) ! equation 1
  •    IF(nr == 2) null = y - COS(x)*SINH(y) ! equation 2
  • END ! of function null()

DATA

• required

  for a least-square fit with less roots (parameters p) than equations (data d)
• for each parameter 1...p the function evaluation NUL=... is called d times
• p is the dimension of vec in Unknown=vec

UnknownIdx

• this root index (or parameter index) can help to debug and organize the callback function
• allows the script to perform index-invariant but timely part evaluations by testing np, e.g.:
  • IF(np == 0) THEN
  • ! sub-expressions valid for >1 parameters could go here
  • ENDIF
• np = 0 : return the base-chi needed to approximate the derivatives ∂chi/∂root.

Const

• optional

  to keep selected parameters constant while iterationg unstable models
• c=2^(j-1) for unknown j to be kept constant, e.g.:
  • SOLVE(..., Const=10, ...) ! constant parameters 2 and 4

Limit

• to limit the next root iteration to root/Limit .... root*Limit
• Limit must be > 1
• default is no Limit (denoted by Limit=0)

MaXIters

• maximum number of iterations
• default is MaXIters = 1000

NumDiff

• to approximate the differential ∂chi/∂p ≈ (chi(p+delta) - chi(p)) / delta. This can be useful if the callback function involves approximation techniques like numeric solutions to differential equations
• default NumDiff = 1E-5

Speed1

• starting value of the iteration acceleration factor on success
• default Speed1 = 1000

SpeedStop

• to stop iterations after multiple failures
• default SpeedStop = 1e-10

ERror

• SOLVE(....., ERROR=999)

Iters

• returns in n the actual number of iterations performed

SpeedHigh

• returns the maximum successful acceleration factor

SpeedLow

• returns the minimum successful acceleration factor

STDdev

• returns a RMS error estimate of Unknown(s)
• dp should have the same dimension as vec in Unknown=vec
• meaningful values only for normally distributed errors