Hi all,
I would like to evaluate double integral of some function which contain one singular point at the origin. I must admit that I am still a newbie when dealing with singular function. I tried to use integral2 command, but I dont know why the singularity (at 0,0) is not detected. I get the answer (q1). I tried to put the singularity on the edge (I split the area) and sum the contribution of both, I get the result (q2+q3) which has more or less the same with q1 but not exactly. Then I tried with symbolic variable, I get the answer (q4), again, quite close with q1 even not exactly. Then, I would like to see the behaviour of the function in the domain, then I tried to plot and discretize the domain. I sum the values on the discretized domain excluding the one at the origin (which is infinite). I thought this summation is equal with the result of the integration. Apparently, its not. This value obviously depends on my discretization size. If I use smaller discretization, the more big values around the singular point will be included in the summation.
So, my question is, is there any correct results from approach that I've done above? What kind of theory I should use to solve this problem? How to tell whether the result obtained by Matlab is correct or not. Now, I am like a blind man. I ask Matlab to solve this problem, I get some results, but I dont know how to tell whether this result is correct or not.
Attached below is my code. Any input/suggestion will highly appreciated.
Regards,
Fredo Ferdian
if true
fun = @(x,y) (1./((x.^2+y.^2).^(0.5)));
format long;
q1 = integral2(fun,-1,1,-1,1);
q2 = integral2(fun,-1,0,-1,1);
q3 = integral2(fun,0,1,-1,1);
isequal(q1,q2+q3)
syms x y;
fx = int((x^2+y^2)^(-1/2),x);
fx1 = subs(fx,x,1);
fx2 = subs(fx,x,-1);
fxa = fx1-fx2;
fy = int(fxa,y);
fy1 = subs(fy,y,1);
fy2 = subs(fy,y,-1);
Value = fy1-fy2;
q4 = eval(Value);
isequal(q4,q1)
discrt = (1-(-1))/100;
[X,Y] = meshgrid(-1:discrt:1);
Coord = [X(:) Y(:)];
Q5 = 1./((Coord(:,1).^2+Coord(:,2).^2).^0.5);
Q51 = sum(Q5(1:(round(length(Q5)/2))-1));
Q52 = sum(Q5((round(length(Q5)/2))+1:end));
q5 = Q51+Q52;
isequal(q5,q1)
Q5A = reshape(Q5,[length(X(:,1)),length(X(1,:))]);
surf(X,Y,Q5A);
colorbar;
end
