use quadgk . the gauss kronrod method doesn't evaluate the function at the

boudaries. If the singularity is inside the interval, you must of cut the interval there

and integrate on the parts separately.

I agree. It depends on the type of improper integral. Basically you can use the tools at:

http://www.mathworks.com/help/symbolic/integration-1.html

adapted to your case. Careful, depending on the took you use, you can get crazy answers based on the machine limit, so make sure what you get makes sense. :-) Best of luck!

Hi Derya,

please provide a reference for the E_\alpha function in your pdf. Just that one is certain about the intended special function there. Did you analyze the nature of the singularities already? Especially those that might be caused by zeroes of the integral in (52)?

It might be also interesting to see the original problem to be solved, i.e. before the transformations.

Unfortunately, the .txt file mentioned in the .pdf is missing.

Besides the \Delta the integrand in (51) seems to be well-behaved in (0,\infty). Due to the oscillating J_0 term, you have to be careful here because of cancellation errors that might turn up using a Gauss rule of Laguerre type, say.

For the proper choice of your quadrature method, you should also analyze whether there are singularities in the complex plane close to the integration interval, and how these singularities are mapped when using coordinate transformations.

Best regards

Herbert

Thank you for your comments. However, I can not finf the correct way to calculate my function in Matlab (without using symbolic calculation). Prof. Herbert H H Homeier, I add the definition of E_alpha function (Mittag Leffler function which is generalization of exp function and is very often used in Fractional Calculus). As you see at the attachement. My functions can not be adapted to symbolic calculation in Matlab. If it is possible, please let you inform me about that. So I can not use the qudgk or any other symbolic commands.  There is alot of way to calculate the singular integral by numerically in Matlab toolbox. However, I dont write a symbolic algoritm. Your suggestions are very precious. 

Thank you for your help

At least for \bar b >0, tau0>0 and alpha>0 (the only case, I am going to discuss), I think, that there are no singularities on the path of integration. Maybe, that the two complex values +/- i (i=sqrt(-1)) pose some difficulties. To evaluate the inner integral, you need a good algorithm to evaluate the Mittag-Leffler function also for large negative z accurately, especially for small positive alpha. Just summing the defining power series is not sufficient here, as it is not sufficient for the related case of calculating exp(-z). For small zeta, the \bar\Delta integral is probably O(\zeta^2), thus, \zeta^2/\bar\Delta should approach a constant. Assuming that the Mittag-Leffler function approaches 0 for large negative z and positive alpha, the \bar\Delta integral should be rather well approximated as int( (1+v^2)^(-2) ,v=0..infinity)=Pi/4. But then, the integrand of the F-integral should be well behaved, and for \bar b >0,  the integrand becomes essentially an exponentially damped oscillation. Provided, that you can evaluate the inner integral accuratately enough also for large \zeta, you can try a Gauss-Laguerre rule for the outer integral. I do not know whether this is available in matlab, however. For the inner integral, you could try a coordinate transformation to a finite interval, followed by automatic integration.

Dear Prof Homeier, I am really thank you for your kindly help. I will analyze you comments as soon as possible and then I will give more detailed information about our calculations in Matlab. I am so sorry forget to mention the arbitrary values of  \bar b=1, \tau=1 and \alpha changes zero to one. In addition we calculate the mittag leffler function in \bar delta integral by using function [e]=mlf(alf,bet,c,fi)  in matlab. mlf comment is ready to use in matlab to calculate the mittag leffler funtion. We see after the calculation of delta integral after changing of variable (v=tanx), the integral is convergent. There is not any divergency at the results. Therefore we consider that the divergency occur in the outer integral. Because after change of variables as tany similar to the inner integral, we can not remove the singularity. The boundries change as 0 to pi/2.  And when we aware of the singulariy we need to consider: how will we calculate the integral, unfortunately without pass into symbolic matlab. If you need more inform about the problem I will be pleased to give more inform and also be happy to hear your suggestions. With best regards

Hi Derya,

I wonder whether _without_ your tan transformation this is a singular problem at all. If you would like to use a transformation to a finite interval, I suggest that you try a different one, at least for the outer integral.

An approach that is relatively easy to implement without the $tan$ transformation is the following: Transform the outer integral into a alternating series by integration between successive zeros of the J_0 function where in each of the finite intervals you use a Gauss-Legendre quadrature to obtain the terms of the series. Then, you can either try to evaluate the resulting series directly, or by using convergence acceleration that - in case of alternating, linearly convergent series (the exponential decay in your problem lets me expect that the series is linearly convergent, btw.)  is a stable process. As convergence acceleration method, you might try the epsilon algorithm, the Levin transformation or e.g. some variant of J transformation that is also a Levin-type transformation.

The latter (eg. the {}_pJ transformation) is rather easy to implement, but I do not know whether this has been done in matlab.

For this approach and Levin-type sequence transformations, see the review in the link below.

There has also been some work of Sidi and others on the combination of quadrature problems and extrapolation methods.

Best regards

Herbert

Article Scalar Levin-Type Sequence Transformations

my functin is 1/(1-x) from 0 to 1.Can anyone find out the integral???