How to Integrate TWO times separately over different unknown parameters (THREE integrals in total) and come up with numerical result? My simulation is endless.

I will try to explain my problem as clear as possible. Thank you for your help in advance!

There are three integrals in my equation, w.r.t parameters "v","d" and "phi", as shown in the picture below. Their limits are [0 Inf], [0,30] and [0 pi] respectively.

NOTE that the latter Two of them are in the index (power) of "e", so they should be calculated first and then go to the former integral (the first one). The integrand of the latter two integrals have the parameter of the first integral "v" so that the result of these two integrals should hold this parameter and send it to the former integral.

is the Bessel Function with first kind and zero order, which can be written as "besselj(0,z)" in MATLAB.

I need numerical result at last. I tried to use function "int()" and "vpaintegral()" but both come up with symbolic results and cannot be converted, because of the multiple unknown parameters. Thus, using functions "integral()" and "integral2()" is the best way.

My idea is, first step calculating the latter two integrals using "@(v) integral2(@(d,phi) func, 0,30,0,pi)". Second step using "integral(@(v) ...)" and call the functional result of the first step.

step1 = @(v)integral2(@(d,phi) 0.0054.*besselj(0,v.*phi.*(30-d)).*phi./(8.*pi.*(1.854-1.848.*cos(phi)).^(3./2)),0,30,0,pi);
step2 = integral(@(v) 3.*exp(-2.5.*10.^(-7).*v.^2-4.53).*besselj(0,0.3.*v).*v.*(exp(step1(v))-1),0,Inf,'ArrayValued',true);

This can also be done by writing a single-line code.

result = integral(@(v) 3.*exp(-2.5.*10.^(-7).*v.^2-4.53).*besselj(0,0.3.*v).*v.*(exp(integral2(@(d,phi) 0.0054.*besselj(0,v.*phi.*(30-d)).*phi./(8.*pi.*(1.854-1.848.*cos(phi)).^(3./2)),0,30,0,pi))-1),0,Inf,'ArrayValued',true);

When running the code, there is a warning showing up repeatly, and the compiling is ENDLESS (keeps running). I guess the issue is in the second step. But if I change the coefficients in this equaion to a simpler form, sometimes the result comes up.

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

I am wondering if anything is wrong or improper here.

Answer 1

[David Hill]

I would try something like:

syms phi d v;
A=vpaintegral(3*exp(-2.5e-7*v^2-4.53)*besselj(0,.3*v)*v*...
    (exp(vpaintegral(vpaintegral((.0054*besselj(0,v*phi*(30-d))*phi)/...
    (8*pi*(1.854-1.848*cos(phi))^1.5),d,0,pi), phi,0,30))-1),v,0,inf);

I am not sure how much run time is needed. I ran it for several minutes with no result.

[Jianming Wang]

Thank you for your reply! I have tried using three "vpaintegral", or three "int" and even mixed style. The result is either endless or a symbolic equation that cannot be converted. I guess these complicated integrals are what MATLAB cannot perfectly solve.

[Paul]

The limits on the integration appear to be backward. Should be d -> 0-30, phi -> 0-pi.

Have you tried gridding the 3D space of v,d, and phi and then using trapezoidal integration, or something like that? How accurate does the result need to be?

[Jianming Wang]

@Paul Hi Paul, can you talk a bit more about the limits backward? Because I couldn't notice any issue there. I haven't tried trapezoidal integration, but I will look into it, or switch to another math calculator for this equation. High accuracy is NOT required, would you like to give any suggestion? Thx.

[Paul]

The original question states: "w.r.t parameters "v","d" and "phi", as shown in the picture below. Their limits are [0 Inf], [0,30] and [0 pi] respectively," which is consistent the code posted in the original question (but not the hand written equation, which shows dd*dphi, but probably should be dphi*dd based on the limits on the inner and outer integrals respectively). But the code posted in this answer shows the integration over d from 0 to pi, and over phi from 0 to 30.

[Jianming Wang]

@Paul Oh I see this issue in David's answer, thank you for reminding this to me and David as well. Fortunately I didn't make this mistake in my question though. I am looking into the trapezoidal integral to see if this method can solve this problem.

[Jianming Wang]

@Paul Hi, I tried to use trapezoidal method but I noticed that this method doesn't allow integration from "0" to "Inf" for parameter "v". I have to look into that to see if there is a solution to deal with the upper limit of "v".

[Paul]

Can't have v cover 0 to inf with trapezoidal, so instead need to integrate over v to some large number where it's believed the integral has converged close enough to the limit. I've been playing around with this problem and I'm seeing the integral over v oscillating. Presumably it's known that this integral converges in the limit as the upper bound on v approaches inf? If it does converge, it seems to be converging rather slowly from what I can tell.