NumericalEvaluation.ml

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

INCLUDE "preamble.ml"

let title _ = <:text<Numerical Evaluation>>

let_service NumericalEvaluation
  (def : diffeq maple)
  (var : name maple)
  (rep : any maple)
  (p_hack : any maple)
  (v_hack : any maple)
  (path : string = "1/4+1/4*i")
  (precision : int = Constants.default_precision) :
  DC.sec_entities * unit with { title = title } =

  <:unit< i := I: >> ;
  let btest =
    <:bool<
      try type(eval(parse($(str:path))),
            'Or'('complex'('numeric'), 'list'('complex'('numeric'))))
      catch "incorrect syntax": false
      end:
    >> in
  let spath = <:symb< eval(parse($(str:path))) >> in
  (* TODO: Passing lists as arguments to a service causes an error. *)
  (* As a hack, we converted the lists params and values to any maple. *)
  (* Now we have to recover their original values. *)
  let params = List.map Obj.magic (CommonTools.symb_list_of_symb p_hack)
  and values =
    List.map
      (fun a -> <:string< $(a) >>)
      (CommonTools.symb_list_of_symb v_hack) in
  let sf_name_as_math =
    CommonTools.subs_expr rep (var :: params) (path :: values) in
  let par1 =
    if not btest then
      DC.warning
        <:par< Unable to evaluate <:imath< $(str:sf_name_as_math) >>:
               $(str:path) is not a valid point or path. >>
    else
      let res =
        <<
          try gfun:-NumGfun:-analytic_continuation(
                  $(def), y(x), $(spath), $(int:precision))
          catch "unable to perform":
            __MyError("the analytic continuation path goes through "
              "a singular point of the differential equation")
          catch "invalid input: expected evalf[] index to be of type":
            __MyError("unable to compute reasonable error bounds")
          catch "step", "invalid path", "evaluation point outside",
                "unable to compute a suitable truncation order":
            __MyError(StringTools[FormatMessage](lastexception[2..-1]))
          end:
        >> in
      <:unit< unassign('i'): >> ;
      if <:bool< hastype($(res), 'specfunc'('string', __MyError)) >> then
        let err = <:string< op($(res)) >> in
        <:par< Error while trying to compute <:imath< $(str:sf_name_as_math) >>:
               $(str:err). >>
      else
        <:par<<:dmath< $(str:sf_name_as_math) \approx $(symb:res) . >>>>

  in
  let par2 =
    (* * * taken from singularExpansions *)
    let eqn = << op(remove(type, $(def), equation)) >> in
    let diff_order = <:int< nops(indets($(eqn), function)) - 1 >> in
    <:unit< cf := coeff($(eqn), diff(y(x), [x $ $(int:diff_order)])): >> ;
    let sing = << sort(convert( {solve(cf)} minus {0}, list)) >> in
    (* * *)
    let rad = << min(op(map(abs, $(sing)))) >> in

    DC.note (
      <:par< (Below, $(t_ent:DC.code "path") may be either a point
             <:imath< z >> >> @:@
      (if <:bool< $(rad) < infinity >> then
        <:par<  with <:imath< |z| < $(symb:rad) >> >>
      else
        <:par<  >>) @:@
      <:par< or a broken-line path <:imath< [z_1, z_2, \ldots, z_n] >> along
             which to perform analytic continuation of the solution of the
             defining differential equation. Each <:imath< z_k >> should be of
             the form $(t_ent:DC.code "x + y*i").) >>
    )
  in
  (DC.section (title ()) (par1 @@@ par2), ())

Generated by GNU Enscript 1.6.5.90.