`PRINT "Hello $1!"`

```
$ feenox hello.fee World
Hello World!
$ feenox hello.fee Universe
Hello Universe!
$
```

Internal symbol

`pi`

equal to libc’s`M_PI`

\pi = \pi

Four times the arc tangent of one

\pi = 4 \cdot \arctan{1}

The root of \tan x=0

x / \tan x = 0

The integral of the gaussian bell, squared

\pi = \left[ \int_{-\infty}^{\infty} e^{-x^2} \, dx \right]^2

The integral of (x^2 + y^2) < 1 inside a unit square

\pi = \int_{-1}^{+1} \int_{-1}^{+1} \left[(x^2 + y^2) < 1\right] \, dy \, dx

The integral of one inside a parametrized circle

\pi = \int_{-1}^{+1} \int_{-\sqrt{1-x^2}}^{+\sqrt{1-x^2}} \, dy \, dx

The Gregoy-Leibniz sum

\pi = 4 \cdot \sum_{i}^{\infty} \frac{(-1)^{i+1}}{2i-1}

The Abraham-Sharp sum

\pi = \sum_{i}^{\infty} \frac{2 \cdot (-1)^i \cdot 3^{\frac{1}{2-i}}}{2i+1}

An integral which is equal to 22/7-\pi

22/7 - \pi = \int_{0}^{1} \frac{x^4 \cdot (1-x)^4}{1+x^2} \, dx

Ramanujan-Sato series

\frac{1}{\pi} = \frac{2 \cdot \sqrt{2}}{99^2} \cdot \sum_{i=0}^{\infty} \frac{(4i)!}{(i^4)!} \cdot \frac{26390 \cdot i + 1103}{396^{4i}}

```
# computing pi in ten different ways
VECTOR piapprox SIZE 10
VAR x y
# the double-precision internal representation of pi (M_PI)
1] = pi
piapprox[
# four times the arc-tangent of the unity
2] = 4*atan(1)
piapprox[
# the abscissae x where tan(x) = 0 in the range [3:3.5]
3] = root(tan(x), x, 3, 3.5)
piapprox[
# the square of the numerical integral of the guassian curve
4] = integral(exp(-x^2), x, -10, 10)^2
piapprox[
# the numerical integral of a circle inscribed in a unit square
5] = integral(integral((x^2+y^2) < 1, y, -1, 1), x, -1, 1)
piapprox[
# the numerical integral of a circle parametrized with sqrt(1-x^2)
6] = integral(integral(1, y, -sqrt(1-x^2), sqrt(1-x^2)), x, -1, 1)
piapprox[
# the gregory-leibniz sum (one hundred thousand terms)
7] = 4*sum((-1)^(i+1)/(2*i-1), i, 1, 1e5)
piapprox[
# the abraham sharp sum (twenty-one terms)
8] = sum(2*(-1)^i * 3^(1/2-i)/(2*i+1), i, 0, 20)
piapprox[
# this integral is equal to 22/7-pi
9] = 22/7-integral((x^4*(1-x)^4)/(1+x^2), x, 0, 1)
piapprox[
# ramanujan-sato series
10] = 1/(2*sqrt(2)/(99^2)*sum(gammaf(4*i)/gammaf(i)^4 * (26390*i + 1103)/(396^(4*i)), i, 0, 2))
piapprox[
PRINT_VECTOR "% .16f" piapprox piapprox(i)-pi
```

```
$ feenox pi.fee
3.1415926535897931 0.0000000000000000
3.1415926535897931 0.0000000000000000
3.1415926535897936 0.0000000000000004
3.1415926535899108 0.0000000000001177
3.1417605332857996 0.0001678796960065
3.1415956548678512 0.0000030012780581
3.1415826535897198 -0.0000100000000733
3.1415926535956351 0.0000000000058420
3.1415926535897931 0.0000000000000000
3.1415927109074255 0.0000000573176324
$
```

You live in a country with a high inflation rate. A retailer offers you a good with two payment options:

- 30% off the full price plus 3 equal monthly installments, or
- full price in 12 equal monthly installments

If the inflation rate is small (large), option a (b) wins.

**Question** What is the monthly inflation rate at which
the two options are equally good (or bad)?

```
# n = number of payments
# d = discount
# r = monthly inflation rate
(n, d, r) = sum((1-d)/n * (1/(1+r))^(i-1), i, 1, n)
present_value= root(present_value(12,0,r) - present_value(3,0.3,r), r, 0, 0.2)
r0
PRINTF "the tipping monthly inflation is %g (that means %.1f%% per month and %.1f%% per year)" r0 100*r0 100*(1+r0)^12
PRINT
PRINTF "with 12 payments no discount with r = %g you pay %.3f" r0+0.01 present_value(12,0,r0+0.01)
PRINTF "with 3 payments and 30%% off with r = %g you pay %.3f" r0+0.01 present_value(3,0.3,r0+0.01)
PRINT
PRINTF "with 12 payments no discount with r = %g you pay %.3f" r0-0.01 present_value(12,0,r0-0.01)
PRINTF "with 3 payments and 30%% off with r = %g you pay %.3f" r0-0.01 present_value(3,0.3,r0-0.01)
```

```
$ feenox inflation.fee
the tipping monthly inflation is 0.0931703 (that means 9.3% per month and 291.2% per year)
with 12 payments no discount with r = 0.10317 you pay 0.617
with 3 payments and 30% off with r = 0.10317 you pay 0.637
with 12 payments no discount with r = 0.0831703 you pay 0.669
with 3 payments and 30% off with r = 0.0831703 you pay 0.648
$
```

Plot the asymptotic behavior of the logistic map

x_{n+1} = r \cdot x \cdot (1-x)

for a range of the parameter r.

```
DEFAULT_ARGUMENT_VALUE 1 2.6 # by default show r in [2.6:4]
DEFAULT_ARGUMENT_VALUE 2 4
= 2^10
steps_per_r = 2^8
steps_asymptotic = 2^10
steps_for_r
static_steps = steps_for_r*steps_per_r
# change r every steps_per_r steps
IF mod(step_static,steps_per_r)=1
= quasi_random($1,$2)
r ENDIF
= 0.5 # start at x = 0.5
x_init x = r*x*(1-x) # apply the map
IF step_static-steps_per_r*floor(step_static/steps_per_r)>(steps_per_r-steps_asymptotic)
# write the asymptotic behavior only
PRINT r x
ENDIF
```

```
$ gnuplot
gnuplot> plot "< feenox logistic.fee" w p pt 50 ps 0.02
gnuplot> quit
$
```

When directly executing FeenoX, one gives a single argument to the executable with the path to the main input file. For example, the following input computes the first twenty numbers of the Fibonacci sequence using the closed-form formula

f(n) = \frac{\varphi^n - (1-\varphi)^n}{\sqrt{5}}

where \varphi=(1+\sqrt{5})/2 is the Golden ratio.

```
# the fibonacci sequence as function
= (1+sqrt(5))/2
phi (n) = (phi^n - (1-phi)^n)/sqrt(5)
fPRINT_FUNCTION f MIN 1 MAX 20 STEP 1
```

```
$ feenox fibo_formula.fee | tee one
1 1
2 1
3 2
4 3
5 5
6 8
7 13
8 21
9 34
10 55
11 89
12 144
13 233
14 377
15 610
16 987
17 1597
18 2584
19 4181
20 6765
$
```

We could also have computed these twenty numbers by using the direct definition of the sequence into a vector \vec{f} of size 20.

```
# the fibonacci sequence as a vector
VECTOR f SIZE 20
i]<1:2> = 1
f[i]<3:vecsize(f)> = f[i-2] + f[i-1]
f[
PRINT_VECTOR i f
```

```
$ feenox fibo_vector.fee > two
$
```

Finally, we print the sequence as an iterative problem and check that the three outputs are the same.

```
static_steps = 20
#static_iterations = 1476 # limit of doubles
IF step_static=1|step_static=2
= 1
f_n = 1
f_nminus1 = 1
f_nminus2 ELSE
= f_nminus1 + f_nminus2
f_n = f_nminus1
f_nminus2 = f_n
f_nminus1 ENDIF
PRINT step_static f_n
```

```
$ feenox fibo_iterative.fee > three
$ diff one two
$ diff two three
$
```

This example illustrates how well FeenoX integrates into the Unix
philosophy. Let’s say one has a function f(t) as an ASCII file with two columns and
one wants to compute the derivative f'(t). Just pipe the function file into
this example’s input file `derivative.fee`

used as a
filter.

For example, this small input file `f.fee`

writes the
function of time provided in the first command-line argument from zero
up to the second command-line argument:

```
end_time = $2
PRINT t $1
```

```
$ feenox f.fee "sin(t)" 1
0 0
0.0625 0.0624593
0.125 0.124675
0.1875 0.186403
0.25 0.247404
0.3125 0.307439
0.375 0.366273
0.4375 0.423676
0.5 0.479426
0.5625 0.533303
0.625 0.585097
0.6875 0.634607
0.75 0.681639
0.8125 0.726009
0.875 0.767544
0.9375 0.806081
1 0.841471
$
```

Then we can pipe the output of this command to the derivative filter. Note that

- The
`derivative.fee`

has the execution flag has on and a shebang line pointing to a global location of the FeenoX binary in`/usr/local/bin`

e.g. after doing`sudo make install`

. - The first argument of
`derivative.fee`

controls the time step. This is only important to control the number of output lines. It does not have anything to do with precision, since the derivative is computed using an adaptive centered numerical differentiation scheme using the GNU Scientific Library. - Before doing the actual differentiation, the input data is interpolated using a third-order monotonous scheme (also with GSL).
- TL;DR: this is not just “current value minus last value divided time increment.”

```
#!/usr/local/bin/feenox
# read the function from stdin
FUNCTION f(t) FILE - INTERPOLATION steffen
# detect the domain range
= vecmin(vec_f_t)
a = vecmax(vec_f_t)
b
# time step from arguments (or default 10 steps)
DEFAULT_ARGUMENT_VALUE 1 (b-a)/10
= $1
h
# compute the derivative with a wrapper for gsl_deriv_central()
VAR t'
(t) = derivative(f(t'),t',t)
f'
# write the result
PRINT_FUNCTION f' MIN a+0.5*h MAX b-0.5*h STEP h
```

```
$ chmod +x derivative.sh
$ feenox f.fee "sin(t)" 1 | ./derivative.fee 0.1 | tee f_prime.dat
0.05 0.998725
0.15 0.989041
0.25 0.968288
0.35 0.939643
0.45 0.900427
0.55 0.852504
0.65 0.796311
0.75 0.731216
0.85 0.66018
0.95 0.574296
$
```

When solving thermal-mechanical problems it is customary to use thermal expansion coefficients in order to take into account the mechanical strains induced by changes in the material temperature with respect to a reference temperature where the deformation is identically zero. These coefficients \alpha are defined as the partial derivative of the strain \epsilon with respect to temperature T such that differential relationships like

d\epsilon = \frac{\partial \epsilon}{\partial T} \, dT = \alpha \cdot dT

hold. This derivative \alpha is
called the *instantaneous* thermal expansion coefficient. For
finite temperature increments, one would like to be able to write

\Delta \epsilon = \alpha \cdot \Delta T

But if the strain is not linear with respect to the temperature—which
is the most common case—then \alpha
depends on T. Therefore, when dealing
with finite temperature increments \Delta T =
T-T_0 where the thermal expansion coefficient \alpha(T) depends on the temperature itself
then *mean* values for the thermal expansion ought to be
used:

\Delta \epsilon = \int_{\epsilon_0}^{\epsilon} d\epsilon^{\prime} = \int_{T_0}^{T} \frac{\partial \epsilon}{\partial T^\prime} \, dT^\prime = \int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime

We can multiply and divide by \Delta T to obtain

\int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime \cdot \frac{\Delta T}{\Delta T} = \bar{\alpha}(T) \cdot \Delta T

where the mean expansion coefficient for the temperature range [T_0,T] is

\bar{\alpha}(T) = \frac{\displaystyle \int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime}{T-T_0}

This is of course the classical calculus result of the mean value of a continuous one-dimensional function in a certain range.

Let \epsilon(T) be the linear thermal expansion of a given material in a certain direction when heating a piece of such material from an initial temperature T_0 up to T so that \epsilon(T_0)=0.

From our previous analysis, we can see that in fig. 1:

\begin{aligned} A(T) &= \alpha(T) = \frac{\partial \epsilon}{\partial T} \\ B(T) &= \bar{\alpha}(T) = \frac{\epsilon(T)}{T-T_0} = \frac{\displaystyle \int_{T_0}^{T} \alpha(T^\prime) \, dT^\prime}{T - T_0} \\ C(T) &= \epsilon(T) = \int_{T_0}^T \alpha(T^\prime) \, dT^\prime \end{aligned}

Therefore,

- A(T) can be computed out of C(T)
- B(T) can be computed either out of A(T) or C(T)
- C(T) can be computed out of A(T)

```
# just in case we wanted to interpolate with another method (linear, splines, etc.)
DEFAULT_ARGUMENT_VALUE 1 steffen
# read columns from data file and interpolate
# A is the instantaenous coefficient of thermal expansion x 10^-6 (mm/mm/ºC)
FUNCTION A(T) FILE asme-expansion-table.dat COLUMNS 1 2 INTERPOLATION $1
# B is the mean coefficient of thermal expansion x 10^-6 (mm/mm/ºC) in going
# from 20ºC to indicated temperature
FUNCTION B(T) FILE asme-expansion-table.dat COLUMNS 1 3 INTERPOLATION $1
# C is the linear thermal expansion (mm/m) in going from 20ºC
# to indicated temperature
FUNCTION C(T) FILE asme-expansion-table.dat COLUMNS 1 4 INTERPOLATION $1
VAR T' # dummy variable for integration
T0 = 20 # reference temperature
T_min = vecmin(vec_A_T) # smallest argument of function A(T)
T_max = vecmax(vec_A_T) # largest argument of function A(T)
# compute one column from another one
(T) := 1e3*derivative(C(T'), T', T)
A_fromC
(T) := integral(A(T'), T', T0, T)/(T-T0)
B_fromA(T) := 1e3*C(T)/(T-T0) # C is in mm/m, hence the 1e3
B_fromC
(T) := 1e-3*integral(A(T'), T', T0, T)
C_fromA
# write interpolated results
PRINT_FUNCTION A A_fromC B B_fromA B_fromC C C_fromA MIN T_min+1 MAX T_max-1 STEP 1
```

```
$ cat asme-expansion-table.dat
# temp A B C
20 21.7 21.7 0
50 23.3 22.6 0.7
75 23.9 23.1 1.3
100 24.3 23.4 1.9
125 24.7 23.7 2.5
150 25.2 23.9 3.1
175 25.7 24.2 3.7
200 26.4 24.4 4.4
225 27.0 24.7 5.1
250 27.5 25.0 5.7
275 27.7 25.2 6.4
300 27.6 25.5 7.1
325 27.1 25.6 7.8
$ feenox asme-expansion.fee > asme-expansion-interpolation.dat
$ pyxplot asme-expansion.ppl
$
```

The conclusion (see this, this and this reports) is that values rounded to only one decimal value as presented in the ASME code section II subsection D tables are not enough to satisfy the mathematical relationships between the physical magnitudes related to thermal expansion properties of the materials listed. Therefore, care has to be taken as which of the three columns is chosen when using the data for actual thermo-mechanical numerical computations. As an exercise, the reader is encouraged to try different interpolation algorithms to see how the results change.

Spoiler alert: they are also highly sensible to the interpolation method used to “fill in” the gaps between the table values.