# GNU Scientific Library examples rewritten

• Difficulty: 020/100
• Author: jeremy theler
• Keywords: VAR, VECTOR, VECTORS, SIZE, DATA, NUMBER, PRINT, SEP, PRINT_FUNCTION, FORMAT, TEXT, SEP, SEPARATOR, PARAMETRIC, MIN, MAX, TYPE, OUTER_STEPS, PHASE_SPACE, INTERPOLATION, MINIMIZE, ALGORITHM, GUESS, GRADIENT, GRADTOL, TOL, IF, ENDIF, FIT, TO, VIA, COLUMNS, FILE_PATH, VERBOSE, integral, random, .=, interpolation, derivative, root, sqrt, func_min, cos, random, sobol, vecmax, vecmin, cspline, cspline_periodic, conjugate_fr, nmsimplex2, done_outer, step_inner, step_outer, in_static, integral,

As wasora heavily relies on the “GNU Scientific Library”: “http://www.gnu.org/software/gsl/” to perform most of its numerical-related tasks, it should be possible to run most of GSL’s examples with a wasora input. Indeed, this section shows how some of the GSL examples found in the documentation can be implemented by just writing a wasora input file instead of coding an compiling and ad-hoc C program. However, as wasora is not (and does not want to be) a programming language, some of these examples may look awkward but it should be kept in mind that they are written as to resemble as much as possible the original GSL examples written in C.

The motley examples are taken from the “GNU Scientific Library v1.16 documentation”: “http://www.gnu.org/software/gsl/manual/html_node/”, and are sorted in increasing complexity (and not by its appearance order in the manual). Each case refers to the original section of the manual where it is discussed. First, the C code from the official GSL distribution is shown. Then, the corresponding wasora input that should accomplish the same computation is listed and finally a terminal where the problem is solved in both ways. If the case involve the generation of figures, these are included afterwards.

The paragraphs below are a combination of texts that belong to the original GSL documentation with explanations written by the author of wasora. Comments and extensions are welcome.

## integration.was

Consider the following integral, which has an algebraic-logarithmic singularity at the original

!bt [ _0^1 x^{-1/2} (x) , dx = 4 ] !et

The original example uses the QAGS integrator, which is not implemented in wasora. However, the regular integral functional (which corresponds to the QAG integrator) does an acceptable job. Currently, wasora does not provide a way to report the estimated errors.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/integration.c

The corresponding wasora input:

# GSL example 17.14, numerical integration
VAR x
PRINT %.18f integral(x^(-1/2)*log(x),x,0,1)
$gcc -o integration integration.c -lgslcblas -lgsl -lm$ ./integration
result          = -4.000000000000085265
exact result    = -4.000000000000000000
estimated error =  0.000000000000135447
actual error    = -0.000000000000085265
intervals =  8
$wasora integration.was -3.999889192698361295$ 

## diff.was

This example estimates the derivative of the function $$f(x) = x^{3/2}$$ at $$x=2$$ and at $$x=0$$. The function $$f(x)$$ is undefined for $$x<0$$ so the derivative at $$x=0$$ is computed using gsl_deriv_forward in the C example. To reproduce this behavior in wasora, a positive value $$p>0$$ is passed as the optional fifth argument (see the derivative functional documentation in the wasora reference). Again, wasora does not provide a way to estimate the error.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/diff.c

The corresponding wasora input:

# GSL example 28.2, differentiation

f(x) := x^(3/2)

PRINT TEXT "f(x) = x^(3/2)"
PRINT TEXT "x = 2.0"
PRINT TEXT "f'(x) = " %.10f derivative(f(x),x,2,1e-8) SEP " "
PRINT TEXT "exact = " %.10f 1.5*sqrt(2) SEP " "
PRINT
PRINT TEXT "x = 0.0"
PRINT TEXT "f'(x) = " %.10f derivative(f(x),x,0,1e-8,1) SEP " "
PRINT TEXT "exact = " %.10f 0 SEP " "
$gcc -o diff diff.c -lgslcblas -lgsl -lm$ ./diff
f(x) = x^(3/2)
x = 2.0
f'(x) = 2.1213203120 +/- 0.0000005006
exact = 2.1213203436

x = 0.0
f'(x) = 0.0000000160 +/- 0.0000000339
exact = 0.0000000000
$wasora diff.was f(x) = x^(3/2) x = 2.0 f'(x) = 2.1213203120 exact = 2.1213203436 x = 0.0 f'(x) = 0.0000000160 exact = 0.0000000000$ 

You may want to test what happens if the derivative at $$x=0$$ is computed without the optional fifth argument $$p$$.

## roots.was

This example uses the function solver gsl_root_fsolver_brent for Brent’s method to solve the following equation

!bt [ x^2 - 5 = 0 ] !et

with solution $$x = \sqrt{5} = 2.236068\dots$$. First, the function to be solved is defined in the following two files demo_fn.h and demo_fn.c:

@@@CODE 020-gsl/demo_fn.h

@@@CODE 020-gsl/demo_fn.c

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/roots.c

The corresponding wasora input:

# GSL example 33.10, root finding
a = 1
b = 0
c = -5
f(x) := (a*x + b)*x + c

PRINT %.7f root(f(x),x,0,5) sqrt(5)
$gcc -o roots roots.c -lgslcblas -lgsl -lm$ ./roots
using brent method
iter [    lower,     upper]      root        err  err(est)
1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
Converged:
6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666
$wasora roots.was 2.2360680 2.2360680$ 

As an illustration, the wasora version (which is significantly smaller than its corresponding C partner) uses the default tolerance for the root finding routines instead of the $$0.001$$ value of the original example, hence the difference in the results.

## min.was

The following program uses the Brent algorithm to find the minimum of the function $$f(x) = \cos(x) + 1$$, which occurs at $$x = \pi$$. The starting interval is $$[0,6]$$. The GSL version uses an initial guess for the minimum of $$x=2$$. However, wasora takes the initial guess at the center of the starting interval.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/min.c

The corresponding wasora input:

# GSL example 34.8, one-dimensional minimization
VAR x
PRINT %.7f func_min(cos(x)+1,x,0,6)
$gcc -o min min.c -lgslcblas -lgsl -lm$ ./min
using brent method
iter [    lower,     upper]       min        err  err(est)
0 [0.0000000, 6.0000000] 2.0000000 -1.1415927 6.0000000
1 [2.0000000, 6.0000000] 3.5278640 +0.3862713 4.0000000
2 [2.0000000, 3.5278640] 3.1748217 +0.0332290 1.5278640
3 [2.0000000, 3.1748217] 3.1264576 -0.0151351 1.1748217
4 [3.1264576, 3.1748217] 3.1414743 -0.0001183 0.0483641
5 [3.1414743, 3.1748217] 3.1415930 +0.0000004 0.0333474
Converged:
6 [3.1414743, 3.1415930] 3.1415927 +0.0000000 0.0001187
$wasora min.was 3.1415914$ 

As seen in the terminal, the functional func_min performs the minimization at once and does not provide a mean to show information about intermediate steps.

## rngunif.was

The following program demonstrates the use of a random number generator to produce uniform random numbers in the range $$[0.0, 1.0)$$. In the GSL version, The numbers depend on the seed used by the generator. The default seed can be changed with the GSL_RNG_SEED environment variable to produce a different stream of numbers. The generator itself can be changed using the environment variable GSL_RNG_TYPE. In wasora, the generator is fixed to the gsl_rng_knuthran2002 method, basically because the author looks up at Knuth as the most important figure of modern computer science. Therefore, to obtain two identical sequences of numbers both with GSL and wasora, we run both examples with the same generator and seed.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/rngunif.c

The corresponding wasora input:

# GSL example 18.13, random number generation
static_steps = 10
PRINT %.5f random(0,1,$1) $ gcc -o rngunif rngunif.c -lgslcblas -lgsl -lm
$GSL_RNG_SEED=123 GSL_RNG_TYPE=knuthran2002 ./rngunif GSL_RNG_TYPE=knuthran2002 GSL_RNG_SEED=123 0.32585 0.53345 0.27282 0.84846 0.64330 0.90798 0.19579 0.63717 0.92833 0.71629$ wasora rngunif.was 123
0.32585
0.53345
0.27282
0.84846
0.64330
0.90798
0.19579
0.63717
0.92833
0.71629
$ ## qrng.was Quasi-random sequences are used by wasora to perform parametric studies. The idea is to sweep a certain portion of the parameter space by choosing test points, which can be linear (type LINEAR), logarithmic (LOGARITHMIC), uniformly chosen at random (RANDOM), randomly distributed using a gaussian distribution (GAUSSIANRANDOM) or follow a certain quasi-random sequence (types SOBOL, NIEDERREITER, HALTON and REVERSEHALTON). By defining a two-dimensional window $$[0,1]\times[0,1]$$ and just printing the chosen parameters, the GSL example can be reproduced in wasora. The C program as implemented in the GSL documentation: @@@CODE 020-gsl/stat.c The corresponding wasora input: # GSL example 19.6 quasi-random sequences PARAMETRIC x y MIN 0 0 MAX 1 1 OUTER_STEPS 1024 TYPE sobol PRINT %.5f x y SEP " " $ gcc -o qrng qrng.c -lgslcblas -lgsl -lm
$./qrng > qrng-gsl.dat$ wasora qrng.was | tee qrng-was.dat | graph -Tpng -C -m -3 -S 2 > qrng.png
$diff -s qrng-gsl.dat qrng-was.dat Files qrng-gsl.dat and qrng-was.dat are identical$ 

## stat.was

A basic example of how to use the statistical functions.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/stat.c

The corresponding wasora input:

# GSL example 21.10, statistics
VECTOR data SIZE 5 DATA 17.2 18.1 16.5 18.3 12.6

N = vecsize(data)
mean = vecsum(data)/N
variance = sum((data(i) - mean)^2/(N-1), i, 1, N)

PRINT TEXT "the dataset is" %g data SEPARATOR ", "
PRINT TEXT "the sample mean is       " %g mean
PRINT TEXT "the estimated variance is" %g variance
PRINT TEXT "the largest value is     " %g vecmax(data)
PRINT TEXT "the smallest value is    " %g vecmin(data)
$gcc -o stat stat.c -lgslcblas -lgsl -lm$ ./stat
The dataset is 17.2, 18.1, 16.5, 18.3, 12.6
The sample mean is 16.54
The estimated variance is 5.373
The largest value is 18.3
The smallest value is 12.6
$wasora stat.was the dataset is, 17.2, 18.1, 16.5, 18.3, 12.6 the sample mean is 16.54 the estimated variance is 5.373 the largest value is 18.3 the smallest value is 12.6$ 

## ode-initval.was

The following program solves the second-order nonlinear Van der Pol oscillator equation

!bt [ u’‘(t) + u’(t) (u(t)^2 - 1) + u(t) = 0 ] !et

This can be converted into a first order system suitable for use with the routines described in this chapter by introducing a separate variable for the velocity, $$v = u'(t)$$,

!bt \begin{align*} u' &= v\cr v' &= -u + \mu v (1-u^2) \end{align*}

!et

The program begins by defining functions for these derivatives and their Jacobian. The main function uses driver level functions to solve the problem. The program evolves the solution from $$(u, v) = (1,0)$$ at $$t=0$$ to $$t=100$$. The step-size $$h$$ is automatically adjusted by the controller to maintain an absolute accuracy of $$10^{-6}$$ in the function values $$(u, v)$$.
The loop in the example prints the solution at the points $$t_i = 1, 2, \dots, 100$$.

It is important to note that wasora does not use GSL’s ordinary differential equations solvers because wasora works with DAEs instead of ODEs (which of course are a particular case of DAEs). The solver is based on the SUNDIALS library IDA. Nevertheless, dynamical systems are one of the main subjects of wasora so this example is particularly relevant. The special variable rel_error is set to $$10^{-6}$$ as in GSL and both the variables min_dt and max_dt are set to one, so as to reproduce the solution only at the points $$t_i=1,2,\dots,100$$.

The C program as implemented in the GSL documentation:

@@@CODE 020-gsl/ode-initval.c

The corresponding wasora input:

# GSL example 26.6, ordinary differential equations
# (solved with sundials, not with gsl)
PHASE_SPACE u v

# solver settings
rel_error = 1e-6
end_time = 100
min_dt = 1   # report results at t=0,1,2,...,100
max_dt = 1

# dynamical system parameter
mu = 10

# initial conditions
u_0 = 1
v_0 = 0

# dynamical system equations
u_dot .= v
v_dot .= -u + mu*v*(1-u^2)

PRINT t u v
$gcc -o ode-initval ode-initval.c -lgslcblas -lgsl -lm$ ./ode-initval > vanderpoolgsl.dat
$wasora ode-initval.was > vanderpoolwas.dat$ pyxplot vanderpool.ppl
$ ## interp.was The following program demonstrates the use of the interpolation and spline functions. It computes a cubic spline interpolation of the 10-point dataset $$(x_i, y_i)$$ where $$x_i = i + \sin(i)/2$$ and $$y_i = i + \cos(i^2)$$ for $$i = 0 \dots 9$$. To keep the original spirit, the data is plotted using GNU ploutils graph progam. The C program as implemented in the GSL documentation: @@@CODE 020-gsl/interp.c The corresponding wasora input: # GSL example 27.7.1, interpolation n = 10 VECTOR fx SIZE n VECTOR fy SIZE n FUNCTION yi(x) VECTORS fx fy INTERPOLATION cspline # remember the example has C indexes! fx(i) = (i-1) + sin(i-1)/2 fy(i) = (i-1) + cos((i-1)^2) # write the original data PRINT TEXT "\#m=-3,S=3" PRINT_FUNCTION yi FORMAT %g # print the interpolated data PRINT TEXT "\#m=3,S=0" PRINT_FUNCTION yi MIN fx(1) MAX fx(n) STEP (fx(n)-fx(1))/100 FORMAT %g $ gcc -o interp interp.c -lgslcblas -lgsl -lm
$./interp > interpgsl.dat$ head interpgsl.dat
#m=0,S=2
0 1
1.42074 1.5403
2.45465 1.34636
3.07056 2.08887
3.6216 3.04234
4.52054 5.9912
5.86029 5.87204
7.32849 7.30059
8.49468 8.39186
$tail interpgsl.dat 9.11 9.56698 9.12 9.58874 9.13 9.61052 9.14 9.63232 9.15 9.65414 9.16 9.67598 9.17 9.69783 9.18 9.71969 9.19 9.74156 9.2 9.76343$ graph -Tpng -C -L GSL < interpgsl.dat > interpgsl.png
$wasora interp.was > interpwas.dat$ head interpwas.dat
#m=-3,S=3
0   1
1.42074 1.5403
2.45465 1.34636
3.07056 2.08887
3.6216  3.04234
4.52054 5.9912
5.86029 5.87204
7.32849 7.30059
8.49468 8.39186
$tail interpwas.dat 8.37751 8.23336 8.46957 8.35571 8.56164 8.49437 8.6537 8.6486 8.74576 8.8162 8.83782 8.99493 8.92988 9.18258 9.02194 9.3769 9.114 9.57568 9.20606 9.77669$ graph -Tpng -C -L wasora < interpwas.dat > interpwas.png
$ ## interpp.was The next program demonstrates a periodic cubic spline with 4 data points. Note that the first and last points must be supplied with the same $$y$$-value for a periodic spline. The C program as implemented in the GSL documentation: @@@CODE 020-gsl/interpp.c The corresponding wasora input: # GSL example 27.7.2, interpolation FUNCTION f(x) INTERPOLATION cspline_periodic DATA { 0.00 0.15 0.10 0.70 0.27 -0.10 0.30 0.15 } # original data PRINT TEXT "\#m=-3,S=3" PRINT_FUNCTION f FORMAT %g # interpolated data PRINT TEXT "\#m=3,S=0" PRINT_FUNCTION f MIN f_a MAX f_b STEP (f_b-f_a)/100 FORMAT %g $ gcc -o interpp interpp.c -lgslcblas -lgsl -lm
$./interpp > interppgsl.dat$ graph -Tpng -C -L GSL < interppgsl.dat > interppgsl.png
$wasora interpp.was > interppwas.dat$ graph -Tpng -C -L wasora < interppwas.dat > interppwas.png
$ ## multimin.was This example program finds the minimum of a paraboloid function. The location of the minimum is offset from the origin in $$x$$ and $$y$$, and the function value at the minimum is non-zero. The initial step-size is chosen as 0.01, a conservative estimate in this case, and the line minimization parameter is set at 0.0001. The program terminates when the norm of the gradient has been reduced below 0.001. In wasora, as the minimization function is algebraic, the gradient can be entered easily as shown. If it was not entered, it would be automatically computed by finite differences (new versions of GSL do the same). The C program as implemented in the GSL documentation: @@@CODE 020-gsl/multimin.c The corresponding wasora input: # GSL example 36.9.1, multidimensional minimization a = 10 b = 20 c = 30 x0 = 1 y0 = 2 f(x,y) := a*(x-x0)^2 + b*(y-y0)^2 + c MINIMIZE f ALGORITHM conjugate_fr { GUESS 5 7 GRADIENT 2*a*(x-x0) 2*b*(y-y0) STEP 1e-2 GRADTOL 1e-3 TOL 1e-4 } IF done_outer PRINT "\# minimum found at" ENDIF IF step_inner=1|done_outer PRINT %5.0f step_outer %.5f x y %10.5f f(x,y) SEPARATOR " " ENDIF $ gcc -o multimin multimin.c -lgslcblas -lgsl -lm
$./multimin 1 4.99629 6.99072 687.84780 2 4.98886 6.97215 683.55456 3 4.97400 6.93501 675.01278 4 4.94429 6.86073 658.10798 5 4.88487 6.71217 625.01340 6 4.76602 6.41506 561.68440 7 4.52833 5.82083 446.46694 8 4.05295 4.63238 261.79422 9 3.10219 2.25548 75.49762 10 2.85185 1.62963 67.03704 11 2.19088 1.76182 45.31640 12 0.86892 2.02622 30.18555 Minimum found at: 13 1.00000 2.00000 30.00000$ wasora multimin.was
0 0.00000 0.00000  120.00000
1 4.99629 6.99072  687.84780
2 4.98886 6.97215  683.55456
3 4.97400 6.93501  675.01278
4 4.94429 6.86073  658.10798
5 4.88487 6.71217  625.01340
6 4.76602 6.41506  561.68440
7 4.52833 5.82083  446.46694
8 4.05295 4.63238  261.79422
9 3.10219 2.25548   75.49762
10 1.20067 -2.49832  435.09973
11 2.19088 1.76182   45.31640
12 0.86892 2.02622   30.18555
13 3.51283 1.49743   98.19448
# minimum found at
13 1.00000 2.00000   30.00000
$ This example shows an important difference. The GSL program prints a line only after successful steps of the optimization method. But wasora prints a line with the instantaneous position each time the objective function $$f(x,y)$$ has to be computed, as the whole input is executed to evaluate the function to be minimized (this example is simple, but wasora was designed to optimize complex computations such as nuclear reactor designs). Therefore, GSL shadows the real number of evaluations (which are directly related to the computational cost of the problem). Analyze with care. ## nmsimplex.was Here is another example using the Nelder-Mead Simplex algorithm to minimize the same example object function, as above. The C program as implemented in the GSL documentation: @@@CODE 020-gsl/nmsimplex.c The corresponding wasora input: # GSL example 36.9.2, multidimensional minimization a = 10 b = 20 c = 30 x0 = 1 y0 = 2 f(x,y) := a*(x-x0)^2 + b*(y-y0)^2 + c MINIMIZE f ALGORITHM nmsimplex2 { GUESS 5 7 GRADIENT 2*a*(x-x0) 2*b*(y-y0) STEP 1 1 TOL 1e-2 } IF in_outer_initial PRINT TEXT "\# x y z" ENDIF IF done_outer PRINT TEXT "\# converged to minimum at" ENDIF IF (not(in_outer_initial)&step_inner=1)|done_outer PRINT %5.0f step_outer %10.3e x y TEXT "f() =" %7.3f f(x,y) SEPARATOR " " ENDIF $ gcc -o nmsimplex nmsimplex.c -lgslcblas -lgsl -lm
$./nmsimplex 1 6.500e+00 5.000e+00 f() = 512.500 size = 1.130 2 5.250e+00 4.000e+00 f() = 290.625 size = 1.409 3 5.250e+00 4.000e+00 f() = 290.625 size = 1.409 4 5.500e+00 1.000e+00 f() = 252.500 size = 1.409 5 2.625e+00 3.500e+00 f() = 101.406 size = 1.847 6 2.625e+00 3.500e+00 f() = 101.406 size = 1.847 7 -1.776e-15 3.000e+00 f() = 60.000 size = 1.847 8 2.094e+00 1.875e+00 f() = 42.275 size = 1.321 9 2.578e-01 1.906e+00 f() = 35.684 size = 1.069 10 5.879e-01 2.445e+00 f() = 35.664 size = 0.841 11 1.258e+00 2.025e+00 f() = 30.680 size = 0.476 12 1.258e+00 2.025e+00 f() = 30.680 size = 0.367 13 1.093e+00 1.849e+00 f() = 30.539 size = 0.300 14 8.830e-01 2.004e+00 f() = 30.137 size = 0.172 15 8.830e-01 2.004e+00 f() = 30.137 size = 0.126 16 9.582e-01 2.060e+00 f() = 30.090 size = 0.106 17 1.022e+00 2.004e+00 f() = 30.005 size = 0.063 18 1.022e+00 2.004e+00 f() = 30.005 size = 0.043 19 1.022e+00 2.004e+00 f() = 30.005 size = 0.043 20 1.022e+00 2.004e+00 f() = 30.005 size = 0.027 21 1.022e+00 2.004e+00 f() = 30.005 size = 0.022 22 9.920e-01 1.997e+00 f() = 30.001 size = 0.016 23 9.920e-01 1.997e+00 f() = 30.001 size = 0.013 converged to minimum at 24 9.920e-01 1.997e+00 f() = 30.001 size = 0.008$ wasora nmsimplex.was
#    x        y         z
1  6.000e+00  6.000e+00 f() = 600.000
2  5.500e+00  5.000e+00 f() = 412.500
3  6.750e+00  2.000e+00 f() = 360.625
4  5.500e+00  1.000e+00 f() = 252.500
5  4.000e+00  3.000e+00 f() = 140.000
6  2.875e+00  5.000e-01 f() = 110.156
7 -1.776e-15  3.000e+00 f() =  60.000
8 -2.500e-01  6.000e+00 f() = 365.625
9 -5.313e-01  1.375e+00 f() =  61.260
10  2.352e+00  7.812e-01 f() =  77.974
11 -1.248e+00  2.477e+00 f() =  85.079
12  1.588e+00  2.564e+00 f() =  39.834
13  1.261e+00  1.651e+00 f() =  33.118
14  1.760e+00  1.804e+00 f() =  36.551
15  7.173e-01  1.828e+00 f() =  31.390
16  9.134e-01  2.131e+00 f() =  30.417
17  7.181e-01  2.088e+00 f() =  30.951
18  1.097e+00  2.060e+00 f() =  30.167
19  1.000e+00  1.962e+00 f() =  30.029
20  1.085e+00  1.948e+00 f() =  30.127
21  9.954e-01  2.043e+00 f() =  30.037
22  1.047e+00  1.986e+00 f() =  30.026
23  1.015e+00  2.019e+00 f() =  30.009
24  9.732e-01  1.984e+00 f() =  30.012
# converged to minimum at
24  9.920e-01  1.997e+00 f() =  30.001
$ ## nlfit.was The following example program fits a weighted exponential model with background to experimental (simulated) data , $$Y = A \exp(-\lambda t) + b$$. The wasora implementation, instead of generating the random gaussian data (which can be nevertheless implemented with the random_gauss function), the data generated by the GSL program is sent to a file which is used to define the objective function. This situation on the one hand is closer to actual applications of wasora to fit real experimental data and, on the other hand, we make sure we use the same set of data in both cases so we expect to find the same results. The C program as implemented in the GSL documentation: @@@CODE 020-gsl/nlfit.c The corresponding wasora input: # GSL example 38.9, non-linear fit # # GRADIENT is optional, wasora will compute it numerically if not provided # VERBOSE activates step-by-step information # # initial guess A = 0 lambda = 1 b = 0 Y(t) := A * exp(-lambda*t) + b FUNCTION data(t) FILE_PATH experimental.dat COLUMNS 2 3 FIT Y TO data VIA A lambda b { DELTAEPSREL 1e-4 DELTAEPSABS 1e-4 # GRADIENT exp(-lambda*t) -t*A*exp(-lambda*t) 1 VERBOSE } IF done_outer PRINT %.5f TEXT "A =" A TEXT "+/-" sigma_A SEP " " PRINT %.5f TEXT "lambda =" lambda TEXT "+/-" sigma_lambda SEP " " PRINT %.5f TEXT "b =" b TEXT "+/-" sigma_b SEP " " ENDIF IF done_outer PRINT_FUNCTION FILE_PATH nlfit.dat data Y ENDIF $ gcc -o nlfit nlfit.c -lgslcblas -lgsl -lm
$./nlfit | tee gsl.out data: 0 6.01339 0.1 data: 1 5.51538 0.1 data: 2 5.26109 0.1 data: 3 4.77746 0.1 data: 4 4.45135 0.1 data: 5 3.9049 0.1 data: 6 3.50439 0.1 data: 7 3.415 0.1 data: 8 3.24274 0.1 data: 9 3.1222 0.1 data: 10 2.83763 0.1 data: 11 2.5347 0.1 data: 12 2.43917 0.1 data: 13 2.38083 0.1 data: 14 2.31609 0.1 data: 15 2.06083 0.1 data: 16 1.94568 0.1 data: 17 1.91413 0.1 data: 18 1.75951 0.1 data: 19 1.66507 0.1 data: 20 1.73793 0.1 data: 21 1.57552 0.1 data: 22 1.52507 0.1 data: 23 1.40961 0.1 data: 24 1.39521 0.1 data: 25 1.41689 0.1 data: 26 1.37604 0.1 data: 27 1.26095 0.1 data: 28 1.28963 0.1 data: 29 1.42267 0.1 data: 30 1.22829 0.1 data: 31 1.19918 0.1 data: 32 1.18999 0.1 data: 33 0.930076 0.1 data: 34 1.22461 0.1 data: 35 1.14738 0.1 data: 36 1.114 0.1 data: 37 1.19512 0.1 data: 38 1.26958 0.1 data: 39 1.06198 0.1 iter: 0 x = 1.00000000 0.00000000 0.00000000 |f(x)| = 117.349 status = success iter: 1 x = 1.64656260 0.01814555 0.64656260 |f(x)| = 76.4594 status = success iter: 2 x = 2.85864925 0.08091443 1.44791355 |f(x)| = 37.6859 status = success iter: 3 x = 4.94892814 0.11943460 1.09463075 |f(x)| = 9.58212 status = success iter: 4 x = 5.02174105 0.10287695 1.03389255 |f(x)| = 5.63076 status = success iter: 5 x = 5.04520407 0.10405524 1.01941629 |f(x)| = 5.44398 status = success iter: 6 x = 5.04535782 0.10404906 1.01924871 |f(x)| = 5.44397 chisq/dof = 0.800996 A = 5.04536 +/- 0.06028 lambda = 0.10405 +/- 0.00316 b = 1.01925 +/- 0.03782 status = success$ cat gsl.out | grep data > experimental.dat
$wasora nlfit.was fititeration: 0 0.000e+00 1.000e+00 0.000e+00 status = success |r|^2 = 16.7263 fititeration: 1 5.040e+00 1.000e+00 2.076e+00 status = success |r|^2 = 6.6945 fititeration: 2 5.033e+00 8.766e-01 2.075e+00 status = success |r|^2 = 6.49716 fititeration: 3 5.021e+00 6.777e-01 2.070e+00 status = success |r|^2 = 6.05048 fititeration: 4 4.986e+00 4.101e-01 2.043e+00 status = success |r|^2 = 4.95003 fititeration: 5 4.686e+00 1.104e-01 1.872e+00 status = success |r|^2 = 4.53229 fititeration: 6 4.805e+00 1.449e-01 1.557e+00 status = success |r|^2 = 1.80061 fititeration: 7 4.922e+00 9.566e-02 1.091e+00 status = success |r|^2 = 1.06197 fititeration: 8 5.036e+00 1.041e-01 1.027e+00 status = success |r|^2 = 0.545438 fititeration: 9 5.045e+00 1.040e-01 1.019e+00 status = success |r|^2 = 0.544395 fititeration: 10 5.045e+00 1.040e-01 1.019e+00 status = success |r|^2 = 0.544395 A = 5.04536 +/- 0.60279 lambda = 0.10405 +/- 0.03157 b = 1.01925 +/- 0.37824 # chisq/dof = 0.00800988 # A = 5.0454 # +/- 0.60279 # lambda = 0.10405 # +/- 0.031571 # b = 1.0192 # +/- 0.37824$ pyxplot nlfit.ppl
\$ 

Note that wasora cannot currently handle uncertainties on the experimental data, so they are assumed to be equal to one (instead of the GSL example where they are identically equal to 0.1 for every point). Therefore, the reported uncertainties on the coefficients are larger in wasora than in the GSL example.