2
$\begingroup$

I would like to retrieve the step size information of an InterpolatingFunction returned from NDSolve while using the Finite Element Method (FEM) (Method->{"FiniteElement"}).

I will illustrate the issue with a very simple NDSolve example where the default options are used:

w2 = 6;
m = 2;
T = 2.0;
sol = NDSolveValue[{q'[t] == \[Zeta][t], \[Zeta]'[t] + w2*Sin[q[t]] ==
     0, q[0] == Pi/3., \[Zeta][0] == 0}, {q, \[Zeta]}, {t, 0, T}]

Usually, I can get the step size information by either doing the following:

sol[[1]]["Coordinates"] // First // Differences

or the following

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
InterpolatingFunctionCoordinates[sol[[1]]]

ISSUE WITH FEM

However, if I calculate the same solution using the FEM option from NDSolve:

w2 = 6;
m = 2;
T = 2.0;
sol = NDSolveValue[{q'[t] == \[Zeta][t], \[Zeta]'[t] + w2*Sin[q[t]] ==
     0, DirichletCondition[{q[t] == Pi/3., \[Zeta][t] == 0}, 
    t == 0]}, {q, \[Zeta]}, t \[Element] Line[{{0}, {T}}], 
  Method -> {"FiniteElement"}]

I can no longer obtain the information I require. For example,

sol[[1]]["Coordinates"] // First // Differences

returns the following:

{NDSolve`FEM`ElementMesh[{{0., 2.}}, {NDSolve`FEM`LineElement["<" 20 ">"]}]}

Now, one could naïvely believe that the mesh contains 20 elements equally spaced. However, I am pretty sure that this is not true. How can the required information be retrieved?

$\endgroup$

1 Answer 1

4
$\begingroup$

sol is 2nd order mesh. You can see this with the folllowing code :

   (sol[[1]])["ElementMesh"] ["MeshOrder"]  

2

In 1D, the elements are segments of lines. Because the order is 2, there are 3 points per segments. The points indices are :

 (sol[[1]])["ElementMesh"] ["MeshElements"]  

{NDSolveFEMLineElement[{{1, 2, 22}, {2, 3, 23}, {3, 4, 24}, {4, 5, 25}, {5, 6, 26}, {6, 7, 27}, {7, 8, 28}, {8, 9, 29}, {9, 10, 30}, {10, 11, 31}, {11, 12, 32}, {12, 13, 33}, {13, 14, 34}, {14, 15, 35}, {15, 16, 36}, {16, 17, 37}, {17, 18, 38}, {18, 19, 39}, {19, 20, 40}, {20, 21, 41}}]}

The corresponding locations are :

(sol[[1]])["ElementMesh"] ["Coordinates"]  


{{0.}, {0.1}, {0.2}, {0.3}, {0.4}, {0.5}, {0.6}, {0.7}, {0.8}, {0.9}, \
{1.}, {1.1}, {1.2}, {1.3}, {1.4}, {1.5}, {1.6}, {1.7}, {1.8}, {1.9}, \
{2.}, {0.05}, {0.15}, {0.25}, {0.35}, {0.45}, {0.55}, {0.65}, {0.75}, \
{0.85}, {0.95}, {1.05}, {1.15}, {1.25}, {1.35}, {1.45}, {1.55}, \
{1.65}, {1.75}, {1.85}, {1.95}}  

EDIT

Some clarifications

The fact that the mesh order is 2 is confusing. It becomes clear (I hope) if you do the same thing with mesh order 1 :

w2 = 6;
m = 2;
T = 2.0;
sol = NDSolveValue[{q'[t] == \[Zeta][t], \[Zeta]'[t] + w2*Sin[q[t]] ==
      0, DirichletCondition[{q[t] == Pi/3., \[Zeta][t] == 0}, 
     t == 0]}, {q, \[Zeta]}, t \[Element] Line[{{0}, {T}}], 
   Method -> {
"FiniteElement", "MeshOptions" -> {"MeshOrder" -> 1}
}];
(sol[[1]])["ElementMesh"]["MeshOrder"]
(sol[[1]])["ElementMesh"]["MeshElements"]
(sol[[1]])["ElementMesh"]["Coordinates"]  

1

{NDSolveFEMLineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 10}, {10, 11}, {11, 12}, {12, 13}, {13, 14}, {14, 15}, {15, 16}, {16, 17}, {17, 18}, {18, 19}, {19, 20}, {20, 21}}]}

{{0.}, {0.1}, {0.2}, {0.3}, {0.4}, {0.5}, {0.6}, {0.7}, {0.8}, {0.9}, {1.}, {1.1}, {1.2}, {1.3}, {1.4}, {1.5}, {1.6}, {1.7}, {1.8}, {1.9}, {2.}}

It seems that in the case of a mesh order 2, the first 21 elements are the same as the 21 elements of order 1, and the new elements beyond 21 (0.05, 0.15 ...)are a addition to implement the order 2 mesh. I have always observed this kind of construction, though it is not documented.

$\endgroup$
4
  • $\begingroup$ This looks like what I am looking for. However I have trouble to interpret the output. On the one hand, we have a list that goes from 0. to 2. in the first 21 elements, but then it goes from 0.05 to 1.95 in the subsequent 20 elements. What does this mean? What are the coordinates in the end? The first 21 elements? $\endgroup$
    – Meclassic
    Nov 16, 2022 at 1:39
  • $\begingroup$ This is interesting. From this, I deduce that the step size is constant with the FEM by default. I would have thought that step size is variable as with the finite differences methods of NDSolve... $\endgroup$
    – Meclassic
    Nov 16, 2022 at 16:23
  • 2
    $\begingroup$ In FEM the mesh is determinated uniquely from geometry (adaptatives meshes are not implemented yet in Mma). Here the "geometry" is a simple segment of line (which is the time from 0 to 2), so it's no wonder that the mesh is regular. There's no reason to refine the mesh anywhere. The main use of FEM is for spatial discretisation, not temporal discretisation like here. Hence the term "step size" is innapropriate : there are no steps. $\endgroup$
    – andre314
    Nov 16, 2022 at 17:34
  • $\begingroup$ It is perfectly possible to create a custom mesh (irregular if you wish) and pass it to NDSolve as an option. $\endgroup$
    – andre314
    Nov 16, 2022 at 17:40

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.