Double integral with singularity

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?

[Walter Roberson]

"it doesn't work" is quite vague. What happens? Does your computer start playing "A Bicycle Built For Two" ?

[Jinsoo]

Sorry, it is quite vague. So I edit it.

[Walter Roberson]

How can you tell that the results are wrong?

[Jinsoo]

I think it's wrong because the components of "Int_GLA" are all the same as 4.467159048042760e-06. It should be the form of something like wave.

Finally, I want to double integration of the following equation:

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

By using the built-in function. So could you please give me some advises.

[Roger Stafford]

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]

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?

Answer 1

[Teja Muppirala]

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.