Marko Riedel

Mr. Marko Riedel

390 Reputation

9 Badges

8 years, 355 days
B.Sc. Computer Science UBC 1994, M.Sc. Computer Science UofT 1996.

MaplePrimes Activity


These are replies submitted by Marko Riedel

@vv 

With all due respect do you think this has a chance of being implemented? It does succeed on a considerable variety of integrals and the code effort is minimal. (You do have to take care to use the correct branch of the logarithm throughout.) Also, you want to stay ahead of your competitors. The Maple quest continues.

@vv 

Greetings once again. The goal being to improve the Maple integrator as much as possible and make it the best it can be as deserved by a first class computer algebra system I would like to present another application of my recursive method from the MSE post. This time we look to evaluate

int((log(x))^n/(x^3-2*x+4), x=0..infinity);

The following Maple session shows differentiating with respect to a parameter not working and Maple giving up on a closed form. The rest speaks for itself. Please do observe that I present this material to encourage an implementation of the recursive method or equivalent for the Maple integration engine.

    |\^/|     Maple 18 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2014
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> int(x^a/(x^3+1), x=0..infinity);    
                 1/3 Pi csc(1/3 Pi a + 1/3 Pi)

> int(x^a/(x^3-2*x+4), x=0..infinity);
                     infinity
                    /               a
                   |               x
                   |          ------------ dx
                   |           3
                  /           x  - 2 x + 4
                    0

> restart; read `cl2.maple`;          
alpha_sum := proc(n)
local poles;
    poles := [[1 + I, 1/2*log(2) + 1/4*I*Pi],
        [1 - I, 1/2*log(2) + 7/4*I*Pi], [-2, log(2) + Pi*I]
        ];
    add(residue(1/(x^3 - 2*x + 4), x = p[1])*p[2]^n,
        p in poles)
end proc

Q := proc(n)
local res;
option remember;
    if n = 0 then return
        simplify(int(1/(x^3 - 2*x + 4), x = 0 .. infinity))
    end if;
    res := -alpha_sum(n + 1)/(n + 1) - add(
        binomial(n + 1, p)*(2*I*Pi)^(n - p)*Q(p),
        p = 0 .. n - 1)/(n + 1);
    simplify(res)
end proc

                            infinity
                           /                 n
                          |            log(x)
           VERIF := n ->  |          ------------ dx
                          |           3
                         /           x  - 2 x + 4
                           0

> Q(4);
memory used=30.4MB, alloc=40.3MB, time=0.66
  139       3   2   4897          4              4
- ---- ln(2)  Pi  - ----- ln(2) Pi  + 9/640 ln(2)  Pi
  1920              76800

                               5
        63       2   3   357 Pi     31       5
     + ---- ln(2)  Pi  + ------- - ---- ln(2)
       1280               10240    1600

> evalf(%);
                          6.866005720

> VERIF(4);
                     infinity
                    /                 4
                   |             ln(x)
                   |          ------------ dx
                   |           3
                  /           x  - 2 x + 4
                    0

> evalf(%);
                          6.866005710

> diff(int(x^a/(x^3-2*x+4), x=0..infinity), a$4);
memory used=65.3MB, alloc=78.3MB, time=1.49
memory used=129.8MB, alloc=86.3MB, time=2.67
                     infinity
                    /           a      4
                   |           x  ln(x)
                   |          ------------ dx
                   |           3
                  /           x  - 2 x + 4
                    0


user@host:~/complex-logint$ math
Mathematica 10.0 for Linux x86 (64-bit)
Copyright 1988-2014 Wolfram Research, Inc.

In[1]:= Integrate[Log[z]^4/(z^3-2*z+4), {z, 0, Infinity}]        

                5          4                 3       2
Out[1]= (5355 Pi  - 9794 Pi  Log[2] + 7560 Pi  Log[2]  - 
 
               2       3                 4              5
>      11120 Pi  Log[2]  + 2160 Pi Log[2]  - 2976 Log[2] ) / 
 
>    153600

In[2]:= Expand[Out[1]]                                           

              5          4               3       2
        357 Pi    4897 Pi  Log[2]   63 Pi  Log[2]
Out[2]= ------- - --------------- + -------------- - 
         10240         76800             1280
 
           2       3              4            5
     139 Pi  Log[2]    9 Pi Log[2]    31 Log[2]
>    --------------- + ------------ - ----------
          1920             640           1600

In[3]:= N[Out[1]]                                                

Out[3]= 6.86601

In[4]:=                                                          
user@host:~/complex-logint$ 

Best regards, Marko Riedel

 

cl2-maple.txt

@vv Thank you very much for this instructive reply. I tried it just now and it works extremely well. We may hope that in the future Maple might perhaps recognize this substitution on its own. Presumably the user consults Maple because they did not know how to solve the integral in the first place. Your formula makes for a nice and simple pattern that should be easy to implement in the integration engine.

@Carl Love To make a long story short this was a combination of late night computing, a detailed write-up and a very simple answer. Thank you ever so much for explaining anyway and for taking the trouble to do the optimization. It looks impressive and I can follow it even though readability has suffered somewhat. In fact I had run into this issue before where I received a warning about iterative appending at your site from another user. You can be sure I won't make this mistake again.

I would like to say just the same that I was surprised to learn about the n squared /quadratic complexity of iterative appending. This seems quite poor and one would think that it ought to have been fixed in several decades of Maple computing. Some programming languages like Perl produce a warning when inefficient constructs are used. Maple itself warns about using old-style lists in place of modern Arrays. I suggest something be done about this. It might be possible to track how sets are grown and warn about this trap. These warnings could be controlled by a global flag.

Lastly I would like to point out that we have an interesting research problem here. I presented two algorithms, enumeration and Burnside. The question is whether there is an algorithm to compute the number of orbits between these two extremes that has complexity of a different order (better) than enumeration. This is related to the open problem of building a constant space iterator over necklaces and bracelets.

Best regards, Marko Riedel

 

 

@Doug Meade The link that I have in my current profile works in Google Chrome but not in Firefox.

@Carl Love Good to see that we are working this out. Another related issue is seen in the following transcript, where we may observe a second command missing the pole / singularity and returning a Taylor series instead of a Laurent series.

user@host:~$ maple
    |\^/|     Maple 18 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2014
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> series(z^2/(z^4 + 2*z^2 + 2)^2, z=2^(1/4)*exp(3*Pi*I/8), 3);
     1/2              2
    2    exp(3/8 I Pi)
1/4 ------------------- +
              2
            %2

    /   3/4              3                 1/4              \
    |  2    exp(3/8 I Pi)  (%1 + 1)       2    exp(3/8 I Pi)|
    |- ---------------------------- + 1/2 ------------------|
    |                3                             2        |
    \              %2                            %2         /

                                 /                    /
          (1/4)                  | 1/2              2 |
    (z - 2      exp(3/8 I Pi)) + |2    exp(3/8 I Pi)  |
                                 |                    |
                                 \                    \

                       6  1/2                  4
        5 exp(3/8 I Pi)  2    + 9 exp(3/8 I Pi)  - 1
    1/2 --------------------------------------------
                              4
                            %2

        1/2              2         2\      1/2              2
       2    exp(3/8 I Pi)  (%1 + 1) |   2 2    exp(3/8 I Pi)  (%1 + 1)
     + -----------------------------| - ------------------------------
                      4             |                  3
                    %2              /                %2

              \
            1 |       (1/4)               2           (1/4)               3
     + 1/4 ---| (z - 2      exp(3/8 I Pi))  + O((z - 2      exp(3/8 I Pi)) )
             2|
           %2 /

       1/2              2
%1 := 2    exp(3/8 I Pi)

                   4
%2 := exp(3/8 I Pi)  + %1 + 1

> evalf(%);
Float(infinity) - Float(infinity) I + (Float(infinity) + Float(infinity) I)

    (z - 0.4550898607 - 1.098684113 I) + (Float(infinity) + Float(infinity) I)

                                      2
    (z - 0.4550898607 - 1.098684113 I)  +

                                        3
    O((z - 0.4550898607 - 1.098684113 I) )

> quit
memory used=2.4MB, alloc=8.3MB, time=0.10
user@host:~

Post Scriptum. I understand about the two different radical forms and how residue needs the right one to work properly. But surely it must be easy to verify that the form that the convert command produced represents a singularity. This would seem to be a matter of straightforward complex algebra.

@Carl Love Thanks for explaining. Converting a complex number from algebraic to polar and the reverse are sufficiently common that one would expect Maple to provide this feature, precisely because an algorithm can be of considerable assistance in the simplifications that are necessary and which can be quite involved, as you have discovered. Then there is the fact that Mathematica succeeded. (Consult the MSE link for the details of how the simplification is carried out.)

@Doug Meade I tried the following URL: http://www.example-test.com/~username and it was rejected by the validator. The URL http://www.example-test.com was accepted.

@Doug Meade I am entering it into the URL field on the profile page. I see that the tilde works for you. I also have a dash in the URL, could that be the problem? (Just tried it again and got the invalid URL message.) The name of the field is "Personal Website URL."

@Doug Meade 

Thanks ever so much. I had looked under profile -- I thought settings referred to something else. How can I find out why it won't accept my home page URL? Is it the tilde in the URL? Is there a document available somewhere that describes the format of URLs permitted at the site?

Thanks!

Marko Riedel

Well in the sense that this is an algorithm which is not the most modern, whichis Goulden-Jackson. What do you think would be a better title?

@Carl Love My apologies, I took your comment to mean that I was wrong to attempt this summation symbolically while what you mean appears to be that the sum is incorrect. Regards.

@Carl Love I was not looking for an integer or a number. I used the initial command precisely with the goal of getting a symbolic answer. As you can see on the last two commands I can sum the polynomial easily myself and without using Maple. My purpose was to discover what form Maple's symbolic answer would take. The three lines should evaluate to the same polynomial. That they don't means we definitely have a problem. If I had not been looking for a symbolic answer I would have used the "add" and "mul" commands. The fact that the value four shows up in the denominator of the rational terms suggests that a pattern matcher is malfunctioning somewhere. Regards.

@Carl Love Thanks for mentioning this, the potential misunderstanding didn't catch my attention as my intent was to make it clear what the provenance and origins of the package are and at the same time communicate the positive nature of the user experience with this software.

@Markiyan Hirnyk 

This is a problem I have sometimes that I do not read the help docs carefully enough; I apologize. Let me add with all due respect that if I knew beforehand that my value is a function of 1+sqrt(2) I would not need the identify command. This is the clinch here.

Marko Riedel

1 2 3 Page 1 of 3