#1418 Problems with $specint

closed
5
2012-12-04
2008-05-16
No

I would like to open a bug report to collect the results of investigations of the function $specint. To get an overview of the problems with $specint I have taken the 99 tabulated Laplace transforms from the website of EqWorld.

55 of the 99 examples fail with the original code. I have divided the problems in the following cases:

1. Maxima has no algorithm:

In most cases Maxima gives an internal symbol like other-defint-to-follow-negtest or arbpow-failed. This is a known problem. In some cases we get a correct noun form of the unevaluated integral. At last there are many problems which gives a wrong answer.

A lot of examples include terms like t^(-1) ... or t^(-1/2) ... Maxima can't calculat these integrals but we know solutions.

A simple type of integral Maxima can't evaluate is the division by the sum of constants:

(%i7) specint(%e^(-s*t)/(x+y),t);
(%o7) other-defint-to-follow-negtest

In this cases the exponential function is hidden in a summation. I have found a correction which works generally and gives the correct result:

(%i7) specint(%e^(-s*t)/(x+y),t);
(%o7) (1/(x+y)*s)

2. Maxima has an algorithm for a special function but don't give the correct result:

We get no results for functions like bessel_k, bessel_y, log, erf, erfc etc. For all these functions Laplace transforms are tabulated. Beside the test of EqWorld I tried to get results for the internal functions %l[n,a](x) - the Laguerre function - or %he[n](x) - the Hermite function. But I dont' get the expected result.

Here the example for the Laguerre function:

(%i6) kill(all);
(%o0) done

(%i1) assume(s>0,n>0),declare(n,integer);
(%o1) [s > 0,n > 0]

(%i2) specint(%e^(-2*t)*%l[n,0](t),t);
Maxima encountered a Lisp error:
Error in MACSYMA-TOP-LEVEL [or a callee]: $N is not of type NUMBER.
Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

After correction of the code (the correction is not included in the diff appended to the next post):

(%i6) kill(all);
(%o0) done
(%i1) assume(s>0,n>0),declare(n,integer);
(%o1) [s > 0,n > 0]
(%i2) specint(%e^(-s*t)*%l[n,0](t),t);
(%o2) (1-1/s)^n/s

In the case of the Laguerre function I have found a bug in the transformation. After correction, Maxima gives the expected result shown above for the Laguerre function. There may be further bugs. Or limitations of the algorithmen prevend the calculation of results. At least, these limitations should be documented for the user.

3. Maxima gets extra factors or terms in the result

A simple example is the bessel_i function. Here we get an additional phase factor:

%e^(%i*%pi*v/(-1)^(v2/2)

in all calculations. This factor vanish when we introduce a small correction to the code. In other cases the problem seems to be more difficult.

Dieter Kaiser

Discussion

  • Dieter Kaiser

    Dieter Kaiser - 2008-05-16

    99 Tests for $specint

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-16

    diff beetween hypgeo-changed.lisp and hypgeo.lisp

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-16

    Logged In: YES
    user_id=2039760
    Originator: YES

    As a first step to improve the code of $specint I would like to present 7 changes:

    1. Function DEFEXEC

    If we cant't find a parameter, we apply $factor to the expression. Now Maxima finds the result for expressions like t/(x+y) --> 1/(s^2(x+y)) where x and y are free of the integration variable.

    2. Function ARBPOW1

    I have specialized the pattern match to be sure that in the expression c*t^v the parameter c is free of the integration variable. This condition now will fail if we enter $specint with expression like u(t) or t^(1/2)*(a+t)^(-1) and with the changes below we get nice and correct noun forms.

    3. Function LT-SF-LOG

    Because we have specialized the pattern match we add at the end of the function as return value a noun form.

    4. Function LT-ARBPOW

    A lot of integrals fail at this point. We add as the return value a noun form.

    5. Function LT-SF-LOG, Condition ONEI

    This is an example how we can avoid additional phase factors. If we use %i directly in the calculation all additional phase factors in the calculations vanish and the results are correct. There are more places we can apply this change to obtain easier results.

    6. Bug in LT-EXP and F35P147

    I have found a bug in the routine lt-exp and f35p147. This bug prevents the calculation of integrals with e.g. sin(2*sqrt(a*t)). $SPECINT gives the result 0.
    Here the output of Maxima after correction:

    (%i6) radcan(specint(%e^(-s*t)*sin(2*sqrt(a*t)),t));
    (%o6) sqrt(%pi)*sqrt(a)*%e^-(a/s)/s^(3/2)

    That is perfectly the tabulated expression and $SPECINT now works for a lot of other integrals too.

    7. Extension of the algorithm

    To show how we can extend the algorithm of $SPECINT to calculate further integrals, I have added code to calculate integrals of the form t^-1*(%e^(-a*t)-%e^(-b*t)). The code works also for integrals like t^-1*sin(a*t). Here an example

    (%i4) specint(%e^(-s*t)*t^-1*sin(a*t),t);
    (%o4) %i*(log(s-%i*a)-log(s+%i*a))/2

    That's equivalent to the tabulated answer atan(a/s).

    With this changes problems 55, 150 and 157 of rtest14.mac will produce different results. In all cases the noun form is improved and now more correct.

    The numbers of correct results of the test file test_eqworld.mac is increased to 87. When Maxima can't evaluate the integral but returns a correct noun form I declared the test as "(OK noun form)". There are 12 remaining examples which fails. This examples mostly include the log or erf function. I think there is something wrong with the mathematic.

    I have added a diff to show the above described changes to the code.

    Hint:

    The test file will stop 5 times and ask for the sign of the internal variable psey. That's a known bug. I have this bug not remarked as an error, because the results are correct. I try to find the reason of the bug.

    Dieter Kaiser

    File Added: diff_hypgeo.txt

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-24

    Logged In: YES
    user_id=2039760
    Originator: YES

    I have collected examples to test only the routine LT-EXP and the code for doing the Laplace transforms. I have found 30 examples tabulated by EqWorld and A&S testing this routine.

    With the original code 17 examples fail. The first changes showed in this thread reduce this to 6 examples. Correcting a bug in LT-EXP all examples work.

    8. Bug in LT-EXP and in the following routines F29P149, F36P147 and F37P147

    The reason for the 6 remaining errors is a bug in the calculation of the Laplace transform in the routines LT-EXP, F29P146, F36P147 and F37P147. The constant term is missing. Additionally it is useful to check v=0 before calling F36P147 and F37P147. For v <>0 we have no Laplace transform. Without this check we can get wrong results.

    BUT:

    For 3 examples I have an additional factor in the result of Maxima which is not present in tabulated result of EqWorld. In one case the factor is %e^(s^2/(s*a)) and in two cases the factor is 4. Furthermore, there is one case Maxima get a different sign of a term. Because the algorithm works very well in all other cases, I think there are errors in the table of EqWorld. But, it is better to verify thesw cases independly.

    Second:

    I have changed code in the routine LT-EXP and F35P147 to get the results for expression with sqrt(t) in the exponent (Point 6). The changed code works very well and gives a lot of additional results involving trigonometric and hyperbolic functions with sqrt(t) as argument. But, I haven't verified why this does work mathematically.

    I have attached a file with an update of the changes.

    Dieter Kaiser

    File Added: diff_hypgeo.txt

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-24

    Logged In: YES
    user_id=2039760
    Originator: YES

    I have attached the test file for the routine LT-EXP.
    File Added: test_hypgeo_LT-EXP.mac

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-26

    Logged In: YES
    user_id=2039760
    Originator: YES

    $SPECINT can not evaluate Laplace transforms for expressions with Logarithmic functions. The reason is that $SPECINT try to integrate the Hypergeometric representation (z-1)*2F1(1,1;2;1-z) of the Logarithmic function. But this doesn't work, because 2F1(1,1;2;1-z) doesn't satisfy the condition for the Laplace transform of a Hypergeometric function pFq.

    9.

    I have extended $SPECINT with an additional routine LT-LOG which calculates the Laplace transform for expressions like c*t^(v-1)*log(a*t). Further extensions can be added to calculate integrals with log(1+a*x), log(a+x) and log(x)^2.

    The formula I have got uses the scaling law f(a*t) -> 1/a*F(s/a) for Laplace transformations. There is a difference between the scaled results I get with my formula and the results of the Maxima function LAPLACE. Are the results of LAPLACE correct? Is something wrong with the scaled formula?

    I have added an updated diff.

    Dieter Kaiser

    File Added: diff_hypgeo.txt

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-26

    Code to calculate Logarithmic functions added

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-31

    Updated diff

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-05-31

    Logged In: YES
    user_id=2039760
    Originator: YES

    I have done further work on $specint. The test file test_eqworld.mac with 99 examples now gives 98 results or a correct noun form. A few examples have to be verified mathematically. Some examples gives factors which doesn't agree with the tabulated results. I am searching for the reasons. These small discrepancies can be very useful to detect some further problems with the used algorithm.

    There remains one problem with the tests from EqWorld. Maxima gets a result, but I can't verify it. This wrong result is difficult to understand. Perhaps we have an error in the routine hgfsimp-exec.

    I have collected the examples of A&S in a test file. We get about 130 integrals. With the orginal code 65 example pass the test. 66 example give a wrong result or not the correct noun form. With the changes up to now 91 examples of A&S pass the test. There are 40 examples remaining which gives an error. For most of the examples we have no algorithm, but Maxima don't give a correct noun form.

    To get the results up to now, I have done the modifications described in this thread and added the following one:

    1'. More general algorithm to handle a constant denominator

    Under point 1 I described a change to handle integrals of the typ "something/(a+b+ ...)", where the denominator is a constant. It is necessary to handle the case of a constant denominator more general. So I have added an algorithm to the function DEFINTEGRATE.

    10. Return noun form in $SPECINT

    It is not possible to construct the noun form easy and in all places of the code. It is more simple to test a flag in the routine $SPECINT and then if necessary to return the noun form. If added at some places code to return the flag 'hyp-return-noun-form. This method could be used at any places of the code.

    11. Routine A*X^M+C and EXECARGMATCH

    The pattern match for the argument of the hypergeometric function is too general. We match cases like a*x^2+b*x+c, but we have no algorithm for such an argument. In some cases Maxima gives a result which is wrong. After specialization of the pattern we get a Lisp error in the routine EXECARGMATCH. We have to return a list if the pattern match fails.

    The testsuite of Maxima has no problems with the changes. We only get three different noun forms (Problems 55, 150 and 157 in rtest14.mac).

    I have attached an updated diff.

    Dieter Kaiser

    File Added: diff_hypgeo.txt

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-06

    Updated diff after testing A&S

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-06

    Logged In: YES
    user_id=2039760
    Originator: YES

    I have finished a first look through about 130 examples of Laplace transforms which are tabulated by A&S. A lot of the integrals are allready tested. But sometimes even know integrals gives new problems, because A&S scale the integrals different from EqWorld.

    A lot of the integrals give a noun form. But as described earlier we get additional results too. Two problems remain open. Maxima gives a result for the functions sin(2*sqrt(k*t))/sqrt(%pi*t) and sinh(2*sqrt(k*t))/sqrt(%pi*t). But the results seems to be wrong. Maxima gets an additional factor erf(sqrt(k).

    To get the results I had to add further modifications to the code of $SPECINT:

    10.' Return noun form in $SPECINT

    To get more noun forms I have added a flag which will be tested in $SPECINT. A&S gives some examples which show that this method isn't general enough. Because we add up the symbols in the routine DISTRDEFEXEC we sometimes get 0 as a result. That is wrong.

    A more general and fast to implement solution is to set a global variable if we have to return a noun form. The disadvantage of this simple method is that Maxima sometimes get a solution for a part of a sum. This solution we loose with this method. Perhaps this problem has to be discussed more deeply.

    For a general solution we have to program a more complex algorithm.

    12. Bug in LTW

    A&S gives examples with the Laguerre and Hermite functions. In both cases the algorithm of $SPECINT don't work. For the Laguerre function it is possible to correct the code. In the routine LTW the Laguerre functions is transformed to a Whittaker W function. Correct is a Whittaker M function. In addition the correct Gamma factors have to be used. The Factorial function don't work generally.

    Up to now I haven't found the reason for the problems with the Hermite function.

    I have attached an updated diff. The testsuite has no problems with the changes.

    Dieter Kaiser
    File Added: diff_hypgeo.txt

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-06

    About 130 tests collected from A&S

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-06

    Logged In: YES
    user_id=2039760
    Originator: YES

    In addition I have attached the examples by A&S.

    File Added: test_A&S.mac

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-06

    Update of the tests collected from EqWorld

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-06

    Logged In: YES
    user_id=2039760
    Originator: YES

    Because of some changes in the code to return noun forms, I had to update the tests of EqWorld too.
    File Added: test_eqworld.mac

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-14

    Logged In: YES
    user_id=2039760
    Originator: YES

    I try to extend $SPECINT to allow Maxima to give more results. I implemented the following five cases for a sum of Exponentials functions:

    7'. Extension of the algorithm

    c * t^-1 * (%e^(-a*t)-%e^(-b*t))
    c * t^(-3/2) * (%e^(-a*t) - %e^(-b*t))
    c * t^-2 * (1 - 2 * %e^(-a*t) + %e^(2*a*t))
    c * t^-1 * (1 - 2 * %e^(-a*t) + %e^(2*a*t))
    c * t^-1 * (1 - %e^(2*a*t))

    The first code I used to implement the first of the above cases wasn't general enough and buggy (I think it worked accidentally). I had to redesign it a bit.

    All tabulated integrals by EqWorld and A&S with Trigonometric and Hyperbolic functions and a factor of t^-1 and t^-2 expand in a sum of the above cases. The sum with a factor t^(-3/2) I have added because A&S has a tabulated result.

    Examples with a factor t^-1 are:

    (%i7) specint(%e^(-s*t)*sin(a*t)/t,t);
    (%o7) %i*log((s-%i*a)/(s+%i*a))/2

    (%i8) specint(%e^(-s*t)*sinh(a*t)/t,t);
    (%o8) -log((s-a)/(s+a))/2

    (%i9) specint(%e^(-s*t)*(1-cos(a*t))/t,t);
    (%o9) log(a^2/s^2+1)/2

    (%i10) specint(%e^(-s*t)*(1-cosh(a*t))/t,t);
    (%o10) log(1-a^2/s^2)/2

    (%i11) specint(%e^(-s*t)*sin(a*t)^2/t,t);
    (%o11) log(4*a^2/s^2+1)/4

    This is an example with a factor t^-2:

    (%i12) specint(%e^(-s*t)*sin(a*t)^2/t^2,t);
    (%o12) -((s+2*%i*a)*log(s+2*%i*a)+(s-2*%i*a)*log(s-2*%i*a)-2*s*log(s))/4

    These examples include the integrals from the bug report [1477965] which now have correct results.

    I have attached a diff. Furthermore I had to update the test files from EqWorld and A&S because we get more results instead of noun forms. The testsuite has no problems with the extension of the code.

    Dieter Kaiser

    File Added: diff_hypgeo.txt

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-06-14

    Update - More results for Trigonometric and Hyperbolic functions

     
  • Raymond Toy

    Raymond Toy - 2008-06-20

    Logged In: YES
    user_id=28849
    Originator: NO

    I have now applied the hypgeo patch. However, as you mention, there are failures in rtest14 now.

    Can you correct that?

    Also, do you want test_eqworld, test_A&S, and lt-exp tests to be part of the testsuite?

     
  • Raymond Toy

    Raymond Toy - 2008-06-20
    • assigned_to: nobody --> rtoy
     
  • Dieter Kaiser

    Dieter Kaiser - 2008-07-18

    Logged In: YES
    user_id=2039760
    Originator: YES

    The suggested changes have been added to the code of $specint. Are any problems known? When not I would like to close the bug report.

    For further work on $specint I would prefer to open a new bug report.

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-07-21

    Logged In: YES
    user_id=2039760
    Originator: YES

    There seems to be no known problems with the changes up to now. For further work on $specint a new bug report can be opened. Closing this bug report.

     
  • Dieter Kaiser

    Dieter Kaiser - 2008-07-21
    • status: open --> closed
     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks