Getting the most accurate bezier curve that plots a sine wave I would like to draw a bezier curve that looks the most like a sine wave that has a single wave length of 1000 pixel and an amplitude of 1, which is 159.15 pixels high (wave length / 2).
Here are the parameters of my bezier curve:

*

*p1 = starting point

*p2 = ending point

*h1 = outer handle linked to p1

*h2 = inner handle linked to p2

Since I have set the wave length at 1000 pixels to begin with, I know my p1 and p2:

*

*p1 = (0,0)

*p2 = (1000,0)

Now, I need to find the coordinates of h1 and h2. But I don't know how to get there by calculation.
By trial and error I get approximately this:

*

*h1 = (477,550)

*h2 = (523,-550)

This is not perfect as no bezier would represent perfectly a sine curve but pretty close (see picture below).
My question is: knowing my starting parameters (wave length in pixels, amplitude in pixels, p1 and p2 coordinates), how can I find mathematically h1 and h2?
It seems that a bezier function with 4 coordinates (p1,h1,p2,h2) uses a cubic equation that should look like this:
$$B(t) = (1-t)^3 . p1 + 3(1-t)^2.t.h1 + 3(1-t).t^2.h2+t^3.p2$$
I don't know what to do with this though since the wave length and the amplitude need to be accounted for so that h1 and h2 can be found.
That said, I only need to find h1 since h2 will always have as its x value the wave length minus the x value of h1; and its y value will be the negative of the y value of h1.
Anyway, here is the visual result I came up with by trial and errors:
The white graph is the Desmos Graph accurate rendering of a sin(x).
Once cropped in Photoshop to make it a single wave length and edited to make it 1000 pixels wide, I used bezier (the black graph) in After Effects to replicate it. In pink you can see the values of h1 and h2 I came up with by trials and errors.
Thanks for your help.

 A: Let's assume for simplicity that we want to draw $\sin(x)$ with $x$ from $0$ to $2\pi$, so an amplitude of $1$.
First note that the "closest" curve is not mathematically well-defined. We have to specify what we want with more details. The assumptions I do in this answer are:

*

*We want the endpoints to be at the right positions.

*We want the curve to be 180° rotationally symmetric around the middle of the two endpoints.

*We want the amplitude to be exactly right.

As you said we can use this formula for a cubic Bézier curve, where $t$ goes from  $0$ to $1$:
$$
B(t) = (1-t)^3 \cdot p_1 + 3(1-t)^2  t \cdot h_1 + 3(1-t)  t^2 \cdot h_2 + t^3 \cdot p_2 \tag{1}\label{B}
$$
If we fix the endpoints and the symmetry around the middle, we have only two unknowns left, which I'll call $u$ and $v$.
$$
\begin{align}
p_1 &= (0, 0) \\
p_2 &= (2\pi, 0) \\
h_1 &= (u,v) \\
h_2 &= (2\pi-u,-v)
\end{align}
$$
Determining $v$
The vertical part of the Bézier curve becomes:
$$
\begin{align}
B_y(t) &= 3(1-t)^2 t v + 3(1-t) t^2 (-v) \\
&= v \cdot (6 t^3 - 9t^2+3t)
\end{align}
$$
Now I'll first set $v=1$ and calculate what amplitude we're getting then.
The maximum of $B_y$ can be calculated by setting its derivative
equal to zero.
$$
\begin{align}
y &= 6 t^3 - 9t^2+3t \tag{2}\label{y} \\
\frac{\mathrm{d}y}{\mathrm{d}t} &= 18 t^2 - 18t+3
\end{align}
$$
We can use the quadratic formula with $a=18$, $b=-18$ and $c=3$.
There are two solutions, one is the minimum and the other is the maximum. The maximum is the one with the lower $t$, because the curve goes from $p_1$ at $t=0$ to $p_2$ at $t=1$. So where the quadratic formula has $\pm$, we use $-$.
$$
\frac{\mathrm{d}y}{\mathrm{d}t} = 0
\quad\implies\quad
t = \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac16 (3 - \sqrt{3}) \tag{3}\label{tmax}
$$
When we plug $\eqref{tmax}$ into $\eqref{y}$, we get $y=\frac{1}{2\sqrt{3}}$.
So to get an amplitude of $1$, we set $v=2\sqrt{3}$.
Determining $u$
To determine $u$, we have to make an other assumption.

*

*The $\sin(x)$ function has a slope of 45° (a derivative of -1 or 1) where it crosses zero. If we want that to be correct at the endpoints, we need the angle between $h_1$, $p_1$ and the x-axis also
to be 45°. We get that if we set $u=v$.



*An other option is to just play around with the value of $u$ and see what looks good. Here is $u=3$:



*We could fix the top at $(\pi/2,1)$. The horizontal part of
the Bézier curve is:
$$
  B_x(t) = 3(1-t)^2  t u + 3(1-t)  t^2 (2\pi-u) + t^3 \cdot 2\pi
  $$
Now if we plug in $\eqref{tmax}$ we have the x coordinate of the top as function of $u$.
$$
  B_x \left(\frac16 (3 - \sqrt{3}) \right) =
  \frac{u}{2 \sqrt{3}}-\frac{4 \pi }{3 \sqrt{3}}+\pi = \frac{\pi}{2} \\
  \implies \\
  u = \left(\frac{8}{3}-\sqrt{3}\right) \pi
  $$



*We could minimize for example:

*

*The maximum Euclidean distance between the curve and the sine wave

*The average squared Euclidean distance between the curve and the sine wave, by doing a line integral over the Bézier curve or over the sine wave. (can have different results for differrent paths)

*The maximum vertical distance between the curve and the sine wave

*The average squared vertical distance between the curve and the sine wave, by integrating $x$ from $0$ to $2\pi$.

These are however, difficult to do symbolically, if possible at all. It can be done numerically with the right software.
Scaling
The values we calculated are for an image with $x \in [0, 2\pi]$ and $y \in [-1, 1]$. We can scale this to any size. Note that in an image, the y-axis is reversed: higher values are lower. If we want to make a sine wave with a width of 1000 pixels, we have to transform the four points:
$$
\begin{align}
x &\quad\to\quad x \cdot \frac{1000}{2\pi} \\
y &\quad\to\quad (1-y) \cdot \frac{1000}{2\pi}
\end{align}
$$
A: So thanks to Paul, I made it and it works great.
Here below is the code I scripted within After Effects and the comments explain how I got there thanks to Paul's help:
//-----------------------------------------------------------//

    // user entry (later: make them variable inputs)

    var wavelength = 1000;
    var amplitude  = 1;

//-----------------------------------------------------------//

    // These values are trivial but necessary for the path function

    const handle_p1  = [0,0];
    const handle_p2  = [0,0];
    const point1     = [0,0];

    var point2       = [wavelength,0];
    var amp          = amplitude * wavelength / (2*Math.PI);

//------------------------------------------------------//

    // handle1 = [u,v], handle2 = [2-u,-v]
    // The bezier with 4 params (point1,handle1,handle2,point2) is cubic:
    // (t) = (1−t)^3 ⋅ p1 + 3(1−t)^2 . t ⋅ h1 + 3(1−t) . t^2 ⋅ h2 + t^3 ⋅ p2
    // If we plugin what we have in terms of y values, we get:
    // (t) = v(6t^3 - 9t^2 + 3t)
    // Let's have v=1 => (t) = 6t^3 - 9t^2 + 3t and optimize this:
    // (t)' = 18t^2 - 18t + 3 = 0
    // discriminant = 108 => roots = (-b ± √108))/2a = (3±√3)/6
    // let's plug those in (t) = 6t^3 - 9t^2 + 3t => y = ± 1/2√3
    // this is y (amplitude) if v=1 but we need v if y=1 => v = 2√3
    // And now we want it if y equals to "amp" => 2√3 * amp

    var v = 2*Math.sqrt(3) * amp; 

//------------------------------------------------------//


    // Find the u value of the Bezier based on the t value found previously
    // ()=3(1−)^2.+3(1−)^2.(2−)+3(2)
    // ((3±√3)/6)= u/(2√3) - 4/(3√3) +  = /2
    // =(8/3−√3)
    // =(8/3−√3). (wavelength/2)

    var u =  (8/3 - Math.sqrt(3)) * wavelength/2; 

//------------------------------------------------------//

    var handle1     = [u,-v];
    var handle2     = [-u,v];

//------------------------------------------------------//

    var points      = [point1, point2];
    var inTangents  = [handle_p2, handle2];
    var outTangents = [handle1, handle_p1];

//------------------------------------------------------//

    createPath(points, inTangents, outTangents, false);

//------------------------------------------------------//

For those who wonder why I use:

*

*handle1   = [u,-v];

*handle2   = [-u,v];

Instead of:

*

*handle1   = [u,v];

*handle2   = [wavelength-u,-v];

It's just because the y coordinates are inverted when programming (v becomes -v) and the x values of the handles are relative to their related points not absolutely positioned against the origin.
Here are various results with various wave length and amplitudes and it nails it all the time.
Thank you all !

