# Double integral with singularity

3 views (last 30 days)
Jinsoo on 11 Jan 2013
Hi, I want to integrate the following equation (singularity at t = \tau)
1./sqrt(4*pi*(t-\tau)).*exp(-(x - \xi).^2./(4*(t-\tau)));
space x = (-10 <= x <= 10), time t = (0 <= t <= 2.1)
variable of integration \xi, \tau
\xi = x, 0 < \tau < t
So I tried as following code.
x = linspace(-10, 10, 41);
t = linspace(0.1, 2.1, 41);
xi = x;
ta = t;
nx = length(x); % number of interval of x
nt = length(t); % number of interval of t
for mt = 1:nt
for mx = 1:nx
GLA = @(ta, xi) 1./sqrt(4*pi*(t(mt)-ta)).*exp(-(x(mx) - xi).^2 ...
./(4*(t(mt)-ta)));
funGLA = quad2d(GLA, 0.05, ta(mt), xi(1), xi(end));
Int_GLA(mt, mx) = funGA;
end
end
However, it leads to wrong results. I don't know what's wrong. Could you help me?
Roger Stafford on 11 Jan 2013
First point. Why have you created these vectors:
xi = x;
ta = t;
You don't use them anywhere except as an integration limit which could be accomplished with different notation, and it seems to me there is a danger of confusing them as parameters to be entered in the function 'GLA' rather than its defined arguments (ta,xi). Also it may indicate a confusion on your part about how 'quad2d' is supposed to work. The arguments of GLA are not restricted to a set of fixed discrete values but are set by 'quad2d' in accordance with its own mysterious internal procedures.
Second point. The point where t(mt) approaches ta in GLA is not a true singularity in the mathematical sense since as ta approaches t(mt) the function approaches zero (L'Hospital's rule). However, it may cause computational difficulties for matlab which may see it as zero divided by zero and produce a NaN. It might be better to redefine GLA so as to make an exception at or very near t(mt) to give a zero there and avoid this difficulty.
You still haven't told us in what way the results are "wrong". This is important information in determining the cause of your trouble.
Teja Muppirala on 11 Jan 2013
Just checking, but you wrote two different variable names "funGLA" and "funGA".
funGLA = quad2d(GLA, 0.05, ta(mt), xi(1), xi(end));
Int_GLA(mt, mx) = funGA; <--------
Could that possibly be your problem?

Teja Muppirala on 11 Jan 2013
I will change my earlier comment into an answer:
Look at these two lines:
funGLA = quad2d(GLA, 0.05, ta(mt), xi(1), xi(end));
Int_GLA(mt, mx) = funGA;
You are calculating funGLA but then you are saying Int_GLA(mt,mx) = funGA
funGA is a variable remaining from an earlier calculation, and so every value of Int_GLA is getting assigned to that same number. You should first correct that mistake, and then see if that helps anything.