Integral2 for circular plot not working

1 回表示 (過去 30 日間)
Ren An Ooi
Ren An Ooi 2020 年 8 月 18 日
コメント済み: Ren An Ooi 2020 年 8 月 19 日
Hi there, I am testing out if integral2 can be used with scatterInterpolant to calculate the surface integral for discrete values.
The function polarfun(t,r) gives me a circular contour plot (of radius 1), and using integral2 gives me the surface integral U of this plot
However, I want to find the surface integral of similar plots whereby I do not have the function. Therefore I used scatteredInterpolant to obtain a continuous function from discrete values, then use integral2 again to find the surface integral U1 (which should be the same as U). Unfortunately, I get a NaN answer instead.
Can someone advice on what's the problem and how do I go about resolving it? Thanks a lot :)
func1 = @(t,r) (1-r).^(1/9) - 0.5/pi.*r.*(1-r).^(1/4).*t.*sin(t);
orient = 90/180*pi;
polarfun = @(t,r) ((1/pi)*( (1+r.*cos(t+pi/2+orient))./(r.^2+1+2*r.*cos(t+pi/2+orient))...
+ (1-r.*cos(t+pi/2+orient))./(r.^2+1-2*r.*cos(t+pi/2+orient)) )) .* func1(t,r);
U = integral2(polarfun,0,2*pi,0,1);
[t,r] = meshgrid(0:0.01:2*pi,0:0.01:1);
[x,y,zi] = pol2cart(t,r,polarfun(t,r));
cont = scatteredInterpolant(x(:), y(:), zi(:), 'linear');
U1 = integral2(@(x,y) cont(x,y), -1,1,-1,1);

回答 (3 件)

Walter Roberson
Walter Roberson 2020 年 8 月 18 日
make sure you put in an extrapolation method, even if it is the value 0 (no value outside the defined area.) The default is nan.
  1 件のコメント
Ren An Ooi
Ren An Ooi 2020 年 8 月 18 日
Thanks for the advice, I have tried both 'linear' and 'nearest' as the Extrapolation Method but it still does not work :(
The scatteredInterpolant function gives me "Warning: Duplicate data points have been detected and removed - corresponding values have been averaged." regardless of whether there's an extrapolation method
Also, the integral2 function gives me "Warning: Non-finite result. The integration was unsuccessful. Singularity likely."
Do you know what is going on? Thanks.

サインインしてコメントする。


Bruno Luong
Bruno Luong 2020 年 8 月 18 日
編集済み: Bruno Luong 2020 年 8 月 18 日
Have you tried to evaluate
>> polarfun(0,1)
ans =
NaN
or
>> cont(1,0)
ans =
NaN
Fix those issue first.
  2 件のコメント
Bruno Luong
Bruno Luong 2020 年 8 月 18 日
This brute truncation will give back a finite U1 result
func1 = @(t,r) (1-r).^(1/9) - 0.5/pi.*r.*(1-r).^(1/4).*t.*sin(t);
orient = 90/180*pi;
polarfun = @(t,r) ((1/pi)*( (1+r.*cos(t+pi/2+orient))./(r.^2+1+2*r.*cos(t+pi/2+orient))...
+ (1-r.*cos(t+pi/2+orient))./(r.^2+1-2*r.*cos(t+pi/2+orient)) )) .* func1(t,r);
U = integral2(polarfun,0,2*pi,0,1);
[t,r] = meshgrid(0:0.01:2*pi,1e-7:0.01:1); % avoid singularity at 0
[x,y,zi] = pol2cart(t,r,polarfun(t,r));
cont = scatteredInterpolant(x(:), y(:), zi(:), 'linear', 'nearest');
U1 = integral2(@(x,y) cont(x,y), -1,1,-1,1)
Ren An Ooi
Ren An Ooi 2020 年 8 月 19 日
Hi Bruno, thanks a lot for this advice :-) I have managed to obtain a finite value after avoiding singularity at 0!

サインインしてコメントする。


John D'Errico
John D'Errico 2020 年 8 月 18 日
編集済み: John D'Errico 2020 年 8 月 18 日
The warning message from scatteredInterpolant tell you there are duplicate data points in your data. Should you be remotely surprised at that? Of course not. you started with a regular grid in r and theta, so 0 to 1 in r.
How many points does that have for r == 0? THINK! TRY IT! Every one of those points at r == 0 are all the same, changing values of theta will not matter, since they all lie at the origin. The origin has r==0.
Interpolation methods cannot handle multiple points at the same location. So it throws an warning message.
As for generating a NaN, take Bruno's advice when generating a function. Is it well-defined everywhere in the domain?
But using scatteredInterpolant and integral2 is easy enough.
>> xy = rand(1000,2)*2 - 1;
>> xy(sum(xy.^2,2) > 1,:) = [];
>> size(xy)
ans =
766 2
That seems about right for the correct remaining fraction of points. I would expect to have roughly (pi/4)*1000 points remaining.
PLOT EVERYTHING THOUGH. Make sure you have what you think you have.
plot(xy(:,1),xy(:,2),'.')
axis equal
Can we use integral2 here? First, let me see if I can compute the area of that circle. I can do it in cartesian coordinates easily enough.
Zfun1 = @(x,y) double((x.^2 + y.^2) <= 1);
integral2(Zfun1,-1,1,-1,1)
ans =
3.14205282442276
The circle does indeed have area pi. The result is not too accurate, but then the kernel for that integral does have a circular singularity at radius 1. Can I do that in polar coordinates while using scatteredInterpolant?
[th,r] = cart2pol(xy(:,1),xy(:,2));
% working in the domain [0,2*pi]. cart2pol returns [-pi,pi].
% I would use integration limits of [-pi,pi] otherwise.
% but I want to work on [0,2*pi]. mod does the shift here.
th = mod(th,2*pi);
zones = ones(size(r));
Zfun2 = scatteredInterpolant(th,r,zones,'linear')
Zfun2 =
scatteredInterpolant with properties:
Points: [766×2 double]
Values: [766×1 double]
Method: 'linear'
ExtrapolationMethod: 'linear'
Here, I'll use scatteredInterpolant to compute the integral in polar coordinates. But something people problably forget that when you integrate in polar coordinates, you need to use the polar differential element as r*dr*dtheta. So do NOT forget that extra factor of r in there.
integral2(@(th,r) r.*Zfun2(th,r),0,2*pi,0,1)
ans =
3.1416
Again, the area of a circle is pi.
How about the volume under a paraboloid of revolution?, so z=x^2+y^2. If we have r go to 1, then I would expect pi/2 as the volume.
>> z = sum(xy.^2,2); % z = x^2 + y^2, z = r.^2 would have worked too.
>> Zfun3 = scatteredInterpolant(th,r,z,'linear');
>> integral2(@(th,r) r.*Zfun3(th,r),0,2*pi,0,1)
ans =
1.57793640036976
>> syms R
>> vpa(int(R.^2 * R,0,1)*2*pi)
ans =
1.5708
Again as expected. We could use nastier functions where I don't know the expected value, but then I would need to check the result to convince you it worked. Be careful to use the correct differential element. To verify your results is important, certainly until you get to the point you can trust your work. And even then check what you get.
Since I'm not sure what exactly you want to do, this should get you started.
  3 件のコメント
Bruno Luong
Bruno Luong 2020 年 8 月 19 日
編集済み: Bruno Luong 2020 年 8 月 19 日
"U2 (integrating in polar coordinates) gives the same answer as U (about 1.81), which is great! However, U1 (integrating in cartesian coordinates) gives an answer slightly different from U (about 2.00). May I ask what causes this difference in value?"
Because U1 you integrate on the rectangle that is larger than the disk. So U1 cannot be compare with U and U2.
(The data outside the disk is extrapolated with "nearest" point).
I have spot this and as well as error of missing "r" in polar integration in your first post that you have fixed.
Ren An Ooi
Ren An Ooi 2020 年 8 月 19 日
Makes sense, thanks!

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeLogical についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by