MATLAB Examples

Arc Length Calculations for a Geometry Function

Contents

This example shows how to create a cardioid geometry using four distinct techniques. The techniques are ways to parametrize your geometry using arc length calculations. The cardioid satisfies the equation $r = 2 \left( 1 + \cos \left( \Phi \right) \right)$.

ezpolar('2*(1+cos(Phi))')

This example shows four approaches to giving a parametrization of the cardioid as a function of arc length.

  • Use the pdearcl function with a polygonal approximation to the geometry. This approach is general, accurate enough, and computationally fast.
  • Use the integral and fzero functions to compute the arc length. This approach is more computationally costly, but can be accurate without you having to choose an arbitrary polygon.
  • Use a parametrization that is not proportional to arc length plus a constant. This approach is simplest, but can yield a distorted mesh that does not give the most accurate solution to your PDE problem.

Polygonal Approximation

The finite element method uses a triangular mesh to approximate the solution to a PDE numerically. So there is no loss in accuracy by taking a sufficiently fine polygonal approximation to the geometry. The pdearcl function maps between parametrization and arc length in a form well-suited to a geometry function. Here is a geometry function for the cardioid.

function [x,y] = cardioid1(bs,s)
% CARDIOID1 Geometry File defining the geometry of a cardioid.
if nargin == 0
  x = 4; % four segments in boundary
  return
end
if nargin == 1
  dl = [0    pi/2   pi       3*pi/2
        pi/2   pi     3*pi/2   2*pi
        1      1      1        1
        0      0      0        0];
  x = dl(:,bs);
  return
end
x = zeros(size(s));
y = zeros(size(s));
if numel(bs) == 1 % bs might need scalar expansion
  bs = bs*ones(size(s)); % expand bs
end
nth = 400; % fine polygon, 100 segments per quadrant
th = linspace(0,2*pi,nth); % parametrization
r = 2*(1 + cos(th));
xt = r.*cos(th); % Points for interpolation of arc lengths
yt = r.*sin(th);
% Compute parameters corresponding to the arc length values in s
th = pdearcl(th,[xt;yt],s,0,2*pi); % th contains the parameters
% Now compute x and y for the parameters th
r = 2*(1 + cos(th));
x(:) = r.*cos(th);
y(:) = r.*sin(th);
end

Plot the geometry function.

pdegplot('cardioid1','EdgeLabels','on')
axis equal

With 400 line segments, the geometry looks smooth.

The built-in cardg function gives a slightly different version of this technique.

Integral for Arc Length

You can write an integral for the arc length of a curve. If the parametrization is in terms of $x \left( u \right)$ and $y \left( u \right)$, then the arclength $s \left( t \right)$ is

$$s \left( t \right) = \int_{0}^{t} \sqrt{\left( \frac{dx}{du} \right)^2
+ \left( \frac{dy}{du} \right)^2} du .$$

So for a given value $s0$, you can find $t$ as the root of the equation $s \left( t \right) = s0$. The fzero function solves this type of nonlinear equation.

For the present example of a cardioid, here is the calculation.

function [x,y] = cardioid2(bs,s)
% CARDIOID2 Geometry file defining the geometry of a cardioid.
if nargin == 0
  x = 4; % four segments in boundary
  return
end
if nargin == 1
  dl = [0    pi/2   pi       3*pi/2
        pi/2   pi     3*pi/2   2*pi
        1      1      1        1
        0      0      0        0];
  x = dl(:,bs);
  return
end
x = zeros(size(s));
y = zeros(size(s));
if numel(bs) == 1 % bs might need scalar expansion
  bs = bs*ones(size(s)); % expand bs
end
cbs = find(bs < 3); % upper half of cardiod
fun = @(ss)integral(@(t)sqrt(4*(1 + cos(t)).^2 + 4*sin(t).^2),0,ss);
sscale  = fun(pi);
for ii = cbs(:)' % ensure a row vector
    myfun = @(rr)fun(rr)-s(ii)*sscale/pi;
    theta = fzero(myfun,[0,pi]);
    r = 2*(1 + cos(theta));
    x(ii) = r*cos(theta);
    y(ii) = r*sin(theta);
end
cbs = find(bs >= 3); % Lower half of cardioid
s(cbs) = 2*pi - s(cbs);
for ii = cbs(:)'
    theta = fzero(@(rr)fun(rr)-s(ii)*sscale/pi,[0,pi]);
    r = 2*(1 + cos(theta));
    x(ii) = r*cos(theta);
    y(ii) = -r*sin(theta);
end
end

Plot the geometry function.

pdegplot('cardioid2','EdgeLabels','on')
axis equal

The geometry looks identical to the polygonal approximation. This integral version takes much longer to calculate than the polygonal version.

Analytic Arc Length

If you are handy with integrals, or have Symbolic Math Toolbox, you can find an analytic expression for the arc length as a function of the parametrization. Then you can give the parametrization in terms of arc length. Here is an approach using Symbolic Math Toolbox.

syms t real
r = 2*(1+cos(t));
x = r*cos(t);
y = r*sin(t);
arcl = simplify(sqrt(diff(x)^2+diff(y)^2));
s = int(arcl,t,0,t,'IgnoreAnalyticConstraints',true)
 
s =
 
8*sin(t/2)
 

So you see that, in terms of arclength s, the parameter t is t = 2*asin(s/8) where s ranges from 0 to 8, corresponding to t ranging from 0 to $\pi$. For s between 8 and 16, by symmetry if the cardioid, t = pi + 2*asin((16-s)/8). Furthermore, you can express x and y in terms of s by the following analytic calculations.

syms s real
th = 2*asin(s/8);
r = 2*(1 + cos(th));
r = expand(r)
 
r =
 
4 - s^2/16
 
x = r*cos(th);
x = simplify(expand(x))
 
x =
 
s^4/512 - (3*s^2)/16 + 4
 
y = r*sin(th);
y = simplify(expand(y))
 
y =
 
(s*(64 - s^2)^(3/2))/512
 

Now that you have analytic expressions for x and y in terms of the arclength s, you can write the geometry function.

function [x,y] = cardioid3(bs,s)
% CARDIOID3 Geometry file defining the geometry of a cardioid.
if nargin == 0
  x = 4; % four segments in boundary
  return
end
if nargin == 1
dl = [0   4   8  12
      4   8  12  16
      1   1   2   2
      0   0   0   0];
  x = dl(:,bs);
  return
end
x = zeros(size(s));
y = zeros(size(s));
if numel(bs) == 1 % bs might need scalar expansion
  bs = bs*ones(size(s)); % expand bs
end
cbs = find(bs < 3); % upper half of cardiod
x(cbs) = s(cbs).^4/512 - 3*s(cbs).^2/16 + 4;
y(cbs) = s(cbs).*(64 - s(cbs).^2).^(3/2)/512;
cbs = find(bs >= 3); % lower half
s(cbs) = 16 - s(cbs); % take the reflection
x(cbs) = s(cbs).^4/512 - 3*s(cbs).^2/16 + 4;
y(cbs) = -s(cbs).*(64 - s(cbs).^2).^(3/2)/512; % negate y
end

Plot the geometry function.

pdegplot('cardioid3','EdgeLabels','on')
axis equal

This analytic geometry looks slightly smoother than the previous versions. However, the difference is inconsequential in terms of calculations.

Geometry Not Proportional to Arc Length

You can write a geometry function where the parameter is not proportional to arc length. This can yield a distorted mesh.

function [x,y] = cardioid4(bs,s)
% CARDIOID4 Geometry file defining the geometry of a cardioid.
if nargin == 0
  x = 4; % four segments in boundary
  return
end
if nargin == 1
  dl = [0    pi/2   pi       3*pi/2
        pi/2   pi     3*pi/2   2*pi
        1      1      1        1
        0      0      0        0];
  x = dl(:,bs);
  return
end
r = 2*(1 + cos(s)); % s is not proportional to arc length
x = r.*cos(s);
y = r.*sin(s);
end

Plot the geometry function.

pdegplot('cardioid4','EdgeLabels','on')
axis equal

The labels are not evenly spaced on the edges because the parameter is not proportional to arc length.

Examine the default mesh for each of the four methods of creating geometry.

subplot(2,2,1)
[p,e,t] = initmesh(@cardioid1);
pdeplot(p,e,t)
title('Polygons')
axis equal
subplot(2,2,2)
[p,e,t] = initmesh(@cardioid2);
pdeplot(p,e,t)
title('Integral')
axis equal
subplot(2,2,3)
[p,e,t] = initmesh(@cardioid3);
pdeplot(p,e,t)
title('Analytic')
axis equal
subplot(2,2,4)
[p,e,t] = initmesh(@cardioid4);
pdeplot(p,e,t)
title('Distorted')
axis equal

While the "Distorted" mesh looks a bit less regular than the other meshes (it has some very narrow triangles near the cusp of the cardioid), all of the meshes appear to be usable.