Trouble evaluating integral where function is made of huge and tiny values

Hi there,

I'm trying to integrate a function which is the product of one tiny and one large value. Analytically, this makes the integrand work out but numerically I am running into problems. Say I have the following parameters:

ChiSq = @(x,y) a.*x.^2 + b.*x.*y + c ;
Chi = @(x,y) sqrt(a.*x.^2 + b.*x.*y + c );
n = 840; %degrees of freedom, a fixed number
a = 1e6;
b = -1.5e5; 
c = 837; %a,b,c are parameters fit from experimental data

I plug these functions into the main function I want to integrate:

Zfun = @(x,y) (Chi(x,y).^(n/2-1)).*exp(-1/2*ChiSq(x,y));

The problem is, I can't integrate Zfun because Chi(x,y)^n/2-1 returns INF and exp(-1/2ChiSq) returns 0. If they could all be evaluated together the integrand would be reasonable, but due to numerical limitations this breaks.

How can I get around this?

UPDATE: I've tried the following but vpaintegral doesn't work with polar coordinates as far as I can tell. I'm confused about if the algorithm would integrate over a circular area if I transform into cartesian coordinates.

n = sym(840);
Zfun = @(x,y) (Chi(x,y).^(n/2-1)).*exp(-sym(1/2)*ChiSq(x,y));
polarZfun = @(theta,r) Zfun(r.*cos(theta),r.*sin(theta)).*r;
I(i) = vpaintegral(vpaintegral(polarZfun, theta, [0 2*pi]), r, [0 big_r]);

[Walter Roberson]

Zfun should probably refer to ChiSquared instead of ChiSq

What are the bounds of integration ?

Could we get some sample a and b values ?

[Are Mjaavatten]

One way to attack Inf/Inf situations like this is to split the expression into parts that remain finite. In your case you could try something like this:

n = 840;
Chi = 315; % Arbirary, but result is neithe very big nor very small
Z = 2^(-n/2)*Chi^(n/2-1)/gamma(n/2);
fprintf('Original expression for Chi = %5.2f, n = %d: %10.4f\n',Chi,n,Z) 
% NaN, as expected

Z = 1/2;
for i = 2:n/2
  Z = Z*Chi/2/(i-1);
end; 
fprintf('Stepwise product    for Chi = %5.2f, n = %d: %10.4f\n',Chi,n,Z)

% Now verify the method for much smaller n:
n = 6;
Z = 2^(-n/2)*Chi^(n/2-1)/gamma(n/2);
fprintf('Original expression for Chi = %5.2f, n =   %d: %10.4f\n',Chi,n,Z)

Z = 1/2;
for i = 2:n/2
  Z = Z*Chi/2/(i-1);
end; 
fprintf('Stepwise product    for Chi = %5.2f, n =   %d: %10.4f\n',Chi,n,Z)

However, for n = 820, this function varies very rapidly, from 0 for Chi = 52 to Inf for Chi = 1700, so unless Chi is restricted to a fairly small interval, the integration of Zfun may be challenging. If you could make do with a smaller vaue for n, things may become easier.

[supernoob]

Walter, thanks for pointing out the typo, I fixed that, as well as put some example a,b values in. I have also simplified the integrand as there are some numerical factors I can account for later. As for integration limits: I transform the function into polar coordinates then integrate out to r=0.3 over all angles (I integrate over a circle in the polar plane). Is it possible to use sym() and vpaintegral to do this? I couldn't find any information on vpaintegral in polar coordinates.

[supernoob]

Are, thanks so much for your suggestions. Unfortunately I can't change the value of n as it is a real parameter in the experimental data. You comment gave me some food for thought, though. It's appreciated.