How to explain a difference in solutions of integrals between MATLAB and MAPLE?

Dear community,

after solving an integral via MAPLE 2022 (see the solution in blue),

I tried to reproduce the calculation in MATLAB R2021b as a numerical and symbolic expression:

%% integral: numerical expression
            
            %empirical coefficents
            vrel = 0.000018;        
            alpha = 0.08;           
            beta = -6.75;           
            h =1;                   
            d =1;                   

            %function
            fun = @(theta) (alpha*beta*cos(theta).^4).*(pi*d-8*beta*vrel*cos(theta).^2).^(-2);
            
            %integral
            F = integral(fun,0,pi/2)

F = -0.0322

%% integral: symbolic expression

            %system parameters
            syms theta alpha beta d vrel

            %function
            fun = (alpha*beta*cos(theta).^4).*(pi*d-8*beta*vrel*cos(theta).^2).^(-2);

            %substiution of empirical coefficents
            fsub = subs(fun,{alpha, beta, d, vrel}, {0.0814,-6.75,1,0.000018});

            %integral
            f = int(fsub,theta,0,pi/2);

            %evaluation for complex number
            c = subs(f, 1+i);
            F = double(c)

F = -0.0328

However, the solution between MATLAB and MAPLE differs by 8 magnitudes. My colleauge and me double-checked for typos but couldn't figure out any. I verified the result from MAPLE via comparison to a simplfied solution for the integral.

I am curious to learn about my mistake(s) when adapting the integral into MATLAB!

[Jiri Hajek]

Hi, I'd suggest to plot the function you are trying to integrete in both software tools. That way, you'll be able to spot if e.g. the software treats the angles in correct units, or operations order is correctly interpreted.

[Tim Hammer]

Hi Jiri,

thanks for your suggestion! I checked it, and it seems that both softwares evaluate the function correctly:

MAPLE:

MATLAB:

            %empirical coefficents
            vrel = 0.000018;        
            alpha = 0.08;           
            beta = -6.75;           
            h =1;                   
            d =1;                   

            %function
            fun = @(theta) (alpha*beta*cos(theta).^4).*(pi*d-8*beta*vrel*cos(theta).^2).^(-2);
            close all
            figure
            fplot(@(theta) fun(theta),[0 1.6])

[Bjorn Gustavsson]

After following @Jiri Hajek's suggestion one gets a smoothly varying curve starting at approximately -0.054 at theta=0 that aproaches zero at pi/2. This curve is everywhere in this interval 0<theta<=1 under the line y = -0.05 + 0.05*x. Therefore the integral has to have a smaller value than the area of the triangle with corners at [0,0], [0,-0.05] and [1,0] which should be -0.05*1/2 or -0.025.

Answer 1

[Jiri Hajek]

Nicely done graphs! So, just by looking at the graphs, it is clear that:

• The integral of your function is a negative value, with lower limit on the value being the area of the graph, i.e. about 1,6x(-0,05)= - 0,08
• Comparing your original results from MATLAB with the above, it seems that both provide good estimates, although it remains to see which one is more precise
• The MAPLE result is a complete nonsense, 8 orders of magnitude off the target.You have to search for error in your MAPLE implementation, not in MATLAB.

[Tim Hammer]

I will try to figure out why my maple solution is incorrect - but this does not belong here.

Thanks Jiri!

Answer 2

[Torsten]

https://www.wolframalpha.com/input?i=integration&assumption=%7B%22C%22%2C+%22integration%22%7D+-%3E+%7B%22Calculator%22%7D&assumption=%7B%22F%22%2C+%22Integral%22%2C+%22rangestart%22%7D+-%3E%220%22&assumption=%7B%22F%22%2C+%22Integral%22%2C+%22integrand%22%7D+-%3E%220.08*%28-6.75%29*%28cos%28x%29%29%5E4%2F%28Pi-8*%28-6.75%29*0.000018*%28cos%28x%29%29%5E2%29%5E2%22&assumption=%7B%22F%22%2C+%22Integral%22%2C+%22rangeend%22%7D+-%3E%22Pi%2F2%22

[Tim Hammer]

Thank you Torsten!

Answer 3

[Alan Stevens]

This is what I get in Maple:

[Tim Hammer]

Thank you Alan!