10.4 Precision
Let us consider the butterfly curve, proposed by the mathematician Temple Fay. This periodic curve is defined by the following equation in polar coordinates:
\[r=e^{\sin \phi} - 2 \cos (4 \phi ) + \sin^5\left(\frac{2 \phi - \pi}{24}\right)\]
Similar to what we did before, we will consider the curve with respect to a starting point \(P\), within a range of variation of the independent parameter \(t \in [t_0, t_1]\).
butterfly_curve(p, t0, t1, n) =
map_division(t -> p+vpol(exp(sin(t))-2*cos(4*t)+sin((2*t-pi)/24)^5, t),
t0,
t1,
n,
false)
To draw this curve, we need to determine the appropriate value of the parameter \(n\). For this, it is important to realize that the period of the butterfly curve is \(24 \pi\), being the entire curve generated for \(t \in [0, 24\pi]\). Since this period is relatively large, the number \(n\) of points to be used should be enough for the curve to be accurately reproduced. However, this raises a difficulty: without knowing the shape of the curve, it is difficult to estimate which is the best value for \(n\). The lower the value of \(n\), the less accurate will be the curve representation, in particular for very high curvatures; the higher the value of \(n\), the greater the need for computational resources. In practice, some experimentation is required to understand what best value is, as evidenced in this figure that shows the butterfly curves drawn for the following increasing values of \(n\):
closed_spline(butterfly_curve(xy(0, 10), 0, 24*pi, 10))
closed_spline(butterfly_curve(xy(10, 10), 0, 24*pi, 20))
closed_spline(butterfly_curve(xy(20, 10), 0, 24*pi, 40))
closed_spline(butterfly_curve(xy(30, 10), 0, 24*pi, 80))
closed_spline(butterfly_curve(xy(40, 10), 0, 24*pi, 160))
closed_spline(butterfly_curve(xy(0, 0), 0, 24*pi, 320))
closed_spline(butterfly_curve(xy(10, 0), 0, 24*pi, 640))
closed_spline(butterfly_curve(xy(20, 0), 0, 24*pi, 1280))
closed_spline(butterfly_curve(xy(30, 0), 0, 24*pi, 2560))
closed_spline(butterfly_curve(xy(40, 0), 0, 24*pi, 5120))
The butterfly curve produced for increasing values of the number \(n\) of points.
As observable, when the number of points is insufficient, the produced curves can be grotesque distortions of reality. For this reason, we should take particular care in choosing the most suitable value for this parameter.
An interesting alternative will be to give the computer the responsibility of adapting the number of points to the necessities of each curve, as it will be seen in the next section.
10.4.1 Adaptive Sampling
The map_division function allows us to generate the curve of any function \(f(t)\) in the interval \([t_0, t_1]\) by applying the function \(f\) to a sequence of \(n\) increments of \(t\) equally spaced \(\{t_0,\ldots,t_1\}\). This sequence is named sampling sequence.
As we have seen, the greater the number of sampling elements, the greater the accuracy of the produced curve. However, the greater the accuracy, the greater the effort expended in computing the function value for all the elements of the sequence. In practice, it is necessary to find a balance between accuracy and computational effort.
Unfortunately, there are several functions for which it is difficult or even impossible to find an appropriate number of sampling elements. As an example, let us try to produce the graph of the function \(f(t)=(t,\frac{1}{2}sin(\frac{1}{t}))\) in the interval \([0.05,0.8]\). The graph of this function is a curve that oscillates with a frequency that is greater the closer it is to the origin. This figure shows a sequence of graphs of this function in which, from left to right, the number of sampling points was progressively increased.
Graph of the function \(f(x)=\frac{1}{2}sin(\frac{1}{x})\) in the interval \([0.05,0.8]\). From left to right, \(8\), \(15\), \(30\) and \(60\) sampling points were employed.
As we can see in the figure, with few sampling points the resulting curve is a rough approximation of the real curve; it is only when the number of sampling points becomes very high that some precision is acquired. However, precision is only relevant in areas of the curve where large variations occur; in the other areas, where the curve evolves smoothly, more sampling points only imply greater processing time, without any significant improvement in the quality of the result. This suggests that we should adapt the number of sampling points along the curve according to its smoothness: in areas where the curvature is higher, we should employ more sampling points, and in areas where the curvature is lower, we should decrease the number of sampling points. This way, the number of sampling points adapts to shape of the curve.
To employ an adaptive sampling we need to define a criterion to classify the smoothness of a curve. As illustrated in this figure, we can start by using two abscissas \(x_0\) and \(x_1\) that we use to calculate the points \(P_0\) and \(P_1\). To know if the curve behaves in a linear fashion between those two points (and thus can be approximated to a line segment), we will calculate the average value \(x_m\) of that interval and the corresponding point \(P_m\) of the curve. If the function is in fact linear in that interval, then the points \(P_0\), \(P_m\) and \(P_1\) are collinear. In practice, that will rarely happen, so we will have to use a concept of approximate collinearity, for which we can use any of the several possible criteria, such as:
The area of the triangle \(P_0,P_m,P_1\) to be smaller than a sufficiently small \(\epsilon\) value;
The angle \(P_0,P_m,P_1\) to be approximately equal to \(\pi\);
The sum of the distance between \(P_0\) and \(P_m\) with the distance between \(P_m\) and \(P_1\) to be approximately equal to the distance between \(P_0\) and \(P_1\).
Adaptive generation of points.
For now, we will choose the first criterion. Its translation into Julia is as follows:
is_collinear(p0, pm, p1, epsilon) =
let a = distance(p0, pm),
b = distance(pm, p1),
c = distance(p1, p0)
triangle_area(a, b, c) < epsilon
end
triangle_area(a, b, c) =
let s = (a+b+c)/2.0
sqrt(s*(s-a)*(s-b)*(s-c))
end
The is_collinear function implements the criterion that allows us to say whether the line segments connecting the points p0, pm, and p1 are a good approximation to the curve between those points. When this criterion does not occur, we can split the interval into two subintervals and analyze each of them separately, concatenating the results.
map_division_adapt(f, t0, t1, epsilon) =
let p0 = f(t0),
p1 = f(t1),
tm = (t0+t1)/2.0,
pm = f(tm)
if is_collinear(p0, pm, p1, epsilon)
[p0, pm, p1]
else
[map_division_adapt(f, t0, tm, epsilon)...,
map_division_adapt(f, tm, t1, epsilon)[2:end]...]
end
end
This adaptive way of producing a sample of a curve’s points allows more accurate results, even when employing a relatively small number of sampling points. This figure shows a sequence of graphs similar to those shown in this figure, but now using an adaptive sampling approach with approximately the same number of sampling points in each case. As we can see, the graphs are substantially more accurate, particularly in areas close to the origin, without a noticeable decrease in the quality of the curve in areas where the change is smoother.
Graph of the function \(f(x)=\frac{1}{2}sin(\frac{1}{x})\) in the interval \([0.05,0.8]\). From left to right, the sampling sequence has, respectively, \(7\), \(15\), \(29\) and \(57\) points.
10.4.2 Exercises 51
10.4.2.1 Question 188
The algorithm used in the map_division_adapt function cannot properly deal with periodic functions in which we have \(f(t_0)=f(t_1)=f(\frac{t_0+t_1}{2})\). In this case, the degenerate triangle formed by these three points has zero area, which makes the algorithm terminate immediately. Modify the function so it receives an additional parameter \(\Delta_t\) that determines the biggest admissible interval between two consecutive values of \(t\).
10.4.2.2 Question 189
The map_division_adapt function is not as efficient as it could be, in particular, because it systematically repeats the application of the function \(f\) at the midpoint: \(f(t_m)\) is calculated to decide whether the recursion continues and, when it continues, it will have to be calculated again on the following function call. Redefine the map_division_adapt function in order to avoid repetitive calculations.