### Quadratic Interpolation

Given a smooth, well-behaved function $f(t)$ that is known to have no more than one root within the range $t_a \le t \le t_b$, we want to find an iterative procedure to quickly narrow in on a value $t_0$ in that range such that $f(t_0) \approx 0$. The function $f$ may be expensive to compute, so we would like to evaluate it as few times as possible. We will use quadratic interpolation, where we approximate $f$ as a parabola.

First find the midpoint time $t_m$ between $t_a$ and $t_b$:

\begin{align}
\\
t_m = \frac {t_b + t_a} {2}
\end{align}

Define a parameter $x$ that has the linear range $-1 \le x \le +1$ over the search interval:

\begin{align}
\\
x = \frac {t - t_m} {t_b - t_m}
\end{align}

A generic parabolic function over this range is

\begin{align}
\\
p(x) = Q x^2 + R x + S
\end{align}

Evaluate $f$ at the left, middle, and right points. These correspond to $x=-1$, $x=0$, and $x=+1$.

\begin{align}
\\
f_a &= f(t_a) \\
f_m &= f(t_m) \\
f_b &= f(t_b)
\end{align}

Now we can relate these 3 values to a unique parabola that passes exactly through the points $(t_a, f_a)$, $(t_m, f_m)$, and $(t_b, f_b)$ as

\begin{align}
\\
p(-1) &= f_a = Q - R + S \\
p(0)  &= f_m = S\\
p(+1) &= f_b = Q + R + S 
\end{align}

Solving this system of equations in terms of the known values $f_a$, $f_m$, and $f_b$, we obtain

\begin{align}
\\
Q &= \frac{f_b + f_a}{2} - f_m \\
R &= \frac{f_b - f_a}{2} \\
S &= f_m
\end{align}

Let's test these formulas:


```python
import math

class Parabola:
    def __init__(self, f, ta, tb):
        tm = (tb + ta)/2
        fa = f(ta)
        fb = f(tb)
        fm = f(tm)
        self.Q = (fb + fa)/2 - fm
        self.R = (fb - fa)/2
        self.S = fm
        
    def Eval(self, x):
        return self.Q*x*x + self.R*x + self.S
    
    def Slope(self, x):
        return 2*self.Q*x + self.R

def Func(t):
    return math.cos(t) - t

ta = 0.6
tb = 0.8
tm = (ta + tb)/2
p = Parabola(Func, ta, tb)
print('Q={:0.6f}, R={:0.6f}, S={:0.6f}'.format(p.Q, p.R, p.S))
print('p(-1)={:0.6f}, f(ta)={:0.6f}'.format(p.Eval(-1), Func(ta)))
print('p( 0)={:0.6f}, f(tm)={:0.6f}'.format(p.Eval( 0), Func(tm)))
print('p(+1)={:0.6f}, f(tb)={:0.6f}'.format(p.Eval(+1), Func(tb)))

```

    Q=-0.003821, R=-0.164314, S=0.064842
    p(-1)=0.225336, f(ta)=0.225336
    p( 0)=0.064842, f(tm)=0.064842
    p(+1)=-0.103293, f(tb)=-0.103293


Now that we have a parabola that passes through the 3 points, let's solve for the value of $x$ such that $p(x)=0$. The quadratic formula gives us two values of $x$:

\begin{align}
\\
x = \frac {-R \pm \sqrt{R^2 - 4QS}} {2Q}
\end{align}

For our purposes, the solution must have a real value and must lie in the range $-1 \le x \le +1$. We expect exactly one such solution. If there are zero or two such values of $x$, our search will return a null value.


```python
def Solve(p):
    if p.Q == 0:
        # This is not a parabola but a straight line
        # p(x) = Rx + S = 0
        # Therefore x = -S/R
        if p.R == 0:
            return None
        x = -p.S / p.R
        if -1 <= x <= +1:
            return x
        return None
    
    # Let u = the quantity inside the square root.
    u = p.R*p.R - 4*p.Q*p.S
    if u <= 0:
        # If u<0, then both solutions are complex-valued.
        # If u=0, the parabola is tangent to the axis, not crossing it.
        return None
    
    ru = math.sqrt(u)
    x1 = (-p.R + ru) / (2 * p.Q)
    x2 = (-p.R - ru) / (2 * p.Q)
    
    if -1 <= x1 <= +1:
        if -1 <= x2 <= +1:
            return None
        return x1
    
    if -1 <= x2 <= +1:
        return x2
    
    return None
```


```python
x = Solve(p)
print(x)
```

    0.39106619198014525


Convert the value $x$ back into a time value $t$.


```python
t = tm + x*(tb-tm)
print(t)
```

    0.7391066191980145


See how close this is to the correct solution.


```python
f = Func(t)
print(f)
```

    -3.595936996025895e-05


We want to iterate and do another parabolic interpolation, but the problem is we don't know where to choose a new time range $\bar{t}_a \le t \le \bar{t}_b$. We want this new range to be as small as possible while definitely bracketing the place where $f(t)$ crosses the $t$-axis. If we knew where that was, we would have already solved the problem and wouldn't need to do this!

We are actually more interested in finding $t_0$ within a sufficiently small window $\Delta t$ that definitely brackets the root, than we are in finding a very small value $\left|{f(t)}\right|$.

We can keep track of the bounding values of $t$ that most closely bracket where the root must be, and iterate with quadratic interpolation around each successive approximation of $t_0$. The trick might be to just create another parabola through the $t$ value we just found by picking a new value of $t_a$ and $t_b$ around it, without worrying about whether the parabola crosses the $t$-axis between the new $t_a$ and $t_b$, so long as it stays inside the original range.


```python

```

b7edc4ef8f295b1fe92acf069438ad37e8160c49
ipynb
Jupyter Notebook
theory/quadratic_interpolation.ipynb
matheo/astronomy
3a1d4ea47a0c04d83bd8ede43dc564e956e999fe
MIT
theory/quadratic_interpolation.ipynb
matheo/astronomy
3a1d4ea47a0c04d83bd8ede43dc564e956e999fe
MIT
theory/quadratic_interpolation.ipynb
matheo/astronomy
3a1d4ea47a0c04d83bd8ede43dc564e956e999fe
MIT
Qwen/Qwen-72B
1. YES
2. YES