How to approximate $\sin(x)$ using Padé approximation?

Question

I need to write a function for $\sin(x)$ using Padé approximation.

I found here a formula, but I don't know how to achive this.

Thanks in advance.

Answer 1

Since you know Taylor expansions, let us try a simple way.

You want to make $$f(x)=\frac{\sum_{i=0}^n a_ix^i } {1+ \sum_{i=1}^m b_ix^i }$$ where $m,n$ are predefined degrees. So, write $$f(x)\left({1+ \sum_{i=1}^m b_ix^i }\right)=\sum_{i=0}^n a_ix^i $$ Use the Taylor expansion of $f(x)$ and identify coefficients.

Let us try with $f(x)=e^x$, $m=2$, $n=3$. So, we have $$\left(1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+\frac{x^5}{120}\right)\left(1+b_1x+b_2x^2+b_3x^3\right)=a_0+a_1x+a_2x^2$$ Expanding, grouping terms and setting coefficients equal to $0$ will lead to the following equations $$a_0=1$$ $$1-a_1+b_1=0$$ $$-a_2+b_1+b_2+\frac{1}{2}=0$$ $$\frac{b_1}{2}+b_2+b_3+\frac{1}{6}=0$$ $$\frac{b_1}{6}+\frac{b_2}{2}+b_3+\frac{1}{24}=0$$ $$5 b_1+20 b_2+60 b_3+1=0$$ from which $$a_0=1 \qquad a_1=\frac 25\qquad a_2=\frac 1{20}$$ $$b_1=-\frac35 \qquad b_2=\frac 3{20}\qquad b_3=-\frac 1{60}$$ and then $$\frac{1+\frac{2}{5}x+\frac{1}{20}x^2 } {1-\frac{3 }{5}x+\frac{3 }{20}x^2-\frac{1}{60}x^3 }$$ is the Padé approximant.

Answer 2

The relation $p_k/q_k=\sin(x)+O(x^N)$ or $p_k=s_N·q_k+O(x^N)$ can be interpreted as one stage in the extended euclidean algorithm starting with $r_0=x^N$, $r_1=s_n$ and propagating the Bezout identities $r_k=u_ks_N+v_kx^N$, $u_0=0, u_1=1$, $v_0=1, v_1=1$ and identifying $p_k=r_k$, $q_k=u_k$. Implementing in Magma CAS which you can try out with the online calculator

PS<z>:=PowerSeriesRing(Rationals());
P<x>:=PolynomialRing(Rationals());
PadeSin:=function(N)
    sN := Truncate(Sin(z+O(z^N)));
    sN := Evaluate(sN,x);
    r0 := x^N; u0:=0;
    r1 := sN;  u1:=1;
    for k in [2..N] do
        if Degree(r1) lt 2 then break; end if;
        q1,r2 := Quotrem(r0,r1); // r2 = r0-q1*r1
        u2 := u0-q1*u1;
        a := Coefficient(u2,0);
        printf "k=%o : p/q = (%o)/(%o)\n",k,r2/a,u2/a;
        r0:=r1; r1:=r2;
        u0:=u1; u1:=u2;
    end for;
    return r1,u1;
end function;
for N in [5,7,13] do "\norder N=", N; PadeSin(N); end for;

gives a result of

order N= 5
k=2 : p/q = (x)/(1/6*x^2 + 1)
order N= 7
k=2 : p/q = (-7/60x^3 + x)/(1/20x^2 + 1)
k=3 : p/q = (x)/(7/360x^4 + 1/6x^2 + 1)
order N= 13
k=2 : p/q = (19/19958400x^9 - 17/138600x^7 + 3/440x^5 - 26/165x^3 + x)/(1/110x^2 + 1)
k=3 : p/q = (-121/2268000x^7 + 601/118800x^5 - 241/1650x^3 + x)/(19/118800x^4 + 17/825x^2 + 1)
k=4 : p/q = (12671/4363920x^5 - 2363/18183x^3 + x)/(121/16662240x^6 + 601/872784x^4 + 445/12122x^2 + 1)
k=5 : p/q = (-2555/25146x^3 + x)/(12671/7604150400x^8 + 2363/31683960x^6 + 3787/1508760x^4 + 818/12573x^2 + 1)
k=6 : p/q = (x)/(73/3421440x^10 + 127/604800x^8 + 31/15120x^6 + 7/360x^4 + 1/6*x^2 + 1)

where you find under N=13, k=4 your cited fraction.