lec20.mws

L-20 MCS 320 Wednesday 12 October 2005

> restart;

1. Truncated Series

Working with series is a way to symbolically compute approximately.

> tf := taylor(f(x),x=a,10);

tf := series(f(a)+D(f)(a)*(x-a)+1/2*`@@`(D, 2)(f)(a)*(x-a)^2+1/6*`@@`(D, 3)(f)(a)*(x-a)^3+1/24*`@@`(D, 4)(f)(a)*(x-a)^4+1/120*`@@`(D, 5)(f)(a)*(x-a)^5+1/720*`@@`(D, 6)(f)(a)*(x-a)^6+1/5040*`@@`(D, 7)(...tf := series(f(a)+D(f)(a)*(x-a)+1/2*`@@`(D, 2)(f)(a)*(x-a)^2+1/6*`@@`(D, 3)(f)(a)*(x-a)^3+1/24*`@@`(D, 4)(f)(a)*(x-a)^4+1/120*`@@`(D, 5)(f)(a)*(x-a)^5+1/720*`@@`(D, 6)(f)(a)*(x-a)^6+1/5040*`@@`(D, 7)(...

> whattype(tf);

series

A series is different from a polynomial, in the sense that a series goes on for ever, whereas a polynomial has a bounded degree.

> dismantle(tf); # Maple has a separate type for serie

SERIES(24)
  SUM(5)
     NAME(4): x
     INTPOS(2): 1
     NAME(4): a
     INTNEG(2): -1
  FUNCTION(3)
     NAME(4): f
     EXPSEQ(2)
        NAME(4): a
  [0]
  FUNCTION(3)
     FUNCTION(3)
        NAME(4): D #[protected, _syslib]
        EXPSEQ(2)
           NAME(4): f
     EXPSEQ(2)
        NAME(4): a
  [1]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 2
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/2
        INTPOS(2): 1
        INTPOS(2): 2
  [2]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 3
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/6
        INTPOS(2): 1
        INTPOS(2): 6
  [3]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 4
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/24
        INTPOS(2): 1
        INTPOS(2): 24
  [4]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 5
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/120
        INTPOS(2): 1
        INTPOS(2): 120
  [5]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 6
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/720
        INTPOS(2): 1
        INTPOS(2): 720
  [6]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 7
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/5040
        INTPOS(2): 1
        INTPOS(2): 5040
  [7]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 8
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/40320
        INTPOS(2): 1
        INTPOS(2): 40320
  [8]
  SUM(3)
     FUNCTION(3)
        FUNCTION(3)
           FUNCTION(3)
              NAME(4): `@@` #[protected, _syslib]
              EXPSEQ(3)
                 NAME(4): D #[protected, _syslib]
                 INTPOS(2): 9
           EXPSEQ(2)
              NAME(4): f
        EXPSEQ(2)
           NAME(4): a
     RATIONAL(3): 1/362880
        INTPOS(2): 1
        INTPOS(2): 362880
  [9]
  FUNCTION(3)
     NAME(4): O #[protected]
     EXPSEQ(2)
        INTPOS(2): 1
  [10]

> s

With Taylor series we can define functions:

> g := z/(1-z-z^2); # a generating functio

g := z/(1-z-z^2)

> n

> fib := taylor(g,z=0,10)

fib := series(z+z^2+2*z^3+3*z^4+5*z^5+8*z^6+13*z^7+21*z^8+34*z^9+O(z^10),z,10)

> ;

The coefficients of the Taylor series of g define the Fibonacci numbers.

> coeff(fib,z^6);

8

> coeff(fib,z^30);

0

The coeff command gives us the coefficient with a certain degree, just as a polynomial.

But the series "fib" is a truncated series, which goes here only to degree 9.  This the reason why we get 0 back as coefficient with z^30.

> coeftayl(g,z=0,30);

832040

The coeftayl computes the coefficients of a Taylor series, without first having to expand a Taylor series up to this given degree.  This is more efficient than first computing the full Taylor series and then selecting the desired coefficient.

> fibfun := n -> coeftayl(g,z=0,n):

> fibfun(30);

832040

This illustrates how we may define new functions from the Taylor expansions of generating functions.

There is also the multivariate Taylor approximation:

> mtaylor(f(x,y),[x=a,y=b],4);

f(a, b)+D[1](f)(a, b)*(x-a)+D[2](f)(a, b)*(y-b)+1/2*D[1, 1](f)(a, b)*(x-a)^2+(x-a)*D[1, 2](f)(a, b)*(y-b)+1/2*D[2, 2](f)(a, b)*(y-b)^2+1/6*(x-a)^3*D[1, 1, 1](f)(a, b)+1/2*(x-a)^2*D[1, 1, 2](f)(a, b)*(...f(a, b)+D[1](f)(a, b)*(x-a)+D[2](f)(a, b)*(y-b)+1/2*D[1, 1](f)(a, b)*(x-a)^2+(x-a)*D[1, 2](f)(a, b)*(y-b)+1/2*D[2, 2](f)(a, b)*(y-b)^2+1/6*(x-a)^3*D[1, 1, 1](f)(a, b)+1/2*(x-a)^2*D[1, 1, 2](f)(a, b)*(...

The above command gave us a fourth order Taylor approximation of any function at (a,b).

2. Approximations of Functions

We may approximate a complicated function by a truncated Taylor series.

As example we take a function defined by an integral:

> intfun := int(exp(sin(x)),x=0..t);

intfun := int(exp(sin(x)), x = (0 .. t))

The independent variable is the t.  The function gives us the area under exp(sin(x)), for x ranging between 0 and t.

To approximate this function, we may expand the integral as a series around t = 2:

> serfun := series(intfun,t=2,3);

serfun := series(int(exp(sin(x)), x = (0 .. 2))+exp(sin(2))*(t-2)+1/2*exp(sin(2))*cos(2)*(t-2)^2+O(t-2^3),t-2,3)

This will lead to a second-degree polynomial in t.  Before we can use this polynomial function, we must convert the series into a polynomial.

> polfun := convert(serfun,`polynom`); # important intermediate ste

polfun := int(exp(sin(x)), x = (0 .. 2))+exp(sin(2))*(t-2)+1/2*exp(sin(2))*cos(2)*(t-2)^2

> p

Only now, we can convert the expression for polfun into a polynomial function.

> pf := unapply(polfun,t);

pf := proc (t) options operator, arrow; int(exp(sin(x)), x = (0 .. 2))+exp(sin(2))*(t-2)+1/2*exp(sin(2))*cos(2)*(t-2)^2 end proc

This is still a symbolic polynomial:

> pf(2.34);

int(exp(sin(x)), x = (0 .. 2))+.34*exp(sin(2))+0.5780000000e-1*exp(sin(2))*cos(2)

> myfun := (n,t) -> evalf(pf(t),n);

myfun := proc (n, t) options operator, arrow; evalf(pf(t), n) end proc

> myfun(30,2.34);

5.02089342977595071797647716681

We just computed a 30-digit approximation of the function at t = 2.34.

> evalf(subs(t=2.34,intfun),30); # exact value with 30 digits

5.01108183462151253570477531240

> (2.34-2)^3; # evaluation of the truncated term O((t-2)^3)

0.39304e-1

We can replace complicated functions by second-degree polynomials, but their accuracy is limited.

3. Formal Power Series

We can define power series formally, without every having to worry about their convergence.

> with(powseries);

[compose, evalpow, inverse, multconst, multiply, negative, powadd, powcos, powcreate, powdiff, powexp, powint, powlog, powpoly, powsin, powsolve, powsqrt, quotient, reversion, subtract, template, tpsf...[compose, evalpow, inverse, multconst, multiply, negative, powadd, powcos, powcreate, powdiff, powexp, powint, powlog, powpoly, powsin, powsolve, powsqrt, quotient, reversion, subtract, template, tpsf...[compose, evalpow, inverse, multconst, multiply, negative, powadd, powcos, powcreate, powdiff, powexp, powint, powlog, powpoly, powsin, powsolve, powsqrt, quotient, reversion, subtract, template, tpsf...

> powcreate(defexp(n) = 1/n!);

This command defines the n-th term, with the exponent x^n for the exponential function.

We can expand to a series:

> serexp := tpsform(defexp,x,10);

serexp := series(1+x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+1/720*x^6+1/5040*x^7+1/40320*x^8+1/362880*x^9+O(x^10),x,10)

This is one of the definitions for the exponential function.

As before, we may convert this series into a function.

> polexp := convert(serexp,`polynom`);

polexp := 1+x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+1/720*x^6+1/5040*x^7+1/40320*x^8+1/362880*x^9

The error between the actual series (and thus the exp() function) and the polynomial is of order x^10.

> funexp := unapply(polexp,x);

funexp := proc (x) options operator, arrow; 1+x+1/2*x^2+1/6*x^3+1/24*x^4+1/120*x^5+1/720*x^6+1/5040*x^7+1/40320*x^8+1/362880*x^9 end proc

> funexp(1.1); exp(1.1);

3.004165230

3.004166024

> truncation_error := .1^10;

truncation_error := 0.1e-9

A more interesting series is a series to compute Pi:

> powcreate(f(n)=(-1)^n/(2*n+1));

> fseries := tpsform(f,n,7);

fseries := series(1-1/3*n+1/5*n^2-1/7*n^3+1/9*n^4-1/11*n^5+1/13*n^6+O(n^7),n,7)

What is interesting here is the sum of the series:

4. Limits

> s := Sum('(-1)^n/(2*n+1)','n'=0..infinity);

s := Sum((-1)^n/(2*n+1), n = (0 .. infinity))

> value(s);

1/4*Pi

> s1 := k -> sum('(-1)^n/(2*n+1)','n'=0..k);

s1 := proc (k) options operator, arrow; sum('(-1)^n/(2*n+1)', 'n' = (0 .. k)) end proc

We can see how fast this series approximates Pi/4 by evaluating this function:

> for i from 1 to 20 do

>  evalf(s1(i),30);

> end do;

.666666666666666666666666666667

.866666666666666666666666666667

.723809523809523809523809523810

.834920634920634920634920634921

.744011544011544011544011544012

.820934620934620934620934620935

.754267954267954267954267954268

.813091483679718973836620895444

.760459904732350552783989316497

.808078952351398171831608364116

.764600691481832954440304016290

.804600691481832954440304016290

.767563654444795917403266979253

.802046413065485572575680772356

.769788348549356540317616256227

.800091378852386843347919286530

.771519950280958271919347857959

.798546977307985298946374884986

.772905951666959657920733859345

.797296195569398682310977761784