Chebyshev.ml

(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)
(* DDMF is distributed under CeCILL-B license. *)

INCLUDE "preamble.ml"

(* Give the closed form of the Chebyshev expansion if this form exist. *)
let closed_form notation def yofz uofk =
  <:unit<
    unassign('y', 'x', 'u', 'n'):
    gfun:-getname($(yofz), 'y', 'x'):
    gfun:-getname($(uofk), 'u', 'n'):
    general_term := Chebyshev[FirstTerms][chebGeneralForm](
                      $(def), y, x, u, n):
    nb_terms := nops(general_term:-general)
  >> ;
  if <:bool< general_term:-isHyper >> then
    (* If the first coefficients are not the initials conditions *)
    (* of the recurrence for example arccos(x). *)
    let s1 = <<
      1/2*coeff(general_term[poly], x, 0) * T[0](x) +
      add(coeff(general_term[poly], x, i) * T[i](x),
          i=1..degree(general_term[poly], x))
    >>
    and s2 =
      (* We deal only the case of the first coefficient *)
      (* (we multiply this coefficient by 1/2).         *)
      if <:bool< eval(general_term:-general[1]:-period, n=0) = 0 >>
      then
        (* We use simplify to find special functions of hypergeometric form. *)
        (* For example the Bessel function *)
        let f1 =
          <<
            simplify(subs(n=0, 1/2*general_term:-general[1]:-gterm)) * T[0](x)
          >>
        and f2 =
          if <:bool< general_term:-general[1]:-gterm <> 0 >> then
            <<
              if (type(general_term:-general[1]:-gterm, `*`) and
                   type(op(1, general_term:-general[1]:-gterm),
                        'complex'('negative'))) then
                  signu := -1
               else signu := 1
               fi;
               signu*Sum( signu*simplify(general_term:-general[1]:-gterm *
                            T[general_term:-general[1]:-period](x)),
                   n = 1..infinity )
            >>
          else << 0 >>
        and f3 =
          if <:bool<
               seq(l[gterm],
                   l = subsop(1 = NULL, general_term:-general)) <>
               seq(0, i=2..nb_terms)
             >>
          then
            <<
              f3maple := 0;
              for l in subsop(1 = NULL, general_term:-general) do
                if (type(l:-gterm, `*`) and
                     type(op(1, l:-gterm),
                        'complex'('negative'))) then
                    signu := -1
                 else signu := 1
                 fi;
                 f3maple := f3maple + signu*Sum(signu*simplify(l[gterm] *
                                T[l[period]](x)), n = 0..infinity)
              od;
              f3maple
            >>
          else << 0 >>
        in << $(f1) + $(f2) + $(f3) >>
      else
        <<
          add(Sum(simplify(l[gterm] * T[l[period]](x)),
                  n = 0..infinity),
              l = general_term:-general)
        >>
    in
    let expn = << subs( Sum(0, n=0..infinity) = 0, $(s1) + $(s2) ) >> in
    (* FIXME: Avoid this Maple hack by using more LaTeX for presentation. *)
    (* That is: don't use Maple's T[...](x) and Sum, but T_{...} and \sum. *)
    <:unit<
      map(sort, map(op,
        select(type, map2(op, 0, indets($(expn), function)), indexed) ))
    >> ;
    Some
      <:par<Chebyshev expansion: <:dmath< $(str: notation) = $(symb: expn). >>>>
  (* Don't display anything if there is no closed form. *)
  else None


let has_a_recurrence def yofz =
  <:unit<
    unassign('y', 'x'):
    gfun:-getname($(yofz), 'y', 'x')
  >> ;
  (* If we don't have the initial conditions of def, the function is not *)
  (* defined over the segment [-1,1]. *)
  (* The set containing the differential equation has 1 element. *)
  (* We want the homogeneous differential equation only. *)
  let def2 = << gfun[diffeqtohomdiffeq]($(def), y(x)) >> in
  def2, <:bool< nops($(def2)) > 1 and type($(def2), set) >>


let prepare_plot def2 rep notation proc =
  <:unit<
    plotfcn := plots[display](plot($(proc), -2..2, numpoints=200,
                                   color=blue,thickness=3)) ;
    xrange, yrange := DynaMoW:-GetViewPort(plotfcn) ;
    series($(rep), x)
  >> ;
  (* FIXME: Justify the choice of the order. *)
  let f_terms = NumericalChebyshev.obj ((def2, notation), 6) in
  f_terms

let_service ChebyshevPlot1
  (def2 : any maple)
  (rep : any maple)
  (notation : string)
  (proc : any maple) :
  DynaMoW.Services.svg * unit =
  let f_terms = prepare_plot def2 rep notation proc in
  let p =
    <<
      plots[display](
        [ plotfcn,
          plot({seq(1/2 * $(f_terms)[0] +
                      add($(f_terms)[l]*orthopoly[T](l,x),l=1..i),
                    i = 0..5)},
               x = -2..2, numpoints=200) ],
        view=yrange)
    >> in
  (<:string< DynaMoW:-PlotToSVG($(p), -2..2, yrange) >>, ())

let_service ChebyshevPlot2
  (def2 : any maple)
  (rep : any maple)
  (notation : string)
  (proc : any maple) :
  DynaMoW.Services.svg * unit =
  let f_terms = prepare_plot def2 rep notation proc in
  let p =
    <<
      plots[display](
        [ plot({seq(1/2*$(f_terms)[0] +
                      add($(f_terms)[l] * orthopoly[T](l,x), l = 1..i) -
                      $(rep),
                    i = 2..5)},
               x = -1..1, numpoints=200) ])
    >> in
  (<:string< DynaMoW:-PlotToSVG($(p)) >>, ())

(* Compute and print the Chebyshev section of the DDMF.  Compute *)
(* the recurrence, the closed form and the approximation of coefficients. *)

let recurrence def yofz rep notation proc =
  let def2, ok_to_continue = has_a_recurrence def yofz in

  if not ok_to_continue then
    failwith
      ("no Chebyshev recurrence available \
       (has_a_recurrence should be tested before calling this service)") ;

  let rec_cheby =
    <<
      Chebyshev:-diffeqtorecchebyshev(
        remove(type, $(def2), equation), y(x), c(n))
    >>
  and cform = closed_form notation def2 << y(x) >> << u(n) >> in
  let unordered_list =
    [ (match cform with Some par -> par | None -> DC.ent_null) ;
      DC.inline_service (NumericalChebyshev.descr (def2, notation) None) ;
      <:par<\
        The coefficients <:imath<c_n>> in the \
        $(t_ent:Glossary.g "Chebyshev expansion") \
        <:imath< $(str: notation) = \sum_{n=0}^\infty c_nT_n(x) >> \
        satisfy the recurrence
        <:dmath<<< op($(rec_cheby)) = 0 >> . >>
      >> ] in

  let plot1 = DC.plot (ChebyshevPlot1.descr (def2, rep, notation, proc) None)
  and plot2 = DC.plot (ChebyshevPlot2.descr (def2, rep, notation, proc) None) in
  let par1 =
    <:par<
      Approximations by successive truncations of the Chebyshev series on
      <:imath<[-2,2]>>: $(plot1)
    >>
  and par2 =
    <:par<
      Errors of approximation by successive truncations of the Chebyshev
      series on <:imath<[-1,1]>>: $(plot2)
    >> in
  DC.unordered_list (unordered_list @ [ par1 ; par2 ])

let title _ = <:text<Chebyshev Expansion over <:imath<[-1,1]>>>>

let_service Chebyshev
  (def : diffeq maple)
  (yofz : fcall maple)
  (rep : any maple)
  (notation : string)
  (proc : any maple) :
  DC.sec_entities * unit with { title = title } =
  DC.section
    (title _req_params)
    (* TODO: instead of this try, one should fix the bugs! *)
    (try (recurrence def yofz rep notation proc)
     with _ -> DC.warning <:par<Some error occurred.>>), ()

Generated by GNU Enscript 1.6.5.90.