Axel Vogt

5936 Reputation

20 Badges

20 years, 251 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

At least for the example values one can use fdiff (for cross checks) 
'eval(H, [a=convert(v, rational),b=convert(w, rational)])';
fdiff(%, y=1e-9);
fdiff(%%, y=0);
    -621.372 .... (100 Digits)
Is it true that lHospital can not be used for H ?

For your actual task: I can not follow your code (too technical for
me), but may be you can use (for your last 'example')

multiseries((exp(y)-1-y)/y^2, y, 8, 'exact_order');

convert(Z, hypergeom): simplify(%); returns zero

Usually factorization into linear or quadratic facors is meant.

But that is equivalent to know the real or complex-paired zeros.

If you additionally allow for parameters a[j], then it means that
you take an arbitrary polynomial and ask for its zeroes.

But you can not expect that Maple will do it. Or do you?

ShellExecute := define_external(
  'ShellExecuteA',
   hwnd::integer[4],
   lpOperation::string,
   lpFile::string,
   lpParameters::string,
   lpDirectory::string,
   nShowCmd::integer[4],
   'RETURN'::integer[4] ,
   LIB="C:\\WINDOWS\\SYSTEM32\\shell32.dll"): # DLL for Win XP

OpenBrowser:=proc(someURL::string :="http://")
  ShellExecute(0, "", someURL, "", "", 1); NULL;
end proc:
# may be "Open" should be used instead of the first, empty string

Then the following opens a new browser window:

OpenBrowser();

and it is possible to do that with a specific URL:

OpenBrowser("https://panopticlick.eff.org/index.php?action=log");


Remark: this can be risky for the user, if s/he has not taken care against
Java programs (which is the default for Firefox, unfortunately allowing for
javascript opens that door).

Had a look. This expression has length ~ 100 pages and depends on 3 parameters
without any assumptions (may be complex constants).

Why do you expect a limit exits and a CAS can find it in reasonable time?

For your underlying question: with(MultiSeries) sometimes is a way.

But if the expression is very lengthy or complicated it may be, that
the system hangs up in intermediate automatic simplifications.

Just try 'simplify' or 'normalize' on your example.

I would not use it that way:
> p := proc() local a,b; a:=2; b:=3; a/b end proc;
> %();

         p := proc() local a, b; a := 2; b := 3; a/b end proc

                                 2/3

> eval(eval(p), [1=5, -1 = 17]);
> %();

           proc() local a, b; a := 2; b := 3; a/b end proc

                                 2/3

> subs(1=5, -1 = 17, eval(p));
> %();

         proc() local a, b; a := 2; b := 3; a^5*b^17 end proc

                              4132485216

> subs(1=Fritz, -1 = TheCat, eval(p));
> %();

     proc() local a, b; a := 2; b := 3; a^Fritz*b^TheCat end proc

                             Fritz  TheCat
                            2      3



I think, that does not matter to take conjugates (conjugate again).

It only 'indicates', that one is looking for complex solutions.

BTW: setting z1 = - z2 reduces the task to a single equation, no?

One gets more than 1 solution, at least in the real case (but forgot
too much about the p function & periodics)

Digits:=15;
[conjugate(1/Pi*I*f1), conjugate(1/Pi*I*f2)];
eval(%, z2 = -z1);
eval(%, z1=y*I);
tst:=eval(%,  [a=0.1, b=2.0]);
eval(tst, y = 1.52847882672565/2): evalf(%);
eval(tst, y = 6.16747170396385/2): evalf(%);
                              -15                       -15
         [0.954929658551373 10   , -0.954929658551373 10   ]


                              -14                       -14
         [0.318309886183791 10   , -0.318309886183791 10   ]

I never liked that behaviour.

For the very question for a work around I would use

  y+piecewise(0<y, a,b);
  convert(%, piecewise, y); #lprint(%);
  codegen[makeproc](%, [y,a,b]);
  codegen[prep2trans](%);

Also the command y + `if`(0<y,a,b) would do and would compile, though I have not
checked for nested situations (=more inputs to piecewise).

A handmade way can results from

  y+piecewise(0<y, a,b);
  convert(%,pwlist,y);

though equality is not described in the help (but you have numerical cases in mind?)

Besides doing it manually it may be worth to look at Heaviside:

  y+piecewise(0<y, a,1 <y, b+a, b);
  convert(%, Heaviside);
  codegen[makeproc](%, [y::numeric,a::numeric,b::numeric]);
  codegen[prep2trans](%);
  p:=Compiler:-Compile(%);

  p(0.5, 2.1, 0.1);
                         2.60000000000000008

Have not looked at 'optimized' variants.

My practical problem once was to generate a (proper) code for Maple's
complex sign, ?csgn, without using complex types (through C++)
 

Your code/sheet shows the first error on "evalhf(MainF(u,v,w));"

It says "Error, lexical scoping is not supported in evalhf".

Entering MainF(u,v,w) one gets the answer 15.6660108548168786.

So it is not a function as expected, it is a constant.

Actually asking for u,v,w the answer is 5.9, 1.0, 6.4.

Entering MainF(_x, _y, _z) one gets "Error, invalid input: MainF expects
its 1st argument, u, to be of type float, but received _x"

That says: you may want to test before asking.

Actually I think, that this very experimental approach to learn
about a CAS like Maple does not make sense and can not be
the right way.

As already said: first ensure that the basic Maple code is correct,
then try to use evalhf. And only finally try to compile or even to
optimize or use complicated constructions.

And if posting: reduce to the very point, do not leave it to the
readers to find that out.


enc:=proc(expr)
# risky: assumes that y,a,b are THE indets of expr ...
local p;
p:= codegen[makeproc](expr, [y::numeric,a::numeric,b::numeric]);
#p:= codegen[prep2trans](p);
Compiler:-Compile(p, [y,a,b]);
end proc;
H:= ((as in your example above))
t:=enc(H);
t(1.0,2,3);
                         4.28730926684034230

@icegood , ok, I have have not recognized it, sorry.

I think the problem with your code is, how Maple understands
- the parameters a and b within the body
- the usage of the modifier, as james 1482 points out

I guess the compiler treats a and b as symbols (=not declared
variables in a C code)

Step-by-Step here would mean

enc:= (...)
local p
p:= unapply(...)
# Compiler:-Compile # deactivate for debuging
end proc;

Now t:= enc(H), and H something simple, say H:=y+a+b.

What do you expect for t(Y,A,B) and what do you actually get?
And why?

Only if you know that it makes sense to compile.

do not have Maple at my hands, now

But the compiler can not handle integrals and as you wrote it
in lower case Maple would come up with incomplete gamma,
I think.

However that may be a function, that the compiler does not know

Anyway I just would try to give the function expression directly
and not by the steps you use, so I can see what happens (and
what are the domains for my inputs)

Moshier's Cephes is a classical numerical library, here a sketch for computing
sin(x), http://www.netlib.org/cephes/doubldoc.html#sin

I do not comment much on range reduction (but integer multiples of Pi will need
increased precision for large values), but sin( Pi/2 - xi ) = cos(xi) allows to
reduce to the interval 0 .. Pi/4, if caring for both trigonometric functions.

Now to handle the Taylor series of sin(x) one can read it as x + x^3*series(x^2)
and E := t -> (sin(t^(1/2))-t^(1/2))/t^(3/2) is that analytic (!) function, that
means that sin(x) = x + x^3*E(x^2). Do not divide by x, but substract it.

Similar for cos(x) = 1 - x^2*someSeries(x^2).

Using Digits = 18 it turns out, that degree = 6 is enough for an approximation,
which reduces computational costs and possible errors through higher degrees:

  M:=numapprox:-chebyshev(E(t), t=0..evalf(Pi/4), 10^(-Digits+0)):
  M:=convert(eval(M, T=orthopoly[T]),horner):
  M:=unapply(evalf[18](M),t);

  Sin:= x -> x + x^3*M(x^2);

For visualization use acer's code (I would add style=point).

For higher precision find out that -Sum((-t)^(-1+k)/(2*k+1)!,k = 1 .. infinity)
is a series presentation for E(t), which is -1/6*hypergeom([1],[2, 5/2],-1/4*t).

As alternating series with decreasing coefficients it allows to estimate an
cut off error after n term.

Since one actually wants x^3*series(x^2) it turns out that n=8 is enough for
usual double precision, beyond one can use Lambert's W function to estimate
it - or just stops in a brute way, if no changes can be found. Or best rely
on Maple to do it.



Edited (after acer's comment): the code for the plot, g = as f,
but using Sin instead of S

xRange:= 0..Pi/4;
yRange:= -17..-15.5;
plot([g,F],xRange,numpoints=201,
   view=[xRange,yRange],color=[pink,blue],
   style=point, symbol=[circle, cross],
   myPlotDefault):
plot(-16,xRange, color=grey):
plots:-display(%%,%);

The relative error is ~ 1 DBL_EPSILON, may be better to use chebpade or minmax.
combine(GAL, radical) assuming x::real;

Almost always I use 15 Digits and often use rationals, and in most
cases I plot the situation:

restart;
Digits:=15;
-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772;
E:=convert(%, rational); #plot(E, x=-6 .. 3);
plot(E, x=-1 .. 2);

This suggestes, there may be 3 real solutions.

RootOf(E,x); [allvalues(%)]; evalf(%);

tells me: there are 2 Reals and 1 complex at x~1/4

eps:=1e-6;
plot(E, x=1/4 - eps .. 1/4 + eps);

shows me, that x=1/4 is not a zero, it is ~ 2e-9.

But you must be aware, that your input is only 10 decimals.
So it may be that is just by the input being not exact enough:

op(E); eval(%, x=1/4); evalf(%); evalf[10](%); 

      -0.6778764205, -7.738746974, -0.005859375000, 8.422482772

In that constellation it is not unlikely that just a rounding error
may be the reason, that one does not see a zero, the result is
by adding the terms.

Edited: be aware, that your factor 7.65 at the log is only 3 decimals.
7.649999997 resp 7.65000000213486 should give a zero in x=1/4
for Digits = 10 resp Digits = 15. Try to check that first.

First 36 37 38 39 40 41 42 Last Page 38 of 93