Numerical Computation Project

1. Parametric Equation

Suppose we are machining a curve on a CNC machine. This curve is typically defined by a parametric equation (e.g., a cubic Bézier curve).

1. Curve Position Definition:

Let the parameter be $u$ (typically $u\in [0,1]$). The position vector is:

$$\mathbf{r}(u)=\left[\begin{array}{c}x(u)\\ y(u)\end{array}\right]$$

2. Velocity Vector (Tangent Vector):

Differentiating with respect to $u$ gives the tangent vector:

$${\mathbf{r}}^{\prime }(u)=\left[\begin{array}{c}{x}^{\prime }(u)\\ y^{\prime }(u)\end{array}\right]$$

3. Instantaneous Speed (The Speed Function):

This is the most critical function. The magnitude of the velocity is:

$$v(u)=\parallel {\mathbf{r}}^{\prime }(u)\parallel =\sqrt{(x^{\prime }(u){)}^2+(y^{\prime }(u){)}^2}$$

Challenge: Even if $x(u)$ and $y(u)$ are simple polynomials, the term under the square root in $v(u)$ makes the function extremely complex. It is usually impossible to find an analytical antiderivative (closed-form solution). This is why Numerical Integration is required.

---

2. The Arc-Length Function

We need to calculate the actual distance traveled (arc length) $s$ when moving from the starting point $u=0$ to a specific parameter $t$. According to the Fundamental Theorem of Calculus:

$$s(t)={\int }_0^{t}v(u)du={\int }_0^t\sqrt{(x^{\prime }(u){)}^2+(y^{\prime }(u){)}^2}du$$

In this project, you need to implement a function calculate_arc_length(t).

Since this cannot be integrated manually, this function must internally call the Adaptive Simpson's Integration method to compute the definite integral.

---

3. Adaptive Simpson's Integration (Method & Derivation)

In numerical control (CNC), machining paths require high precision. A simple fixed-step integration (e.g., dividing the curve into 1000 equal segments) is inefficient: it wastes computation on straight lines and loses precision on sharp corners.

To solve this, we use Adaptive Simpson's Quadrature. This method automatically detects where the function changes rapidly (high curvature) and refines the grid locally using a "Divide and Conquer" strategy.

3.1 The Foundation: Simpson's 1/3 Rule

Simpson's rule improves upon the Trapezoidal rule (which fits straight lines) by fitting a quadratic polynomial (parabola) through three points.

For an arbitrary interval $[a,b]$:

1. Let the step size be half the interval width: $h=\frac{b-a}{2}$.

2. Let the midpoint be $c=\frac{a+b}{2}$.

The integral is approximated by the area under the parabola passing through $(a,v(a))$, $(c,v(c))$, and $(b,v(b))$:

$$S(a,b)=\frac{h}{3}[v(a)+4v(c)+v(b)]$$

Theoretical Error:

From Taylor series analysis, the error term $E$ for Simpson's rule depends on the fourth derivative of the function $v^{(4)}(\xi )$ and the step size to the 5th power:

$$E(h)=-\frac{1}{90}{h}^5{v}^{(4)}(\xi )$$

3.2 Derivation of the Error Estimator (Richardson Extrapolation)

In a practical CNC system, we do not know the fourth derivative $v^{(4)}$. Therefore, we cannot calculate the error directly. Instead, we use the technique of Richardson Extrapolation to estimate the current error by comparing two different approximations.

Assumption: For a sufficiently small interval $[a,b]$, the fourth derivative $v^{(4)}$ is approximately constant. Let $K=-\frac{1}{90}{v}^{(4)}(\xi )$.

Step 1: The Coarse Estimate ($S_1$)

We calculate the integral using one single Simpson step over the full interval $[a,b]$ (step size $h$).

Let $I$ be the true integral value. The coarse approximation relates to the true value as:

$$I=S_1+K\cdot h^5\dots \text{(Eq. 1)}$$

Step 2: The Fine Estimate ($S_2$)

We split the interval $[a,b]$ into two sub-intervals: $[a,c]$ and $[c,b]$. We apply Simpson's rule to each half separately and sum them up.

• The step size for each half is now $\frac{h}{2}$.
• The error for the left half is $K\cdot (\frac{h}{2}{)}^5$.
• The error for the right half is $K\cdot (\frac{h}{2}{)}^5$.

The total error for $S_2$ is the sum of these two errors:

$$E(S_2)\approx K{\left(\frac{h}{2}\right)}^5+K{\left(\frac{h}{2}\right)}^5=2\cdot K\frac{h^5}{32}=\frac{1}{16}K{h}^5$$

Thus, the fine approximation relates to the true value as:

$$I=S_2+\frac{1}{16}K{h}^5\dots \text{(Eq. 2)}$$

Key Insight: By halving the step size, the error of Simpson's rule drops to 1/16th of its original value.

Step 3: Extracting the Error

We now have a system of two linear equations with two unknowns ($I$ and $K{h}^5$). We want to solve for the error of the fine estimate, which is $E_2=\frac{1}{16}K{h}^5$.

Subtract (Eq. 1) from (Eq. 2):

$$I-I=(S_2+\frac{1}{16}K{h}^5)-(S_1+K{h}^5)$$$$0=(S_2-S_1)+K{h}^5(\frac{1}{16}-1)$$$$0=(S_2-S_1)-\frac{15}{16}K{h}^5$$

Rearranging to solve for the term $\frac{1}{16}K{h}^5$ (which is exactly $E_2$):

$$S_2-S_1=15\cdot \left(\frac{1}{16}K{h}^5\right)$$$$S_2-S_1=15\cdot E_2$$$$E_2=\frac{S_2-S_1}{15}$$

3.3 The Recursive Algorithm

We define a tolerance TOL (e.g., ${10}^{-6}$). The algorithm is recursive:

1. Compute Estimates: Calculate $S_1$ (one-step) and $S_2$ (two-step).

2. Evaluate Error: Calculate the estimated error magnitude: $|E|=\frac{|S_2-S_1|}{15}$.

3. Check Condition:

   • Success ($|E|<\text{TOL}$): The interval is smooth enough.

     • Bonus Accuracy: Instead of returning $S_2$, we return $S_2+E$. This simple addition (adding the estimated error back to the result) increases the order of accuracy from $O(h^5)$ to $O(h^7)$. This is known as Lyness's or Boole's Rule.

   • Failure ($|E|\ge \text{TOL}$): The function is varying too rapidly (high error).

     • Recursion: We simply cannot trust the result here. We split the interval in half and call the function recursively on the left side $[a,c]$ and the right side $[c,b]$.

     • Tolerance Halving: We pass TOL/2 to the sub-calls. This ensures that the sum of errors from the sub-intervals ($TOL/2+TOL/2$) does not exceed the total allowed error ($TOL$).

Fail-Safes:

In the actual code implementation, we include two critical checks to prevent infinite loops:

• Max Depth: Stop recursion if the depth exceeds a limit (e.g., 20 levels).

• Min Interval: Stop if the interval width $|b-a|$ becomes smaller than machine epsilon.

4. The Inverse Problem (Finding Parameter $t$)

Goal: We want the tool head to move at a constant velocity $V_{const}$.

Assume that at physical time $\tau$, the tool head should have traveled a target distance $L_{target}=V_{const}\times \tau$.

Problem: We know the target distance $L_{target}$, and we need to find the corresponding curve parameter $t^{\ast }$ such that:

$$s(t^{\ast })=L_{target}$$

This requires solving the integral equation:

$${\int }_0^{t^{\ast }}\sqrt{(x^{\prime }(u){)}^2+(y^{\prime }(u){)}^2}du=L_{target}$$

This is a Root-Finding Problem. We define a function $F(t)$:

$$F(t)=s(t)-L_{target}=({\int }_0^{t}v(u)du)-L_{target}=0$$

We need to find $t$ such that $F(t)=0$. The standard method for this is the Newton-Raphson Method.

Derivation of Newton's Iteration:

The iteration formula is $t_{k+1}=t_k-\frac{F(t_k)}{F^{\prime }(t_k)}$.

1. The Numerator $F(t_k)$:

   $$F(t_k)=s(t_k)-L_{target}$$

   Here, $s(t_k)$ must be computed by calling your Adaptive Simpson's Integration, passing the speed function $v(u)$ as the integrand.

2. The Denominator $F^{\prime }(t_k)$:

   According to the Fundamental Theorem of Calculus, the derivative of an accumulation function is simply the integrand itself!

   $$F^{\prime }(t)=\frac{d}{dt}({\int }_0^{t}v(u)du-L_{target})=v(t)$$

   Therefore, the denominator is simply the Instantaneous Speed $\sqrt{x^{\prime }(t{)}^2+y^{\prime }(t{)}^2}$.

5. Final Iterative Algorithm Formula

Combining the root-finding logic with the numerical integration:

$$t_{new}=t_{old}-\frac{\text{Current Arc Length}-\text{Target Distance}}{\text{Instantaneous Speed}}$$

Expressed mathematically with the explicit function call, including the integrand $v$:

$$t_{k+1}=t_k-\frac{\text{AdaptiveSimpson}(v,0,t_k,\text{TOL})-L_{target}}{\parallel {\mathbf{r}}^{\prime }(t_k)\parallel }$$