Is there any easy way to calculate integral with cylinder coordinates in Matlab numerically?

11 views (last 30 days)
Dag Borlaug
Dag Borlaug on 19 Jan 2022
Edited: John D'Errico on 19 Jan 2022
Is there any easy way to calculate integral with cylinder coordinates in Matlab numerically? (3 dimensional space)
Can I use: integral3 for this work?
  1 Comment
Torsten
Torsten on 19 Jan 2022
If you include the transformation between carthesian and cylindrical coordinates on your own, you can use integral3, of course.
An easier way ? I think: no.

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 19 Jan 2022
Edited: John D'Errico on 19 Jan 2022
It seems pretty easy. Go back to calc 101. Maybe 102. WTP?
For example, what is the integral of z*(x^2+y^2), over a unit cylinder. So a cylinder of unit radius, centered at the origin, where lets say the cylinder has base at z==0, and top at z==1. Just making things up, so I have no idea what the solution is. We can trivially do it analytically. (Actually, pencil and paper should suffice too.)
syms z r theta
x = r*cos(theta);y = r*sin(theta);
Kernel = z*(x^2 + y^2);
int(int(int(r*Kernel,theta,[0,2*pi]),r,[0,1]),z,[0,1])
ans = 
I hope you remembered the differential element has volume r*dr*dtheta*dz. Consider that extra factor of r in there, multiplied by Kernel. That is important.
But can we do it numerically? Again, of course. Make sure you use the dotted operators, for exponentiation and multiplication. Your function needs to be properly vectorized.
x = @(r,theta) r.*cos(theta);
y = @(r,theta) r.*sin(theta);
Kernel = @(r,theta,z) z.*(x(r,theta).^2 + y(r,theta).^2);
integral3(@(r,theta,z) r.*Kernel(r,theta,z),0,1,0,2*pi,0,1)
ans = 0.7854
That should look very much like pi/4.
A simple Monte Carlo should also work. I'll just use rejection to do the sampling. And I'll sample in the first quadrant only, then multiply by 4.
n = 1e8;
x = rand(n,1);
y = rand(n,1);
z = rand(n,1);
keep = sqrt(x.^2 + y.^2) <= 1;
x = x(keep,1);
y = y(keep,1);
z = z(keep,1);
4*sum(z.*(x.^2 + y.^2))/n
ans = 0.7854
But is it easy? Easy is difficut to quantity, since easy for one person is not necessarily easy for another.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!