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
$
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 = r*x*(1-x) # apply the map
x
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
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
.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.#!/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,
# 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
= 20 # reference temperature
T0 = vecmin(vec_A_T) # smallest argument of function A(T)
T_min = vecmax(vec_A_T) # largest argument of function A(T)
T_max
# 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.