The ArrayInterpolation command takes a finite set of distinct data points given by xdata and ydata (or xydata), and interpolates to approximate the y-values corresponding to the points given in xvalues.  It considers an interpolant function $f$ such that $f\left(x\right)=y$ for all respective pairs $\left(x\,y\right)$ in xdata and ydata (or xydata). Such a function can be constructed using one of various methods (see below).  It then computes and returns $f\left(x_i\right)$ for all $x_i$ in xvalues.

The focus of the ArrayInterpolation command is the performance of quick and efficient data resampling and table lookup.  To actually compute and return interpolants, functions such as CurveFitting[Spline] and CurveFitting[RationalInterpolation] can be used instead.

The ArrayInterpolation function can interpolate numeric data in $n$ dimensions, where $n$ is any positive integer.

The list of independent coordinates of the data points, given in xdata, can be input in a number of different ways.  xdata can be:

(preferred if $n=1$) a Vector, DataSeries, list, or one-dimensional Array of strictly increasing x-coordinates. The data set will then have size $a_1$, where $a_1$ is the length of xdata.

(preferred if $1<n$) a list of $n$ Vectors, lists, or one-dimensional Arrays, one for each dimension of the data.  The $j$th Vector, list, or Array in the input must contain, in increasing order, all of the possible $j$th coordinates of the data points.  In this case, the block of data points will be assumed to lie on an $a_1$ by $a_2$ by ... by $a_n$ grid, where $a_j$ is the length of the $j$th Vector or Array in the input.  The $p$th coordinate of the data point at index $\&lsqb;j_1,j_2,\; ...,j_n\&rsqb;$ (where $1\le j_i\le a_i$) will be equal to the $a_p$th element of the $p$th Array in the input.

a list of $n$ Arrays of size $a_1$ by $a_2$ by ... by $a_n$, where the $j$th array contains the $j$th independent coordinate of each of the $a_1$ by $a_2$ by ... by $a_n$ data points.  The coordinates must form a proper "grid" of values, and must be sorted in strictly increasing order along each dimension.  More formally, $\mathrm{op}\left(p\&comma;\mathrm{xdata}\right)\&lsqb;j_1,j_2,\; ...,j_n]\; -\mathrm{op}\left(p\&comma;\mathrm{xdata}\right)\&lsqb;k_1,k_2,\; ...,k_n\&rsqb;$ must be zero if $j_p=k_p$, and must be positive if $k_p<j_p$.

The preferred methods minimize memory usage and execution time by avoiding unnecessary storage and verification of redundant data.  In all cases, xdata must contain real values of type numeric.

The list of dependent coordinates of the data points, given in ydata, must be input as an Array (or a Matrix, Vector, or list for appropriate values of $n$) of size $a_1$ by $a_2$ by ... by $a_n$, so that the value of $\mathrm{ydata}\&lsqb;j_1,j_2,\; ...,j_n\&rsqb;$ corresponds to the element in xdata of index $\&lsqb;j_1,j_2,\; ...,j_n\&rsqb;$. Values in ydata may be complex, but must be of type complex(extended_numeric). The output will not contain complex values unless ydata does.

As an alternate form of input, a single structure xydata containing all coordinates of the data points can be entered.  It can be formatted in one of the following ways:

an Array, DataFrame, or Matrix of size $a_1$ by $a_2$ by ... by $a_n$ by ($n+1$), giving the independent and dependent coordinate(s) of each of $a_1$ by $a_2$ by ... by $a_n$ data points as an ordered ($n+1$)-tuple.  The first n elements in each ($n+1$)-tuple represent the independent coordinates of each point, and must adhere to the same restrictions as above (a proper "grid" must be formed, and the independent coordinates must be sorted in strictly increasing order along each dimension).  The $n+1$st coordinate in each ($n+1$)-tuple then represents the dependent coordinate of the respective data point.

a list of $n+1$ Arrays, Vectors, Matrices, or lists of size $a_1$ by $a_2$ by ... by $a_n$, where the $j$th array contains the $j$th independent coordinate of each of the $a_1$ by $a_2$ by ... by $a_n$ data points for $1\le j\le n$, and the $n+1$st Array contains the dependent coordinates of each point.  As above, the independent coordinates must adhere to certain restrictions (a proper "grid" must be formed, and the independent coordinates must be sorted in strictly increasing order along each dimension).

For multidimensional data, these methods are not recommended, since space is wasted storing the full grid of independent coordinates instead of a list of all the possible coordinates in each dimension. In both cases, the values representing independent coordinates must be of type numeric, but the values representing dependent coordinates can be of type complex(extended_numeric).

The list of values to interpolate at, given in xvalues, may be input in one of the following formats:

for one-dimensional data, a single numeric value, or a Vector, list, or one-dimensional Array of numeric values can be input.  The output will be returned in a format matching the format of the input.

for multidimensional data, an Array or Matrix of size $u_1$ by $u_2$ by ... by $u_k$ by $n$ of numeric values can be input.  It must contain the $n$ coordinates of each of $u_1$ by $u_2$ by ... by $u_k$ values to interpolate at, with the value of $\mathrm{xvalues}\&lsqb;j_1,j_2,\; ...,j_k,p\&rsqb;$ giving the $p$th coordinate of the respective point.  The output will be returned in an array of size $u_1$ by $u_2$ by ... by $u_k$ containing the interpolated results.

alternatively, a list of $n$ Vectors, lists, or one-dimensional Arrays can be input. The $j$th Vector, list, or Array in the input will be assumed to contain all of the possible $j$th coordinates of the values to interpolate at.  In this case, interpolation will be performed on an $a_1$ by $a_2$ by ... by $a_n$ block of points, where $a_j$ is the length of the $j$th Vector or Array in the input.  The output will then be returned in a Vector, Matrix, list, or Array of size $a_1$ by $a_2$ by ... by $a_n$.

If any of the data points in xvalues lie outside the rectangular bounding box specified by the input, then extrapolation will be performed to approximate their corresponding y-values.  The method by which extrapolation is performed can be controlled by using option extrapolate; see below.

This routine has separate numeric methods for handling hardware and software floats.  The decision about which routine to use can be controlled by setting the UseHardwareFloats environment variable.  If UseHardwareFloats remains unset, then hardware floats are used if and only if $\mathrm{Digits}\le \mathrm{evalhf}\left(\mathrm{Digits}\right)$, in which case all software floats in the input will be converted to hardware floats.

Only computations involving numeric floating-point data are supported by this routine. If the input does not contain floating-point data, an error will be thrown.

For optimal performance, all rtables in the input should be Fortran order with rectangular storage (the default). Otherwise, a conversion will take place.  All rtables in the output will be Fortran order rtables with rectangular storage.

This function is part of the CurveFitting package, so it can be used in the short form ArrayInterpolation(..) only after executing the command with(CurveFitting).  However, it can always be accessed through the long form of the command by using CurveFitting[ArrayInterpolation](..).

If the option method = <name> is given, then one of the following interpolation methods are used to compute the interpolant $f$ and evaluate $f\left(x_i\right)$ for each point $x_i$ in xvalues:

method = nearest: Perform nearest neighbor interpolation.  Given a point $x_i$ in xvalues, $f\left(x_i\right)$ is defined to be $y$, where $\left(x\&comma;y\right)$ is the data point such that the Euclidean distance $‖x-x_i‖$ is minimized.

method = lowest: Perform lowest neighbor interpolation.  Given a point $x_i$ in xvalues, $f\left(x_i\right)$ is defined to be $y$, where$\left(x\&comma;y\right)$ is the data point such that $x_i-x$ is non-negative in all coordinates, but the Euclidean distance $‖x-x_i‖$ is minimized.

method = highest: Perform highest neighbor interpolation.  Given a point $x_i$ in xvalues, $f\left(x_i\right)$ is defined to be $y$, where $\left(x\&comma;y\right)$ is the data point such that $x-x_i$ is non-negative in all coordinates, but the Euclidean distance $‖x-x_i‖$ is minimized.

method = linear: Perform n-dimensional linear interpolation (lerping).  In the one-dimensional case, $f$ is a piecewise-linear function passing through each data point $\left(x\&comma;y\right)$ in the input.  In the multidimensional case, $f$ is the tensor product of $n$ such piecewise linear functions, one for each dimension.  $f\left(x_i\right)$ is computed by performing linear interpolation along the first dimension, then along the second dimension, and so on.

method = cubic: Perform piecewise cubic Hermite interpolation. In the 1-dimensional case, $f$ is a piecewise-cubic function passing through each data point$\left(x\&comma;y\right)$ in the input. In this case, $f\left(x\right)$ = $f_i\left(x\right)$ if $x$ lies in the interval $\left[x_i\&comma;x_{i+1}\right]$, where each $f_i$ is a cubic polynomial such that $f_i\left(x_i\right)=y_i$ and $f_i\left(x_{i+1}\right)=y_{i+1}$ for all data points $\left(x_i\&comma;y_i\right)$ in the input (where $i$ ranges from $0$ to $k$). The coefficients of the functions $f_i$ are determined locally by assigning slopes $s_i$ to each data point $x_i$ and solving for the unique cubic function $f_i\left(x\right)$ determined by the additional constraints that $f_i\text{'}\left(x_i\right)=s_i$and$f_i\text{'}\left(x_{i+1}\right)=s_{i+1}$.  This forces $f$ to be continuously differentiable ($\mathrm{C1}$).  The $s_i$ themselves are computed using Bessel's method: $s_i$ is the slope at $x_i$ of the parabola passing through $x_{i-1},y_{i-1}$, $x_i,y_i$, and $x_{i+1},y_{i+1}$. In the multidimensional case, $f$ is the tensor product of $n$ such spline functions, one for each dimension.

method = spline: Perform spline interpolation.  By default, natural cubic spline interpolation is used. In the 1-dimensional case, $f$ is a piecewise-cubic function passing through each data point$\left(x\&comma;y\right)$ in the input. In this case, $f\left(x\right)$ = $f_i\left(x\right)$ if $x$ lies in the interval $\left[x_i\&comma;x_{i+1}\right]$, where each $f_i$ is a cubic polynomial such that $f_i\left(x_i\right)=y_i$ and $f_i\left(x_{i+1}\right)=y_{i+1}$ for all data points $\left(x_i\&comma;y_i\right)$ in the input (where $i$ ranges from $0$ to $k$).  The coefficients of the functions $f_i$ are selected such that $f$ is twice continuously differentiable ($\mathrm{C2}$), that is,  $f_i\text{'}\left(x_{i+1}\right)=f_{i+1}\text{'}\left(x_{i+1}\right)$and $f_i\text{'}\text{'}\left(x_{i+1}\right)=f_{i+1}\text{'}\text{'}\left(x_{i+1}\right)$.  In addition, the "natural" condition of the spline specifies that $f\text{'}\text{'}\left(x_0\right)=0$ and $f\text{'}\text{'}\left(x_{\mathrm{k}}\right)=0$. In the multidimensional case, $f$ is the tensor product of $n$ such spline functions, one for each dimension. Using method=spline will produce a smoother interpolant than method=cubic ($\mathrm{C2}$ instead of $\mathrm{C1}$), but is more expensive to set up and more prone to numerical instability because each segment of the spline is determined globally by the positions of all other points in the data set.

method=linear is used by default.

If the options degree=d and endpoints=e are given, where d is a positive integer and e is one of natural, notaknot, or periodic, then spline interpolation will be performed using the provided degree and endpoint conditions.  See Spline Continuity and End Conditions for details. These options only affect the result if method=spline is used.  In the multidimensional case, the same degree and endpoint conditions are used for the splines generated in each dimension. The defaults are degree=3 and endpoints=natural, in which case natural cubic spline interpolation will be performed.

If splines of an even degree are being used, the option knots=data forces the use of a spline function where the spline knots are positioned on the nodes.  See Spline Continuity and End Conditions for details.  The default method, knots=default, defines the spline knots at the midpoints of the nodes when even degree splines are used.  This option has no effect when other methods are used.

If the option uniform=true is given, then ArrayInterpolation assumes that the data points are sampled over a grid of uniformly spaced points in each dimension. In other words, if $a_{i,j}$ is the $j$th possible coordinate in the $i$th dimension, then $d_i=a_{i,j}-a_{i,j-1}$ is assumed to be constant over all possible $j$, given any fixed value of $i$. This gives a considerable speedup when the input contains uniform data, since it allows ArrayInterpolation to use a fast lookup algorithm when evaluating the interpolant at the specified points. The default is uniform=false, in which case ArrayInterpolation uses a slower but more general binary search algorithm to perform interpolation.  Using the uniform=true option with non-uniform data may produce incorrect results.

If the option verify=false is given, then ArrayInterpolation skips the various checks it performs to ensure correctly formatted input.  This can decrease the time required to solve large problems, but will prevent the function from detecting any errors in the input.  If the input is improperly sorted, contains Arrays indexed from values other than 1, contains non-rectangular or C order rtables, or is otherwise formatted incorrectly, ArrayInterpolation may return incorrect results or throw an unexpected error.

If the option extrapolate=e is given, where e is of type extended_numeric or truefalse, then one of the following possible extrapolation methods will be used to compute $f\left(x\right)$ if $x$ lies outside the bounding box specified by the input:

extrapolate = true: Perform extrapolation using the closest valid branch of the interpolating function. In the case of method=lowest and method=highest, this is not be defined for some points, in which case undefined will be returned.

extrapolate = false: Do not extrapolate.  An error will be thrown if any point in xvalues lies outside the bounding box specified by the input.

extrapolate = e, where e is of type extended_numeric : Define $f\left(x\right)$ to be e if $x$ lies outside the bounding box specified by the input.  e is commonly zero or undefined.

extrapolate=true is used by default.

If the option container=c is given, where c is an appropriately sized and typed rtable, then the result is written into c. c must be of the correct size and datatype to match the output of the routine. If ydata (and thus c) have only real values, then with this option, no additional memory is allocated to store the result; this is a programmer-level feature that can be used to reduce memory usage and decrease the time spent by Maple's garbage collector. If nonreal values are present, then the output is still copied into c, but there are no memory savings.  The default is container=false, in which case Maple creates and returns a new rtable containing the result.

$\mathrm{with}\left(\mathrm{CurveFitting}\right)\&colon;$

$\mathrm{with}\left(\mathrm{plots}\right)\&colon;$

$\mathrm{Times}≔\left[0.\&comma;0.01\&comma;0.02\&comma;0.03\&comma;0.04\&comma;0.05\&comma;0.06\&comma;0.07\&comma;0.08\&comma;0.09\&comma;0.10\right]\&colon;$

$\mathrm{Amplitudes}≔\left[-0.6\&comma;0.\&comma;0.4\&comma;0.6\&comma;0.3\&comma;-0.1\&comma;-0.2\&comma;0.\&comma;0.1\&comma;-0.3\&comma;-0.6\right]\&colon;$

$\mathrm{pointplot}\left(\mathrm{Times}\&comma;\mathrm{Amplitudes}\right)$

$\mathrm{NewTimes}≔\left[\mathrm{seq}\left(0.001i\&comma;i=0..100\right)\right]\&colon;$

$\mathrm{NewAmplitudes}≔\mathrm{ArrayInterpolation}\left(\mathrm{Times}\&comma;\mathrm{Amplitudes}\&comma;\mathrm{NewTimes}\right)\&colon;$

$\mathrm{pointplot}\left(\mathrm{NewTimes}\&comma;\mathrm{NewAmplitudes}\right)$

$\mathrm{NewAmplitudes}≔\mathrm{ArrayInterpolation}\left(\mathrm{Times}\&comma;\mathrm{Amplitudes}\&comma;\mathrm{NewTimes}\&comma;\mathrm{method}=\mathrm{spline}\right)\&colon;$

$\mathrm{pointplot}\left(\mathrm{NewTimes}\&comma;\mathrm{NewAmplitudes}\right)$

$\mathrm{NewAmplitudes}≔\mathrm{ArrayInterpolation}\left(\mathrm{Times}\&comma;\mathrm{Amplitudes}\&comma;\mathrm{NewTimes}\&comma;\mathrm{method}=\mathrm{spline}\&comma;\mathrm{endpoints}=\mathrm{periodic}\right)\&colon;$

$\mathrm{pointplot}\left(\mathrm{NewTimes}\&comma;\mathrm{NewAmplitudes}\right)$

$\mathrm{Ranges}≔\left[\left[\mathrm{seq}\left(1..10\right)\right]\&comma;\left[\mathrm{seq}\left(1..16\right)\right]\right]\&colon;$

$\mathrm{Img}≔\mathrm{evalf}\left(\mathrm{Matrix}\left(\left[\left[255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\right]\&comma;\left[255\&comma;255\&comma;208\&comma;93\&comma;22\&comma;0\&comma;16\&comma;98\&comma;210\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\right]\&comma;\left[255\&comma;196\&comma;16\&comma;0\&comma;0\&comma;0\&comma;0\&comma;0\&comma;19\&comma;200\&comma;255\&comma;255\&comma;54\&comma;54\&comma;212\&comma;255\right]\&comma;\left[255\&comma;61\&comma;34\&comma;156\&comma;231\&comma;255\&comma;230\&comma;154\&comma;26\&comma;58\&comma;255\&comma;255\&comma;61\&comma;0\&comma;59\&comma;255\right]\&comma;\left[255\&comma;5\&comma;205\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;196\&comma;7\&comma;255\&comma;255\&comma;230\&comma;91\&comma;7\&comma;255\right]\&comma;\left[255\&comma;20\&comma;235\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;235\&comma;36\&comma;255\&comma;255\&comma;255\&comma;233\&comma;32\&comma;255\right]\&comma;\left[255\&comma;125\&comma;129\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;133\&comma;186\&comma;255\&comma;255\&comma;217\&comma;89\&comma;147\&comma;255\right]\&comma;\left[255\&comma;255\&comma;125\&comma;50\&comma;91\&comma;99\&comma;93\&comma;40\&comma;22\&comma;121\&comma;93\&comma;58\&comma;10\&comma;140\&comma;255\&comma;255\right]\&comma;\left[255\&comma;255\&comma;255\&comma;216\&comma;131\&comma;75\&comma;14\&comma;5\&comma;7\&comma;12\&comma;59\&comma;121\&comma;214\&comma;255\&comma;255\&comma;255\right]\&comma;\left[255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\&comma;255\right]\right]\right)\right)\&colon;$

$\mathrm{listdensityplot}\left(\mathrm{Img}\&comma;\mathrm{style}=\mathrm{PATCHNOGRID}\&comma;\mathrm{axes}=\mathrm{none}\&comma;\mathrm{range}=0..255\right)$

$\mathrm{NewRanges}≔\frac{\left[\left[\mathrm{seq}\left(3..52\right)\right]\&comma;\left[\mathrm{seq}\left(3..82\right)\right]\right]}{5.0}\&colon;$

$\mathrm{NewImg}≔\mathrm{ArrayInterpolation}\left(\mathrm{Ranges}\&comma;\mathrm{Img}\&comma;\mathrm{NewRanges}\&comma;\mathrm{extrapolate}=255\right)\&colon;$

$\mathrm{listdensityplot}\left(\mathrm{NewImg}\&comma;\mathrm{style}=\mathrm{PATCHNOGRID}\&comma;\mathrm{axes}=\mathrm{none}\&comma;\mathrm{range}=0..255\right)$

$\mathrm{NewImg}≔\mathrm{ArrayInterpolation}\left(\mathrm{Ranges}\&comma;\mathrm{Img}\&comma;\mathrm{NewRanges}\&comma;\mathrm{extrapolate}=255\&comma;\mathrm{method}=\mathrm{cubic}\right)\&colon;$

$\mathrm{listdensityplot}\left(\mathrm{NewImg}\&comma;\mathrm{style}=\mathrm{PATCHNOGRID}\&comma;\mathrm{axes}=\mathrm{none}\&comma;\mathrm{range}=0..255\right)$

$f≔\left(i\&comma;j\right)\mapsto {\left(3-\sin \left(i\right)\right)}^2-{\left(3-j\right)}^2\&colon;$

$\mathrm{plot3d}\left(f\&comma;0..10\&comma;0..10\&comma;\mathrm{axes}=\mathrm{normal}\right)$

$v≔\mathrm{Array}\left(\left[0\&comma;1.5\&comma;3.5\&comma;5\&comma;6\&comma;8\&comma;9\right]\right)\&colon;$

$w≔\mathrm{Array}\left(\left[0\&comma;3\&comma;5\&comma;6\&comma;6.5\&comma;7.5\right]\right)\&colon;$

$y≔\mathrm{Matrix}\left(7\&comma;6\&comma;\left(a\&comma;b\right)\mapsto \mathrm{evalf}\left(f\left(v\left[a\right]\&comma;w\left[b\right]\right)\right)\right)\&colon;$

$\mathrm{pointplot3d}\left(\left[\mathrm{seq}\left(\mathrm{seq}\left(\left[v\left[i\right]\&comma;w\left[j\right]\&comma;y\left[i,j\right]\right]\&comma;j=1..6\right)\&comma;i=1..7\right)\right]\&comma;\mathrm{axes}=\mathrm{normal}\&comma;\mathrm{symbol}=\mathrm{sphere}\right)$

$\mathrm{a1}≔\mathrm{Matrix}\left(50\&comma;50\&comma;\left(i\&comma;j\right)\mapsto \frac{i}{5}\right)\&colon;$

$\mathrm{a2}≔\mathrm{Matrix}\left(50\&comma;50\&comma;\left(i\&comma;j\right)\mapsto \frac{j}{5}\right)\&colon;$

$A≔\mathrm{ArrayTools}\left[\mathrm{Concatenate}\right]\left(3\&comma;\mathrm{a1}\&comma;\mathrm{a2}\right)\&colon;$

$B≔\mathrm{ArrayInterpolation}\left(\left[v\&comma;w\right]\&comma;y\&comma;A\&comma;\mathrm{method}=\mathrm{linear}\right)\&colon;$

$\mathrm{matrixplot}\left(\mathrm{Matrix}\left(B\right)\&comma;\mathrm{axes}=\mathrm{normal}\right)$

$B≔\mathrm{ArrayInterpolation}\left(\left[v\&comma;w\right]\&comma;y\&comma;A\&comma;\mathrm{method}=\mathrm{nearest}\right)\&colon;$

$\mathrm{matrixplot}\left(\mathrm{Matrix}\left(B\right)\&comma;\mathrm{axes}=\mathrm{normal}\right)$

$B≔\mathrm{ArrayInterpolation}\left(\left[v\&comma;w\right]\&comma;y\&comma;A\&comma;\mathrm{method}=\mathrm{spline}\right)\&colon;$

$\mathrm{matrixplot}\left(\mathrm{Matrix}\left(B\right)\&comma;\mathrm{axes}=\mathrm{normal}\right)$

$B≔\mathrm{ArrayInterpolation}\left(\left[v\&comma;w\right]\&comma;y\&comma;A\&comma;\mathrm{method}=\mathrm{spline}\&comma;\mathrm{degree}=5\right)\&colon;$

$\mathrm{matrixplot}\left(\mathrm{Matrix}\left(B\right)\&comma;\mathrm{axes}=\mathrm{normal}\right)$

$f≔x\mapsto \sin \left(x+\mathrm{I}\right)$

$f≔x\mapsto \sin \left(x+\mathrm{I}\right)$

$\mathrm{vx}≔\left[\mathrm{seq}\left(-3..3\&comma;0.5\right)\right]\&colon;$

$\mathrm{vy}≔\left[\mathrm{seq}\left(\mathrm{evalf}\left(f\left(i\right)\right)\&comma;i=-3..3\&comma;0.5\right)\right]\&colon;$

$\mathrm{complexplot}\left(\mathrm{vy}\&comma;\mathrm{style}=\mathrm{point}\&comma;\mathrm{scaling}=\mathrm{constrained}\right)$

$A≔\left[\mathrm{seq}\left(i\&comma;i=-\mathrm{evalf}\left(\mathrm{\pi }\right)..\mathrm{evalf}\left(\mathrm{\pi }\right)\&comma;\mathrm{numelems}=100\right)\right]\&colon;$

$B≔\mathrm{ArrayInterpolation}\left(\mathrm{vx}\&comma;\mathrm{vy}\&comma;A\&comma;\mathrm{method}=\mathrm{spline}\right)\&colon;$

$\mathrm{complexplot}\left(\mathrm{convert}\left(B\&comma;\mathrm{list}\right)\&comma;\mathrm{scaling}=\mathrm{constrained}\right)$

$f≔\left(i\&comma;j\right)\mapsto \mathrm{sqrt}\left(i\right)+\mathrm{sqrt}\left(j\right)\&colon;$

$v≔\mathrm{Array}\left(\left[0\&comma;2\&comma;4\&comma;5\&comma;8\&comma;10\&comma;12\right]\right)\&colon;$

$w≔\mathrm{Array}\left(\left[-6\&comma;-5\&comma;-4\&comma;-3\&comma;-2\&comma;-1\right]\right)\&colon;$

$y≔\mathrm{Matrix}\left(7\&comma;6\&comma;\left(a\&comma;b\right)\mapsto \mathrm{evalf}\left(f\left(v\left[a\right]\&comma;w\left[b\right]\right)\right)\right)\&colon;$

$\mathrm{a1}≔\mathrm{Matrix}\left(50\&comma;50\&comma;\left(i\&comma;j\right)\mapsto \frac{i}{5}\right)\&colon;$

$\mathrm{a2}≔\mathrm{Matrix}\left(50\&comma;50\&comma;\left(i\&comma;j\right)\mapsto \frac{j}{5}\right)\&colon;$

$A≔\mathrm{ArrayTools}\left[\mathrm{Concatenate}\right]\left(3\&comma;\mathrm{a1}\&comma;\mathrm{a2}\right)\&colon;$

$B≔\mathrm{ArrayInterpolation}\left(\left[v\&comma;w\right]\&comma;y\&comma;A\&comma;\mathrm{method}=\mathrm{spline}\right)\&colon;$

$\mathrm{complexplot}\left(\mathrm{convert}\left(B\&comma;\mathrm{list}\right)\&comma;\mathrm{axes}=\mathrm{normal}\right)$

$\mathrm{xdata}≔\left[\mathrm{Matrix}\left(7\&comma;6\&comma;\left(i\&comma;j\right)\mapsto v\left[i\right]\right)\&comma;\mathrm{Matrix}\left(7\&comma;6\&comma;\left(i\&comma;j\right)\mapsto w\left[j\right]\right)\right]\&colon;$

$B≔\mathrm{ArrayInterpolation}\left(\mathrm{xdata}\&comma;y\&comma;A\&comma;\mathrm{method}=\mathrm{linear}\right)\&colon;$

$\mathrm{v1}≔\mathrm{Array}\left(\left[-10\&comma;-9\&comma;-8\&comma;-7\right]\right)\&colon;$

$\mathrm{v2}≔\mathrm{Array}\left(\left[-1\&comma;0\&comma;1\&comma;2\&comma;3\right]\right)\&colon;$

$\mathrm{v3}≔\mathrm{Array}\left(\left[-3\&comma;-2\&comma;0\&comma;6\right]\right)\&colon;$

$\mathrm{v4}≔\mathrm{Array}\left(\left[2\&comma;3\&comma;4\&comma;5\&comma;6\&comma;7\right]\right)\&colon;$

$f≔\left(i\&comma;j\&comma;h\&comma;k\right)\mapsto k\cdot \mathrm{sqrt}\left(i\right)+\mathrm{sqrt}\left(h\right)-j\&colon;$

$y≔\mathrm{Array}\left(1..4\&comma;1..5\&comma;1..4\&comma;1..6\&comma;\left(i\&comma;j\&comma;h\&comma;k\right)\mapsto \mathrm{evalf}\left(f\left(\mathrm{v1}\left[i\right]\&comma;\mathrm{v2}\left[j\right]\&comma;\mathrm{v3}\left[h\right]\&comma;\mathrm{v4}\left[k\right]\right)\right)\right)\&colon;$

$A≔\left[\mathrm{Array}\left(\left[\mathrm{seq}\left(i\&comma;i=-10..-7\&comma;0.2\right)\right]\right)\&comma;\mathrm{Array}\left(\left[\mathrm{seq}\left(i\&comma;i=-1..3\&comma;0.2\right)\right]\right)\&comma;\mathrm{Array}\left(\left[\mathrm{seq}\left(i\&comma;i=-3..6\&comma;1.5\right)\right]\right)\&comma;\mathrm{Array}\left(\left[\mathrm{seq}\left(i\&comma;i=2..7\&comma;0.4\right)\right]\right)\right]\&colon;$

$B≔\mathrm{ArrayInterpolation}\left(\left[\mathrm{v1}\&comma;\mathrm{v2}\&comma;\mathrm{v3}\&comma;\mathrm{v4}\right]\&comma;y\&comma;A\right)\&colon;$

$A≔\mathrm{evalf}\left(\mathrm{Vector}\left(\left[\mathrm{seq}\left(0..50000\right)\right]\right)\right)\&colon;$

$B≔\mathrm{evalf}\left(\mathrm{Vector}\left(\left[\mathrm{seq}\left(i^2\&comma;i=0..50000\right)\right]\right)\right)\&colon;$

$C≔\mathrm{evalf}\left(\mathrm{Vector}\left(\frac{\left[\mathrm{seq}\left(0..200000\right)\right]}{4}\right)\right)\&colon;$

$\mathrm{time}\left(\mathrm{ArrayInterpolation}\left(A\&comma;B\&comma;C\right)\right)$

$0.066$

$\mathrm{time}\left(\mathrm{ArrayInterpolation}\left(A\&comma;B\&comma;C\&comma;\mathrm{verify}=\mathrm{false}\right)\right)$

$0.018$

$\mathrm{time}\left(\mathrm{ArrayInterpolation}\left(A\&comma;B\&comma;C\&comma;\mathrm{verify}=\mathrm{false}\&comma;\mathrm{uniform}=\mathrm{true}\right)\right)$

$0.012$

$\mathrm{time}\left(\mathrm{ArrayInterpolation}\left(A\&comma;B\&comma;C\&comma;\mathrm{verify}=\mathrm{false}\&comma;\mathrm{uniform}=\mathrm{true}\&comma;\mathrm{method}=\mathrm{cubic}\right)\right)$

$0.015$