Numerical Computing with Geometry

Rob Knapp
Paritosh Mokashi

Wolfram Research

Join the Conversation #WolframTechConf

Example of Workflow

The only additional steps you typically need for doing numerical computations with geometric regions is to define a region and indicate how that region is a part of the solution you seek.

The region can be anything that is RegionQ.  In this presentation we'll show for the most part regions coming from symbolic representations, but you can also get regions using Import, or even constructing meshes manually.

The most common way that a region is a part of the solution is in the definition of the domain for the problem variables.  The notation  vars∈Ω is used to indicate that the variables vars are restricted to being elements of the region Ω.

Define a region symbolically

Ω=Annulus[]

Of course, one still needs to specify what is to be solved.  For this example, I’ll choose Laplace’s equation:

Define a differential equation

pde = Laplacian[u[x,y],{x,y}]==0;

Depending on the type of problem you wish to solve, you may have to give more information about the interaction of the region with the solution.  PDEs such as this need boundary conditions for uniqueness.

Define a boundary condition that the solution is equal to x×y at every point where both x and y are positive

Γ=DirichletCondition[u[x,y]==x y,x>0 && y> 0];

... and now you can ask NDSolve to solve the problem

Solve the problem

usol = NDSolveValue[{pde,Γ},u,{x,y}∈Ω]

The next stage is to visualize the solution.  The best way of doing this depends on the problem you want to solve, and so on, and so forth.

The plotting functions typically work with the variables restricted to a region as well, so that is a convenient way to visualize this solution

Visualize the solution

Plot3D[usol[x,y],{x,y}∈Ω,PlotRange->All]

And then you may want to refine the solution depending on what you see .... part of the purpose of this presentation is to give you the knowledge to do this.

Under the hood

For solutions with geometric regions, the general symbolic specification of a region needs to be simplified so that the numerical computations we all know and love can work with the simpler form.  Usually this simplification is done in one of two ways:

Discretizing the region into smaller, simpler subregions where the solution of the subproblem on each is simpler (typically the set of subregions is called a mesh)

Expressing the region in terms of symbolic formulae that may be used in combination with numerical computations.

Symbolic Analysis

What is done is very problem dependent, but using some general methods, this can reduce the problem to one on which standard methods can be used

Numerically find the minimum of a function over the region Ω

NMinimize[{ Sin[3(x^2+y^2)]+Cos[7 x y]},{x,y}∈Ω]

Plot3D[Sin[3(x^2+y^2)]+Cos[7 x y],{x,y}∈Ω]

Get a description of the region in terms of inequalities

imp = Region`RegionProperty[Ω,{x,y},"ImplicitDescription"]

Solve the problem expressed in terms of inequality constraints

NMinimize[{Sin[3(x^2+y^2)]+Cos[7 x y],imp},{x,y}]

It is not always so easy.

Numerically approximate the integral of  (x^2)+(y^2)over the region Ω

approx = NIntegrate[x^2+y^2,{x,y}∈Ω]

This was (chosen to be) a relatively simple integral that can be done exactly with just a little thought:

Exactly integrate (x^2)+(y^2)over the region Ω

exact = Integrate[x^2+y^2,{x,y}∈Ω]

2 Pi Integrate[r^3,{r,1/2,1}]

Show the error:

exact - approx

The second way the exact result was computed gives one idea for how it might be done.  However, as nice and simple as that solution was for this particular example, that approach does not work in general.  An approach that does work in general is to use the cylindrical algebraic decomposition (CAD) on the region

Get a symbolic description of the region in terms of inequalities and find the CAD

Reduce[imp,{x,y}]

RegionPlot[Ω,Epilog->{Thickness[0.01],Map[Line[{{#,-1},{#,1}}]&,{-1,-1/2,1/2,1}]}]

Compute the integral using CAD

t1=NIntegrate[x^2+y^2,{x,-1,-1/2},{y,-Sqrt[1-x^2],Sqrt[1-x^2]}]

t2=NIntegrate[x^2+y^2,{x,1/2,1},{y,-Sqrt[1-x^2],Sqrt[1-x^2]}]

t3 = NIntegrate[x^2+y^2,{x,-1/2,1/2},{y,-Sqrt[1-x^2],-Sqrt[1- 4x^2]/2}]

t4 = NIntegrate[x^2+y^2,{x,-1/2,1/2},{y,Sqrt[1-4x^2]/2,Sqrt[1-x^2]}]

t1+t2+t3+t4

%-exact

Note this problem can also be solved by discretization:

Numerically approximate the integral of  (x^2)+(y^2)over the region Ω using the finite element method

NIntegrate[x^2+y^2,{x,y}∈Ω,Method->"FiniteElement"]

%-exact

Discretization

You can actually see the simplification of the region used in the first example:

Get the mesh from the solution

mesh = usol["ElementMesh"]

mesh["Wireframe"]

This is a fairly coarse mesh that is used by default to try and get at least a reasonable solution.

You can specify that a finer mesh be used

Solve the problem with a finer mesh

usolf = NDSolveValue[{pde,Γ},u,{x,y}∈Ω,Method->{"FiniteElement","MeshOptions"->{MaxCellMeasure->0.001}}]

Plot3D[usolf[x,y],{x,y}∈Ω,PlotRange->All]

mesh = usolf["ElementMesh"]

mesh["Wireframe"]

Discretization: A cautionary tale:

A tempting thing to do when you have many problems to solve with the same geometry is to discretize the geometry once and use that repeatedly.  This is certainly a reasonable thing to do, but must be used with caution.

Discretize the region

(Subscript[Ω, d])=DiscretizeRegion[Ω,Method->"Continuation"]

Compute the integral on the discrete region

dapprox = NIntegrate[x^2+y^2,{x,y}∈(Subscript[Ω, d])]

dapprox-exact

So this is much less accurate than the solution over the symbolic region -- once you have discretized a region, the solution can only be as good as the discretization you actually obtained.