# Chebyshev.ml

(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)

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
<<
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.