Numerical Integration [Matlab]

Prev: outline image object
Next: I need help with fin. timeseries and tomonthly

From: Michael on 28 Jul 2010 17:53

I am trying to perform numerical integration on a lot of data in order to come up with some calibration plots of sorts. I am needing to integrate Planck's black body emission multiplied by the camera's spectral response with respect to wavelength. Now Planck's equation is a function of wavelength and temperature and the spectral response is a function of wavelength. I have digitized the spectral response graph given to me by the camera manufacturer so my data points are unequal. I developed some code to perform this numerical integration on uneven intervals using two of Simpson's rules where applicable and the trapezoidal when necessary. I am taking a whole bunch of different temperatures at which to perform the integration.
When I run my code all I get is a linear set of answers. My logic tells me that when you integrate something non-linear you get something non-linear.
I'm new to these forums and Matlab doesn't come natural to me so I appreciate your patience.

Any suggestions of what could be going wrong?
Thanks!

From: John D'Errico on 28 Jul 2010 20:01

"Michael " <Michael.Newton(a)utah.edu> wrote in message <i2q8s1$ga3$1(a)fred.mathworks.com>...
> I am trying to perform numerical integration on a lot of data in order to come up with some calibration plots of sorts. I am needing to integrate Planck's black body emission multiplied by the camera's spectral response with respect to wavelength. Now Planck's equation is a function of wavelength and temperature and the spectral response is a function of wavelength. I have digitized the spectral response graph given to me by the camera manufacturer so my data points are unequal. I developed some code to perform this numerical integration on uneven intervals using two of Simpson's rules where applicable and the trapezoidal when necessary. I am taking a whole bunch of different temperatures at which to perform the integration.
> When I run my code all I get is a linear set of answers. My logic tells me that when you integrate something non-linear you get something non-linear.
> I'm new to these forums and Matlab doesn't come natural to me so I appreciate your patience.
>
> Any suggestions of what could be going wrong?
> Thanks!

The crystal ball is so cloudy. It just keeps on saying
"Probable user error."

Seriously, how are we expected to guess what you
have done wrong? Yes, you have likely made a mistake,
but you know that already, and without seeing anything
that you have done, no more can be said.

John

From: Michael on 29 Jul 2010 14:44

Okay, sorry for not posting my code yesterday. I was not on my computer with the code. I was simply wondering if there was something fundamentally wrong first. Here is my code as I have it now. I have checked the preliminary calculations up until performing the uneven numerical integration and everything seems to be okay.
Also, I have the wavelength and Quantum Efficiency data saved in a mat file as I had to digitize that seperately.

%%%%%%%%%%%%%%%%%%%%

h=6.6260693*10^-34; %Planck's Constant
e=1.60217653*10^-19; %elementary charge
c=2.99792458*10^8; %speed of light in a vacuum
kb=1.3806505*10^-23; %Boltzmann's constant
wavr=WavR; %wavelength data for red
n=length(wavr); %Define number of data points
Tred=linspace(500,10000,n); %create linearly spaced vector with Temp data (K)
CRred=e.*wavr.*QER./(h.*c.*100); %Calculate spectral response from Quantum Efficiency
IprimeR=zeros(n,n); %allocating size of matrix
for i=1:n %Loop to create matrix where each column represents diff. Temp.
for j=1:n
IprimeR(i,j)=(2*h*c^2/(wavr(i)^5))/(exp(h*c/(wavr(i)*kb*Tred(j)))-1);
end
end

IprimCRred=zeros(n,n);
for k=1:n %each column represents data for diff. Temps
for l=1:n
IprimCRred(l,k)=IprimeR(l,k).*CRred(l); %the inside of the integral
end
end

IntenseR=zeros(1,n);
for m=1:n
IntenseR(m)=uneven(n,wavr,IprimCRred(:,m));
end
%%%%%%%%%%%%%%%%%
Here are the functions uneven.m and its constituents
%%%%%%%%%%%%%%%%%
%Function for numerical integration with unevenly spaced data
%incorporates Trapezoidal, Simpson's 1/3, and Simpson's 3/8 rules where
%possible.
%Must have functions: Trap.m, Simp13.m, and Simp38.m

function Int=uneven(n,x,f) %n=number of data points,
%x=data that's being integrated with respect to, f=the inside integral data
h=x(2)-x(1);
k=1;
sum=0;
for j=1:n-1
hf=x(j+1)-x(j);
if abs(h-hf)<.0000001
if k==3
sum=sum+Simp13(hf,f(j-2),f(j-1),f(j));
k=k-1;
else
k=k+1;
end
else
if k==1
sum=sum+Trap(hf,f(j),f(j+1));
else
if k==2
sum=sum+Simp13(hf,f(j-1),f(j),f(j+1));
else
sum=sum+Simp38(hf,f(j-2),f(j-1),f(j),f(j+1));
end
k=1;
end
end
h=hf;
end
Int=sum;
end

%Function for numerical integration using Simpson's 1/3 rule
function sum=Simp13(h,f0,f1,f2)
sum=2.*h.*(f0+4.*f1+f2)./6;
end

%Function for numerical integration using Simpson's 3/8 rule
function sum=Simp38(h,f0,f1,f2,f3)
sum=3.*h.*(f0+3.*(f1+f2)+f3)./8;
end

%Function for numerical integration using trapezoidal rule (single segment)
function sum=Trap(h, f0,f1)
sum=h.*(f0+f1)./2;
end
%%%%%%%%%%%%%%%

Okay that should be it. Thanks for your help!

From: Michael on 29 Jul 2010 18:29

So, just out of curiosity I used MATLAB's built in trapz function in the same fashion as using my code above. This yielded very similar numbers to those I've already gotten; following the exact same trend when the two results were plotted together.
Either nothing is going wrong (highly doubt this) or it could be in the way I am passing the data to either of these functions. Or maybe I am screwing something up before this. I'm really not sure.

From: John D'Errico on 30 Jul 2010 09:28

"Michael " <Michael.Newton(a)utah.edu> wrote in message <i2svc0$4am$1(a)fred.mathworks.com>...
> So, just out of curiosity I used MATLAB's built in trapz function in the same fashion as using my code above. This yielded very similar numbers to those I've already gotten; following the exact same trend when the two results were plotted together.
> Either nothing is going wrong (highly doubt this) or it could be in the way I am passing the data to either of these functions. Or maybe I am screwing something up before this. I'm really not sure.

Since you have not posted your data (nor do I wish you
would do this on the newsgroup) I can't look at anything
to know what is the problem.

I would tend to guess, since you say that trapz does
similar things, I would guess it is a result of numerical
instability issues with a nasty looking function. But it
is no more than a wild guess at this point.

John

|
Pages: 1
Prev: outline image object
Next: I need help with fin. timeseries and tomonthly