Does anyone know where could I find a code (in Matlab or Mathematica, for example) for he Stokes equation in 2D? It has been solved numerically by so many people and referenced in so many paper that I guess someone has had the generous (and in science, appropriate) idea to share it somewhere. Thanks a lot.


3 Answers 3


In Mathematica you can do this on multiple levels. First, let's just try to solve it. We need a geometry, an equation and boundary conditions:

reg = ImplicitRegion[0 <= x <= 2 && 0 <= y <= 0.5 && ! (x >= 1 && y <= 0.1) && ! (x >= 1 && y >= 0.4), {x,y}];
rp = RegionPlot[reg, AspectRatio -> Automatic]

enter image description here

Next, we need the equation and boundary conditions:

stokesFlowOperator = {
Div[{{-1, 0}, {0, -1}}.Grad[u[x, y], {x, y}], {x, y}] + Derivative[1,0][w][x, y], 
Div[{{-1, 0}, {0, -1}}.Grad[v[x, y], {x, y}], {x, y}] + Derivative[0,1][w][x, y], 
Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]};

bc = {DirichletCondition[u[x, y] == 4*0.3*y*(0.5 - y)/(0.41)^2,x == 0.], 
   DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 0 < x < 2], 
   DirichletCondition[w[x, y] == 0., x == 2]};

Top-level solution:

OK, this was the setup now to solve the equation we use NDSolve.

{xVel, yVel, pressure} = NDSolveValue[{stokesFlowOperator == {0, 0, 0}, bc}, {u, v, w}, {x, y} \[Element] reg, Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, w -> 1}}];

And visualize the velocity vector field over the region:

Show[rp, StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] reg]]

enter image description here



An issue raised in the comments, here is an alternative way to visualize the same stream plot:

rmf = RegionMember[reg];
Show[rp, StreamPlot[{xVel[x, y], yVel[x, y]}, {x, 0, 2}, {y, 0, 0.5}, 
  AspectRatio -> Automatic, 
  RegionFunction -> Function[{x, y}, rmf[{x, y}]]]]

which gives the same result.

Low-level access:

Now, say you want to drill in deeper, you certainly can do that and there are many ways you can do that explained here. I show one way. We use NDSolve as an equation pre-processor and extract the finite element data from it. The fact that you can do that means that you have access to every stage of the solution.

Use NDSolve as an equation pre-processor:

{ndstate} = 
 NDSolve`ProcessEquations[{stokesFlowOperator == {0, 0, 0}, bc}, {u, 
   v, w}, {x, y} \[Element] reg, 
  Method -> {"FiniteElement", 
    "InterpolationOrder" -> {u -> 2, v -> 2, w -> 1}}]

This will return you an NDSolve state object.We load the finite element package:


And extract and inspect the finite element data:

femd = ndstate["FiniteElementData"];
sd = ndstate["SolutionData"][[1]];
{"BoundaryConditionData", "FEMMethodData", "PDECoefficientData", \
"Properties", "Solution"}

So the equation pre-procesing created a couple of data structures. Let's extract and use them:

bcd = femd["BoundaryConditionData"];
methodData = femd["FEMMethodData"];
pded = femd["PDECoefficientData"];

These are then the key steps in this low level approach:

dpde = DiscretizePDE[pded, methodData, sd];
dbc = DiscretizeBoundaryConditions[bcd, methodData, sd];

Get the system matrices:

{load, stiffness, damping, mass} = dpde["SystemMatrices"];

You could look at them, modify them or what not.


enter image description here

To apply the boundary conditions use:

DeployBoundaryConditions[{load, stiffness, None, None}, dbc]

Solve the equations:

stationarySolution = LinearSolve[stiffness, load];

Create an interpolating function:

mesh = methodData["ElementMesh"];
offsets = methodData["IncidentOffsets"];
xVel = ElementMeshInterpolation[{mesh}, 
   stationarySolution[[offsets[[1]] + 1 ;; offsets[[2]]]]];

And visualize, say, the x-velocity in a contour plot:

ContourPlot[xVel[x, y], {x, y} \[Element] mesh, 
 AspectRatio -> Automatic, PlotRange -> All, 
 ColorFunction -> ColorData["TemperatureMap"], Contours -> 10]

enter image description here

For an overview of the finite element method in NDSolve have a look here.

Hope this helps.

  • $\begingroup$ Thanks a lot! Do you have any idea why Show[rp, StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] reg]] does not work in my version? It's $\endgroup$ Commented Mar 26, 2015 at 16:09
  • $\begingroup$ If I am not mistaken this was fixed in 10.0.1. $\endgroup$
    – user21
    Commented Mar 26, 2015 at 18:13
  • $\begingroup$ @usumdelphini, I have a 10.0.2 where it works. So just grab the update, if for some reason you can not get the update, I put in a workaround. $\endgroup$
    – user21
    Commented Mar 26, 2015 at 22:24

You can of course implement a simple algorithm in Matlab or similar programs. Or you could use a professional implementation using modern numerical methods that is part of one of the existing finite element packages. For example, here is the Stokes tutorial program of deal.II. (Disclaimer: I wrote this program.) I'm sure there are similar programs for libMesh, FEniCS, etc.


Well, implementing the Stokes equations, as you said, is rather easy. I would suggest you to do it yourself, just by reading the first 20 pages of chapter 7 on Peyret's (1987) "Computational Methods For fluid flow" you will be able to implement such a simple model (just drop the convective term).

  • $\begingroup$ I thank you for your suggestion, but if I asked the question is because I do not want to spend time writing the code. Mainly not because I wouldn't be able to understand how to do it, but because I do not want to deal with debugging, logic and syntax issues at all. So my question is: isn't there a code in MatLab or Mathematica that I could use and modify? $\endgroup$ Commented Mar 21, 2015 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.