How to prevent MATLAB from replacing a predefined integral with besselj?

I am asked via a partial differential equations class to verify that . This needs to be done by comparing the graph of besselj(0,x) over a pre-defined interval of x to the integral as stated, and see how close the graph of the integral is to MATLAB's default besselj(0,x) function.

However, MATLAB replaces the numeric computation of the integral with besselj(0,x) during the running of the code (I put the integrand in a separate function file and named it besselintegrand, then used the integral command in the file I will submit for the class), making the graphical comparison useless as it is computing besselj(0,x) twice and plotting both, resulting in the exact same graph for both the integral and besselj(0,x), instead of computing besselj(0,x) and then computing the integral and then graphing both functions as different entities. It is as if it is ignoring the integral, as something in MATLAB seems to go "Oh, this integral looks like besselj(0,x), let's just compute that instead of actually computing the integral."

Why do I want MATLAB to not default the integral to besselj(0,x)? Because I am asked to comment on the accuracy of the integral in comparison to MATLAB's default function. I can't get that if MATLAB replaces the integral with the default function that I am trying to compare the integral to!

Is there a way to prevent MATLAB from defaulting the computation of the integral to besselj, so that the actual output of the integral can be compared to besselj graphically? Or do I have to get sneaky and use some trigonometric identities to "trick" MATLAB into not recognizing the integral as besselj(0,x)? If there is a way to prevent MATLAB from doing this in the first place?

A classmate has tried using the symbolic math toolbox, and the system literally rewrites his integral in the output to besselj(0,x) so even the symbolic math toolbox has this behavior.

[Walter Roberson]

I would expect this for symbolic toolbox when int() was used on a symbolic expression or on a function handle. If you are not using int() or vpaintegral() then it would not be possible for besselj to appear.

[Galen Baumgartner]

I see, thanks for explaining. So, this is happening because I am defining the integrand as a function, and that causes MATLAB to compare the integrand with its existing default functions.

So, would quad cause the same problem? Does quad accept "vectorized" arrays? Or is there some other method that I am unaware of? Or is there some other command that could accomplish the same thing?

[Walter Roberson]

No, that explanation does not work.

I do not see how you could have the results you indicate if you do not use int() or vpaintegral(). You would not get this result with a numeric integral with integral() or the older quadgk() or similar, which is the approach that the assignment is expecting you to take.

[Galen Baumgartner]

Thanks again for clarifying. I was not using the symbolic toolbox but I was using integral() on a function handle, and I think my classmate used the symbolic toolbox and int(), so that may be why we had this problem.

I will rewrite my code without the symbolic toolbox and not as a function handle and see if I can compute using quadgk, it may take a day or two before I can get to it.

[Walter Roberson]

integral() and function handle do not generate that result.

For any given x

F = @(phi) cos(x.*sin(phi))
integral(F, 0, pi) /pi

[Galen Baumgartner]

I just tried your code, this is what I did:

x_bess=(0:0.01:5);
J_Bessel=besselj(0,x_bess)
integrand=@(phi) cos(x_bess.*sin(phi))
J_integral=(1/pi).*integral(integrand,0,pi,'ArrayValued',true)

As it is an array of x values, I needed to use the matrix operators. If you look at the numeric output, J_Bessel has, element by element, exactly the same values as J_integral, indicating that MATLAB is defaulting the integral command to besselj.

EDIT:

Unless the integral is EXACTLY equal to besselj and I am overthinking this... but I talked with my professor, and he says that there is a difference between the integral and besselj, and that MATLAB is screwing up for some reason. Not only that, he showed how using MathCAD he could demonstrate the difference between the integral and the bessel function!

I am very stuck, I do not know how to use quadgk to integrate an array (if that is what it takes), and because the integration is done over phi and not x but there is an array of x-values, quadgk is saying the matrix dimensions must agree, but I am struggling to figure out how to get the dimensions to agree because quadgk cannot use 'ArrayValued'.

[Walter Roberson]

max(abs(J_Bessel - J_integral))
ans =

      4.44089209850063e-16

That is round-off error that all numeric integration is subject to. You can adjust the 'AbsTol' and 'RelTol' parameters of, but you will get the same result. The maximum difference, shown above, is 4*eps() of the relevant entry (index 125) of J_Bessel, which is not much difference.

Indeed, if you check further, it is the besselj() evaluation that is further off.

      [maxdiff, idx] = max(abs(J_Bessel - J_integral))
      v = vpa(besselj(0,sym(x_bess(idx),'d')))
      db = v-J_Bessel(idx)
      di = v-J_integral(idx)
      simplify(abs(db) > abs(di))

you will see a result of symtrue, indicating that the besselj(0, x_bess(idx)) calculation is further away from the true version calculated symbolically to (default) 32 decimal places, than the integral is away from that extended value.

[David Goodmanson]

Hi Galen,

your code appears to be doing exactly what it should.

x = (0:0.01:5);
J_Bessel = besselj(0,x)
integrand = @(phi) cos(x.*sin(phi))
J_integral = (1/pi).*integral(integrand,0,pi,'ArrayValued',true)
figure(1)
plot(x,J_integral,x,besselj(0,x))   % overlay
grid on

format compact
% show the difference, a la mathcad
JminusJ = J_integral - besselj(0,x);
JminusJ(1:10) 
max(JminusJ)

ans =
   1.0e-15 *
  Columns 1 through 9
         0    0.2220         0   -0.1110         0         0         0         0   -0.1110
  Column 10
         0
ans =
   4.4409e-16

There are some very small differences, showing that 'integral' is working reallly well and that the residuals are probably a lot smaller than you thought.

[Galen Baumgartner]

I understand, but the professor won't accept this!

I guess I have no choice but to do as the professor recommended and just use MathCAD. I have the feeling that he is expecting us to find a much greater difference between the integral and besselj than what is currently shown.

EDIT:

I just did this in MathCAD... it is the same result.

I will just give this to my professor and state the tiny error difference, and hope that is what he is looking for, even though I am very certain it isn't.

[Walter Roberson]

I understand, but the professor won't accept this!

If that is the case, then the professor is wrong. Numeric integration will never be perfect, should never be expected to be perfect, and that should be part of the lesson.

You can switch to symbolic work to show that the integral is accurate to less than one bit in floating point precision:

x = 0:0.01:5;
J_Bessel = besselj(0,sym(x,'d'));
syms phi
integrand(phi) = cos(x.*sin(phi));
J_integral = (1/pi).*vpaintegral(integrand,phi,0,pi);

subplot(2,2,1)
plot(x, J_Bessel, '-b*');
title('besselj()')
subplot(2,2,2)
plot(x, J_integral, '--k.');
title('integral')
subplot(2,2,3)
plot(x, J_Bessel, '-b*', x, J_integral, '--k.')   % overlay
title('both together')
subplot(2,2,4)
JminusJ = J_integral - J_Bessel;
plot(x, JminusJ)
title('difference')

format compact
% show the difference, a la mathcad

fprintf('with differences calculated in symbolic domain:\n');

JminusJ(1:10) 
[maxdiff, maxidx] = max(abs(JminusJ));
fprintf('Max symbolic difference is %.16g at index %d\n', maxdiff, maxidx)
eps_factor = JminusJ(maxidx) / eps(x(maxidx));
fprintf('Number of epsilon is %g\n', eps_factor);

fprintf('with differences calculated in floating point domain:\n');
JminusJ = double(J_integral) - double(J_Bessel);

JminusJ(1:10) 
[maxdiff, maxidx] = max(abs(JminusJ));
fprintf('Max floating point difference is %.16g at index %d\n', maxdiff, maxidx)
eps_factor = JminusJ(maxidx) / eps(x(maxidx));
fprintf('Number of epsilon is %g\n', eps_factor);

[Walter Roberson]

If he is expecting to see a larger numeric difference, then you could always use integral() with a worse 'reltol' parameter.

[Galen Baumgartner]

I tried this:

syms x_bess phi
J_Bessel=besselj(0,x_bess)
J_integral=(1/pi).*int(cos(x_bess.*sin(phi)),phi,0,pi)
J_err=abs(J_Bessel-J_integral)

Running this code shows that J_integral defaults to a multiple of Jo(x). This isn't what it is supposed to do.

As a result, I can conclude that the system is using the default function of Jo(x) to approximate the integral, which gives the very tight tolerance results seen above. Clearly, this is not what the teacher is looking for.

[Walter Roberson]

>> J_integral=(1/sym(pi)).*int(cos(x_bess.*sin(phi)),phi,0,sym(pi))
J_integral =
besselj(0, x_bess)

As a result, I can conclude that the system is using the default function of Jo(x) to approximate the integral

No, you cannot conclude that at all !! int() is from the symbolic toolbox and does an analytic conversion when it can. integral() is a purely numeric tool that never ever converts anything to the besselj() function !! integral() does purely numeric approximation. It just happens to be doing a good job of numeric approximation.

I just tested with much reduced reltol and abstol, and it still does a quite good job of numeric approximation.

[Galen Baumgartner]

Look at the output once you run the code. You will see in J_integral some coefficient multiple of pi times Jo(x), instead of the integral argument. Because the output shows this, and substituting numeric values for the variables and evaluating the result produces exactly the same J_integral vector as is found via numeric integration via the integral command, I hope you can understand my confusion and frustration, and why my professor believes that this is happening and will not accept this answer.

[Galen Baumgartner]

I will try another tactic: use Euler's identity and evaluate the integral via that. Maybe that will be enough of a difference? I will get back once I have tried it.

EDIT: Nope, same difference....

EDIT AGAIN: Actually, the output no longer defaults to multiples of Jo(x) using the symbolic toolbox! I think that might be justification enough for my professor to say that the integral is besselj... I hope. I can check with him on Monday.

[Walter Roberson]

The coefficient multiple you were seeing was due to the approximation of 1/π in floating point. When the floating point number gets converted to symbolic form to multiply by the expression, the symbolic toolbox does not recognize the 1/pi floating point approximation as being the same as the reciprocal of the symbolic π that it has calculated, so what you are seeing in the output is π/(approximate pi)

The output I showed where I used 1/sym(pi) instead of 1/pi tells the symbolic calculator to use the transcendental π instead of the floating-point approximation, and that exactly balances the transcendental π that is calculated with the symbolic integral.

[David Goodmanson]

Hi Galen,

If you don't use syms, then there is no chance of defaulting to anything. 'integral' happens to be good at minimizing numerical error, hence the very small values. I don't know what point your prof is trying to make, but if he really wants to see some errors you can do

phipoints = 10;
phi = linspace(0,pi,phipoints);
x = 0:.1:10;
[xmat,phimat] = meshgrid(x,phi);
integrand = (1/pi)*(cos(xmat.*sin(phimat)));
J_integral = trapz(phi,integrand);
figure(1)
J_bessel = besselj(0,x);
plot(x,J_bessel,x,J_integral,'o')
grid on
J_err = max(abs((J_integral-J_bessel)))
J_err =
   3.0488e-04

This method sets up an array with x across, phi down (the equivalent of 'arrayvalued') and integrates down each column in phi. It uses trapz, the fairly primitive trapezoidal method.

Increasing phipoints improves the accuracy. At phipoints = 20 you are already down to an error of 3e-16, which is comparable to 'integral' and a lot better than I anticipated.