Double integral with singularity

Question

Jinsoo 2013-1-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/58490-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?

6 个评论
显示 4更早的评论隐藏 4更早的评论

Roger Stafford 2013-1-11

在 MATLAB Online 中打开

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 2013-1-11

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 2013-1-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/58490-double-integral-with-singularity#answer_70779

在 MATLAB Online 中打开

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.