Two-dimensional functions

  • Difficulty: 010/100
  • Author: jeremy theler
  • Keywords: FUNCTION, :=, FILE_PATH, DATA, INTERPOLATION, nearest, rectangle, PRINT_FUNCTION, MIN, MAX, STEP, NUMBER, FILE, OUTPUT_FILE,

These examples illustrate the facilities wasora provides to interpolate two-dimensional functions. As shown in section ref{007-functions}, multidimensional functions may be defined using algebraic expressions, inline data, external files, vectors or dynamically-loaded routines. The main focus of this section is to illustrate the difference between algebraic and pointwise-defined two-dimensional functions, and to further illustrate the difference between a nearest-neighbor interpolation (based on a \(k\)-dimensional tree structure) and a rectangle interpolation (based on finite-elements-like shape functions for quadrangles).

paraboloid.was

This input defines an algebraic function \(f(x,y)\) and uses PRINT_FUNCTION to dump its contents (as three columns containing \(x\), \(y\) and \(f(x,y)\) as shown in the terminal mimic) to the standard output. The range is mandatory because \(f(x,y)\) is defined by an algebraic expression.

As far as the function \(f(x,y)\) is concerned, \(a\) and \(b\) are taken as parameters. Even though they can change over time, the value they take is the value they have when \(f(x,y)\) is evaluated. In this case, \(f(x,y)\) is evaluated when the instruction PRINT_FUNCTION is executed, so the output is written with \(b=2\).

a = 1
b = 1

f(x,y) := (x/a)^2 + (y/b)^2

b = 2

# range is mandatory as f(x,y) is algebraically-defined
PRINT_FUNCTION f MIN -b -b MAX b b STEP b/20 b/20
$ wasora paraboloid.was > paraboloid.dat
$ head paraboloid.dat
-2.000000e+00   -2.000000e+00   5.000000e+00
-2.000000e+00   -1.900000e+00   4.902500e+00
-2.000000e+00   -1.800000e+00   4.810000e+00
-2.000000e+00   -1.700000e+00   4.722500e+00
-2.000000e+00   -1.600000e+00   4.640000e+00
-2.000000e+00   -1.500000e+00   4.562500e+00
-2.000000e+00   -1.400000e+00   4.490000e+00
-2.000000e+00   -1.300000e+00   4.422500e+00
-2.000000e+00   -1.200000e+00   4.360000e+00
-2.000000e+00   -1.100000e+00   4.302500e+00
$ gnuplot paraboloid.gp
$ 
paraboloid.png
paraboloid.png
paraboloid2d.png
paraboloid2d.png

nearest.was

This input defines a two-dimensional pointwise-defined function \(g(x,y)\) using data inlined in the input file using the DATA keyword. As can be seen, expressions are allowed. However, they are evaluated at parse-time so references to variables should be avoided, as they will default to zero. For that end, numbers defined with the NUMBER keyword should be used. By default, pointwise-defined multidimensional functions are interpolated by the nearest neighbor algorithm (i.e. default is INTERPOLATION nearest). When calling PRINT_FUNCTION with no range, the definition points are printed. When a range is given, the function gets evaluated at the grid points.

In this case, first the function \(g(x,y)\) is dumped into a file called g_def.dat with no range. Then, the function is dumped into a file called g_int.dat using a range and thus resulting in an interpolated output using nearest neighbors.

FUNCTION g(x,y) DATA {
0    0    1-1
0    1    1-0.5
0    2    1
1    0    1
1    1    1+0.25
1    2    1
2    0    1-0.25
2    1    1+0.25
2    2    1+0.5
}

# print g(x,y) at the definition points
PRINT_FUNCTION g FILE_PATH g_def.dat

# print g(x,y) at the selected range
# by defaults wasora interpolates using nearest neighbors
PRINT_FUNCTION g FILE_PATH g_int.dat MIN 0 0 MAX 2 2 STEP 0.05 0.05
$ wasora nearest.was
$ gnuplot nearest.gp
$ 
nearest.png
nearest.png
nearest2d.png
nearest2d.png

rectangle.was

This time, a function \(h(x,y)\) is defined by reading point-wise data from a file. This file is g_def.dat which are the definition points of the function \(g(x,y)\) of the previous example. Note that when reading function data from a file, no expressions are allowed. Function \(h(x,y)\) is interpolated using the rectangle method. The interpolated data is written in a file called h_int.dat, in which the function \(h(x,y)\) is evaluated at the very same points the function \(g(x,y)\) of the previous example was.

# h(x,y) is equal to g(x,y) at the definition points,
# but it is interpolated differently
FUNCTION h(x,y) INTERPOLATION rectangle FILE_PATH g_def.dat

# print h(x,y) at the selected range
PRINT_FUNCTION h FILE_PATH h_int.dat MIN 0 0 MAX 2 2 STEP 0.05 0.05
$ cat g_def.dat
0.000000e+00    0.000000e+00    0.000000e+00
0.000000e+00    1.000000e+00    5.000000e-01
0.000000e+00    2.000000e+00    1.000000e+00

1.000000e+00    0.000000e+00    1.000000e+00
1.000000e+00    1.000000e+00    1.250000e+00
1.000000e+00    2.000000e+00    1.000000e+00

2.000000e+00    0.000000e+00    7.500000e-01
2.000000e+00    1.000000e+00    1.250000e+00
2.000000e+00    2.000000e+00    1.500000e+00

$ wasora rectangle.was
$ gnuplot rectangle.gp
$ 
rectangle.png
rectangle.png
rectangle2d.png
rectangle2d.png

scattered.was

The last example of two-dimensional interpolation involves a pointwise-defined function~\(s(x,y)\) using scattered data, i.e. not necessarily over a rectangular grid. In this case—and with no information about any underlying finite-element-like mesh—wasora can use either a nearest-neighbor interpolation or a Shepard-like inverse distance weighting. For the original Shepard method, the only parameter than can be tweaked is the exponent~\(p\) of the distance in the weight~\(w_i=1/d_i^p\). For the modified Shepard algorithm, the radius~\(r\) of the nearest neighbors taken into account is to be provided. These nearest neighbors are found using a \(k\)-dimensional tree, that is a very efficient way of doing this task. For complex functions, all the alternatives should be investigated taking into account accuracy and code speed.

# scattered multidimensional data may be interpolated
# using a nearest-neighbor approach
FUNCTION n(x,y) INTERPOLATION nearest DATA {
0     0    0
1     0    1
0     1    2
-0.5  0.5  3
-1    -1   2
0.75  0    1.5
0.25  0.25 1
}

# another way of giving the same set of data
VECTOR datax SIZE 7 DATA 0 1 0 -0.5 -1 0.75 0.25
VECTOR datay SIZE 7 DATA 0 0 1  0.5 -1 0    0.25
VECTOR dataz SIZE 7 DATA 0 1 2  3    2 1.5  1

# using shepard's interpolation method
FUNCTION s(x,y) VECTORS datax datay dataz INTERPOLATION shepard SHEPARD_EXPONENT 4

# or using shepard's modified algorithm
FUNCTION m(x,y) VECTORS datax datay dataz INTERPOLATION modified_shepard SHEPARD_RADIUS 2

# print the definition points
PRINT_FUNCTION n FILE_PATH n_def.dat

# print the different functions at the selected range
PRINT_FUNCTION n s m FILE_PATH n_int.dat MIN -1 -1 MAX 1.5 1.5 STEP 0.05 0.05
$ wasora scattered.was
$ gnuplot scattered.gp
$ 
scattered1a.png
scattered1a.png
scattered2a.png
scattered2a.png
scattered3a.png
scattered3a.png
scattered1b.png
scattered1b.png
scattered2b.png
scattered2b.png
scattered3b.png
scattered3b.png
scattered1c.png
scattered1c.png
scattered2c.png
scattered2c.png
scattered3c.png
scattered3c.png
scattered2d1.png
scattered2d1.png
scattered2d2.png
scattered2d2.png
scattered2d3.png
scattered2d3.png

The figures illustrate how the multidimensional data interpolation scheme work for pointwise defined functions over scattered data. Nearest neighbors give constant values for each voronoi triangle whilst Shepard-based algorithms provide continuous surfaces.

compwater.was

This example shows an extension of the example about saturated water in section ref{007-functions} by giving properties of compressed water as a function of pressure \(p\) and temperature \(T\). The file compwater.txt contains some properties of water as a function of temperature and pressure in a rectangular grid over the pressure-temperature space. The enthalpy is not contained in the file, but it can be computed from pressure \(p\), the internal energy \(u(p,T)\) and the specific volume \(v(p,T)\) as

!bt [ h(p,t) = u(p,T) + p v(p,T) ] !et

FUNCTION v(p,T) FILE_PATH compwater.txt COLUMNS 2 1 4 INTERPOLATION rectangular
FUNCTION u(p,T) FILE_PATH compwater.txt COLUMNS 2 1 5 INTERPOLATION rectangular
FUNCTION h(p,T) = u(p,T) + p*v(p,T)

PRINT_FUNCTION h MIN 1e5 300 MAX 200e5 1000 STEP 2e5 4
$ wasora compwater.was > compwater.dat
$ gnuplot compwater.gp
$ 
compwater.png
compwater.png

The figure shows the enthalpy of compressed water as a continuous function of \(p\) and \(T\) and the discrete experimental data.