In this worksheet we use rational polynomials to approximate functions

1. Approximation with Rational Functions

For example, let us consider the function sin(x).

> restart;

> n := 3: # degree numerator

> m := 3: # degree denominator

> N := n+m: # number of terms in Maclaurin series

We start with a Maclaurin series:

> ms := taylor(sin(x),x=0,N+1);

> pms := convert(ms,polynom);

> q := (sum(a['i']*x^'i','i'=0..n))/(1 + sum(b['i']*x^'i','i'=1..m));

> err := pms*denom(q) - numer(q);

> eqs := {}:

> for i from 0 to N do

> eqs := {op(eqs),coeff(err,x,i) = 0};

> end do:

> eqs;

> vars := {seq(a[k],k=0..n),seq(b[l],l=1..m)};

> sols := solve(eqs,vars);

> assign(sols);

> q;

Maple has a built-in function to construct Pade approximations:

> r4 := numapprox[pade](sin(x),x,[n,m]);

We see that the error of the rational approximation is less than that of the Maclaurin series:

> plot([sin(x),r4,pms],x=-Pi..Pi);

2. Continued Fractions

To obtain a form of the approximation that is efficient to evaluate, we use continued fractions.

> r4;

> codegen[cost](r4);

> cr4 := convert(r4,`confrac`,x);

> codegen[cost](cr4);

Thus we traded in 4 multiplications for one extra division.

We obtain the continued fraction representation by repeated polynomial division:

> a := numer(r4);

> b := denom(r4);

> rest1 := rem(a,b,x,'quot1');

> quot1;

Sanity check:

> a = expand(b*quot1 + rest1);

Thus we can write the quotient as

> quot1 + rest1/b;

Because deg(rest) < deg(b), we work with the equivalent 1/(b/rest1).

> rest2 := rem(b,rest1,x,'quot2');

Sanity check:

> b = expand(rest1*quot2 + rest2);

Thus we have the expression

> fr4 := quot1 + 1/(quot2 + rest2/rest1);

The expression we obtained is equivalent with what Maple gave us:

> normal(fr4 - cr4);

> codegen[cost](fr4);

We can save one multiplication by multiplying the quotient in the second term of fr4 by 200/3

> dd := op([2,1],fr4)*200/3;

> nfr4 := op(1,fr4) + (200/3)/dd;

> codegen[cost](nfr4);