How to calculate the integral of a function with a spline in it

Hello Everybody,

So i have a function called FB, and it is the product of a couple of functions

nB = spline(Wavelength,B3);
emi = @(x) 1.989e-6*(10^9*x).^2-0.002;
FB = @(x,T) nB(x).*emi(x)./(x.^5.*(exp((h.*c)./(x.*k.*T))-1))
%integration:
RaBG(Counts) = integral(@(x)FB(x,T),400e-9,720e-9)

nB is a function i must adquire from data. When fitting it with a polinomial function i get a small error which i believe is making me get wrong results, though the code works. I'm trying to fit it with a spline but haven't been able to get the code to work. I also tryied calling FB as this, which didn't work:

FB = @(x,T) ppval(nB,x).*emi(x)./(x.^5.*(exp((h.*c)./(x.*k.*T))-1));

Can anyone please help me?

[Andrew Newell]

If your data are evenly spaced, you could skip the spline fit and integrate using one of the Newton-Cotes formulae.

[Renan Kops]

Some of my data is evenly spaced so I can use this, good idea!

Any alternatives for the data that isn't evenly spaced?

[John D'Errico]

I'd like to chime in here,but without knowing what wavelength, B3, and G3 are, it is impossible to give a useful answer.

Answer 1

[Andrew Newell]

For evenly or unevenly spaced data, you could use the trapezoidal rule (MATLAB function trapz).

Interpolation is also reasonable. How exactly isn't the code working?

I don't see any values assigned to c, k or T. Are you doing that earlier? If so, it would make sense to define

FB = @(x) ppval(nB,x).*emi(x)./(x.^5.* (exp((h.*c)./(x.*k.*T))-1));

[Renan Kops]

Yes i am assigning values to all the other variables earlier, i didn't copy the entire code. When i call FB with ppval the code actually runs, but the results are way off, completely wrong.

I ran the code with the Newton-Cotes (Simpsons) method and it worked well with the evenly spaced data, I will try to use the trapz function with the ones that are unevenly spaced, thank you!

[John D'Errico]

Andrew is correct hee. You just need to evaluate the spline correctly. spline does not return a function you can evaluate directly as NB(x). You need to use ppval as he did.

Note that using trapz will also increase the error.

[Andrew Newell]

John, I'm not so sure that trapz will increase the error. Will an accurate integral of the spline be an accurate integral of the unknown function? Perhaps someone has analyzed this in the literature.

[Renan Kops]

Actually a correction, the value of T is assigned after the function is called. As i want to evaluate the integral for different values of T, i vary it using a 'while' loop.

I think it will help if I paste the rest of the code:

   emi = @(x) 1.989e-6*(10^9*x).^2-0.002;
   h = 6.6260693e-34;
   c = 2.99792e+08;
   k = 1.3806488e-23;
   nB = spline(Wavelength,B3);
   nG = spline(Wavelength,G3); 
   T = 500;
   Tp = ones(31,1)*500;  
   FB = @(x,T) ppval(nB,x).*emi(x)./(x.^5.*(exp((h.*c)./(x.*k.*T))-1));
   FG = @(x,T) ppval(nG,x).*emi(x)./(x.^5.*(exp((h.*c)./(x.*k.*T))-1));
   Counts = 1;
      while T<=2500
       RaBG(Counts) = integral(@(x)FB(x,T),400e-9,720e-9)/integral(@(x)FG(x,T),400e-9,720e-9);
       Tp(Counts) = T;
       T = T+50;
       Counts = Counts+1;
    end

So, as i mentioned, when calling FB like this the code works but the value returned is absurd, way different than what it should be.

[Andrew Newell]

Without the data, I can only guess. However, I notice that in the interval you're integrating over, you're dividing over the term

(x.^5.*(exp((h.*c)./(x.*k.*T))-1))

which is very small, so maybe you're dealing with numerical errors. Maybe you could rescale it so x is of order 1.

[Renan Kops]

Well it is possible, i will try to rescale it.

But, if that was the case, shouldn't i get wrong answers when integrating a 10th degree polinomial or with the Newton-Cotes method? Because both these methods worked and gave me a nice approximation of the real answer (For now i am trying to replicate the results of a paper to check if my code is performing as expected).

Thank you for all the help Andrew.

[Andrew Newell]

I'm glad to know that the alternative methods are working!