integral function doesn't work well！

I have eidted my own function like this:

      function y = integral_Joe(a,b)
      if b<0 || b<a
          error('Error. Check the integral limit, they are not correct!');
      else
          fun1 = @(x)x.^3.*exp(x)./(exp(x)-1).^2;
          fun2 = @(x)x.^3./(exp(x)-1) + x.^3./(exp(x)-1).^2;
          if a*b < 0
              y = integral(fun1,a,0) + integral(fun2,0,b);
          elseif a*b >= 0
              y = integral(fun2,a,b);
          end
      end

when I input in the command window:

    integral_Joe(0,10)

    ans =

        7.1503

the function works well,but when I give the following sentence:

    integral_Joe(0,1e8)

    ans =

      8.1049e-17

the answer is obviously wrong!

Answer 1

[David Goodmanson]

Hi Joe, this is interesting and probably worth a second question on this website for the following reason. Here is a function called f3a. How this was arrived at is in the boring details below. But:

fun3a = @(x) x.^3.*exp(-x)./(1-exp(-x));

integral(@(x) fun3a(x), 0,1e3)    ans  = 6.4939
integral(@(x) fun3a(x), 0,1e6)    ans  = 6.4939
integral(@(x) fun3a(x), 0,7e7)    ans  = 6.4939
integral(@(x) fun3a(x), 0,8e7)    ans = 1.8879e-12            
integral(@(x) fun3a(x), 0,2e8)    ans = 2.1804e-39
integral(@(x) fun3a(x), 0,4e8)    ans = 9.8628e-86
integral(@(x) fun3a(x), 0,6e8)    ans = 1.4116e-132
integral(@(x) fun3a(x), 0,8e8)    ans = 1.2613e-179

however

integral(@(x) fun3a(x), 0,inf)    ans  = 6.4939

That's the issue. Some smart individual could probably reverse engineer this and find out what's going on.

======== boring details ==========

Ignoring the selection logic in front you have several functions,

fun1 =  @(x) x.^3.*exp(x)./(exp(x)-1).^2;
fun3 =  @(x) x.^3./(exp(x)-1);
fun4 =  @(x) x.^3./(exp(x)-1).^2;

I split fun2 into its two halves, with fun2 = fun3 + fun4.

Algebraically, fun2 = fun1.

As far as numerical problems it turns out that fun3 and fun4 are similar, so for simplicity fun4 is dropped. Now in fun1 the exp(x) in both numerator and denominator cause a problem for large x because the function blows up before Matlab has a chance to divide them. Better to multiply numerator and denominator by decreasing exponentials and get

fun1a = @(x) x.^3.*exp(-x)./(1-exp(-x)).^2;

For comparison

integral(@(x) fun1(x),  0, 100)    ans = 7.2123
integral(@(x) fun1(x),  0,1000)    ans = NaN

integral(@(x) fun1a(x), 0, 100)    ans = 7.2123
integral(@(x) fun1a(x), 0,1000)    ans = 7.2123

and fun1 is not considered further.

fun3 does not have the same num-den problem with exponentials as does fun1 but putting it into the form

fun3a = @(x) x.^3.*exp(-x)./(1-exp(-x));

avoids all numerical questions about increasing exponentials.

As far as numerical problems it turns out that fun3a and fun1a are similar, so the simpler function fun3a is kept for evaluation.