# How to use Finite elements to solve an initial value ODE with NDSolve?

This is on 11.3, windows 7

I have not used Mathematica FEM much at all. So sorry for this basic question on using it to solve a basic second order initial value ODE.

I want to use NDSolve but force it to use FEM to solve a time dependent initial value ODE. (spring/damped system).

My understanding is that, to use FEM, one just needs to change the initial conditions from y[0]==0,y'[0]==0 to use DirichletCondition[y[t] == 0, t == 0] and NeumannValue[0, t == 0], and then use NDSolve as before, but also add the Method -> {"FiniteElement"} as an option.

Is this how one tells NDSolve to use FEM? I am doing something wrong in what follows, since I get wrong answer from NDSolve when I did the above. So I think my initial conditions are not specified correctly.

## Problem

Solve y''[t] + y'[t] + 3 y[t] == Sin[t] with y[0]==0,y'[0]==0

## NDSolve

ClearAll[y,t];
ode  = y''[t]+y'[t]+3y[t]==Sin[t];
ic   = {y[0]==0,y'[0]==0};
sol  = NDSolve[{ode,ic},y,{t,0,20}];
Plot[Evaluate[y[t]/.sol],{t,0,20},AxesOrigin->{0,0},PlotRange->All]

I want to get same solution as above, but want to force NDSolve to use FEM.

## NDSolve with FEM

ClearAll[y,t];
ic1  = DirichletCondition[y[t]==0,t==0];
ic2  = NeumannValue[0,t==0];(*this is not even needed*)
ode  = y''[t]+y'[t]+3y[t]==Sin[t]+ic2;
sol  = NDSolve[{ode,ic1},y,{t,0,20},Method->{"FiniteElement"}];

Plot[Evaluate[y[t]/.sol],{t,0,20},AxesOrigin->{0,0},PlotRange->All]

Which is not correct. I noticed that I can't just write

ClearAll[y,t];
ic  = {y[0]==0,y'[0]==0};
ode = y''[t]+y'[t]+3y[t]==Sin[t];
sol = NDSolve[{ode,ic},y,{t,0,20},Method->{"FiniteElement"}]

As this gives an error. So that is why I changed initial conditions to use DirichletCondition

I think my error is in the "boundary" conditions settings. But I do not know where and how to fix it. The problem is that removing NeumannValue[0,t==0] still gives the same solution. My understanding is that NeumannValue==0 is the default always, that is why removing it makes no change to the solution.

I also understand that FEM is typically used for static problems (or steady state, no time dependence). So the DirichletCondition and NeumannValue typically used on space and not on time. So here I am treating "time" as "space" since I do not know what else to do.

So how to solve the above using NDSolve (or NDSolveValue) but using FEM?

• your last case "works" for me (ie no error), however the result is wrong, simply ignoring the derivative b.c. (v10.1) May 10 '18 at 17:17
• @george2079 thanks. On 11.3, I get error message for last case, also error when I try to plot it. May 10 '18 at 17:20
• Initial values are not boundary conditions. The purpose of FEM is spatial discretization. May 13 '18 at 5:17

Well, I'm not quite familiar with FEM theory, but according to this comment from user21:

It is important to realize that NeumannValue[0, whatever] does not contribute anything ever. It is taken out at parser level. Now, let's assume you have NeumannValue[something, whatever] and DirichletCondition[u==someting, whatever] then the DirichletCondition will trump the NeumannValue.

So ic2 in your 2nd sample is simply ignored, and the actual b.c.s are

$$y(0)=0, \ y'(20)=0$$

This can be verified by

ic1 = DirichletCondition[y[t] == 0, t == 0];
ode = y''[t] + y'[t] + 3 y[t] == Sin[t];
sol = NDSolve[{ode, ic1}, y, {t, 0, 20},
Method -> {"FiniteElement", MeshOptions -> MaxCellMeasure -> 0.001}][[1]];
bctraditional = {y[0] == 0, y'[20] == 0};

Plot[Evaluate[y[t] /. {sol, soltraditional}], {t, 0, 20}, AxesOrigin -> {0, 0},
PlotRange -> All, PlotStyle -> {Automatic, {Red, Dashed}}]

So, how to circumvent this? The only solution I can think out at the moment is transforming the ODE to a 1st order system so the Neumann condition becomes Dirichlet condition and won't be ignored anymore:

odemodified = z'[t] + y'[t] + 3 y[t] == Sin[t];
ic2modified = DirichletCondition[z[t] == 0, t == 0];
odeauxiliary = z[t] == y'[t];
sol = NDSolve[{odemodified, odeauxiliary, ic1, ic2modified}, {y, z}, {t, 0, 20},
Method -> {"FiniteElement"}];

Plot[Evaluate[y[t] /. sol], {t, 0, 20}, AxesOrigin -> {0, 0}, PlotRange -> All]

BTW, though I've transformed the ODE manually here, it can be done automatically with the solutions under this post.

As to the 3rd sample, it fails because "FiniteElement" method can't handle b.c. y'[0] == 0. When "FiniteElement" is chosen, Neumann b.c. and Robin b.c. can only be introduced with NeumannValue, at least now. (This is frustrating I should say. See this post for example. )

• Why do you think this is frustrating. y'[0]==0 is only equivalent to a NeumannValue in a very special case. If y'[0]==0 would map to a NeumannValue that would confuse people that understand what NeumannValue does. I am not going to re-implement that. We had that and it caused more confusion than it helped. May 12 '18 at 9:57
• @user21 By "frustrating" I'm refering to the "Neumann b.c. and Robin b.c. can only be introduced with NeumannValue" part i.e. manually transforming Neumann and Robin b.c. in… Er… traditional form to a NeumannValue may be not easy (as we've encountered in the last link). If it's confusing to make NDSolve do automatic transforming, perhaps we can have a function like ToNeumannValue? May 12 '18 at 10:51
• How would a function ToNeumannValue look like? What would be it's arguments? The point is that a NeumannValue depends on the equation it is associated with y'[0] does not depend on the equation. And that is by design. May 13 '18 at 5:15
• @user21 How about ToNeumannValue[boundaryCondition, associatedEquation]? May 13 '18 at 5:20
• @user21 The 2 arguments are the boundary condition in traditional form and the associated equation, the output is the equation with a NeumannValue[…] inside, for example ToNeumannValue[y'[0] == 0, y''[t] + y'[t] + 3 y[t] == Sin[t]] outputs y''[t] + y'[t] + 3 y[t] == Sin[t] + NeumannValue[0, t == 0]. (OK, perhaps we should call it ToEquationWithNeumannValue. ) May 13 '18 at 5:44

Indeed, one can solve this ODE with finite elements, but currently, the deployment of the boundary conditions has to be done by hand. I use piecewise linear finite elements, because I am more familiar with them than with second order elements.

Let's start by setting up the ODE and its boundary conditions:

Needs["NDSolveFEM"]
ν = 1;
β = 1;
γ = 3;
f = Sin;
dir = 2.;
neu = 0.;
ode = ν y''[t] + β y'[t] + γ y[t] == f[t];
ic = {y[0] == dir, y'[0] == neu};
n = 229;
L = 20;

Using variables β, γ, f, dir, neu, etc. allows us to to see how the following code can be generalized without leaving OP's example.

Now, let's generate a random 1D mesh and use Mathematica's finite element facilities to get the matrices for the weak formulation of our system:

SeedRandom[20180511];
R = ToElementMesh[
(# - #[[1, 1]]) (L/(#[[-1, 1]] - #[[1, 1]])) &@
Accumulate[RandomReal[{0.1, 1}, {n, 1}]],
"MeshElements" -> Line[Partition[Range[n], 2, 1]]
];

vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{y}, {t}}];
sd = NDSolveSolutionData[{"Space"} -> {R}];

cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{{\[Nu]}}},
"MassCoefficients" -> {{1}},
"ConvectionCoefficients" -> {{{\[Beta]}}},
"ReactionCoefficients" -> {{\[Gamma]}},
];

mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];

The usual process of NDSolve for FEM would be to issue a call to DiscretizedBoundaryConditionDataand DeployBoundaryConditions to interweave the matrices for the boundary conditions with the system matrix. That's what we have to do by hand, now. First, let's grab the system matrices as they are without the boundary conditions deployed.

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

Constraining the first degree of freedom (the value of y[0] at the left boundary) to dir can be done by ignoring the first row of the system matrix stiffness. Constraining y'[0] implies that the second degree of freedom (y[h] with h being the diameter of the first mesh cell) must be set to dir + h neu. However, we must not delete the second row of stiffness for it delivers the defining equation for the third degree of freedom. So, we have n-2 values of y to determine, but we have left n-1 equations. This can resolved by testing the weak formulation of the ODE only by those functions that also vanish at the right boundary of the domain. This erases the last row of A. Moreover, we have to add a certain correction to the right hand side that stems from our knowledge of the values of the solution on the first two mesh vertices. Here is how we get the corrected linear system. Since it is banded, we simply can employ LinearSolve with method specialized for banded matrices:

n = Length[stiffness];
A = stiffness[[2 ;; -2, 3 ;;]];
b -= With[{h = R["Coordinates"][[2, 1]] - R["Coordinates"][[1, 1]]},
stiffness[[2 ;; -2]].SparseArray[{1 -> dir, 2 -> dir + h neu}, n]
];
yFEM = Join[{dir, dir + h neu}, LinearSolve[A, b, Method -> "Banded"]];

Finally, let's see how the solution compares to the solution obtained from NDSolve's ODE solver:

g1 = ListLinePlot[Transpose[{Flatten[R["Coordinates"]], yFEM}],
PlotRange -> All,
AxesOrigin -> {0, 0},
PlotStyle -> Directive[ColorData[97][2], Dashed, Thick]
];
ClearAll[y, t];
sol = NDSolve[{ode, ic}, y, {t, 0, 20}];
g2 = Plot[Evaluate[y[t] /. sol], {t, 0, 20},
AxesOrigin -> {0, 0},
PlotStyle -> Directive[Thick],
PlotRange -> All
];
Show[g2, g1]

It's almost perfect, isn't it?

## Word of warning

Using this approach (time-discretization with piecewise-linear functions tested against piecewise-linear functions) for parabolic PDE's will quite likely disappoint you: This discretization has the tendency to get instable if the largest time step is not significanly smaller than the square of the smallest mesh cell diamater in the space domain. This is why Petrov-Galerkin schemes (either piecewise-linear functions tested against piecewise-constant functions or piecewise-constant functions tested against piecewise-linear functions) were invented.

The original problem is an initial value problem, where you specify $y(0)$ and $y'(0)$.

Most FEMs are used for boundary value problems, where you should specify all boundary conditions, not just one: here the boundary is $\{0\}\cup \{10\}$ so you should have one Dirichlet or Neumann (or Robin) condition at $0$, and another one at $10$. You can see for instance that

NDSolveValue[{D[y[t], t, t] + D[y[t], t] + 3*y[t] - Sin[t] ==
NeumannValue[500, t == 10], DirichletCondition[y[t] == 0, t == 0]}, y, {t, 0, 10}, Method -> "FiniteElement"]

perfectly works.

Or, if you really want to solve the IVP with FEM (and not a BVP) you should use a least square process, according to Daniel Nunez:

When considering an IVP, the differential operator is either non-self adjoint or nonlinear. It is never self-adjoint. Hence, the only F.E. method that can guarantee positive definite coefficient matrices for all IVPs is least squares process.

• I am sorry, but I am not following you. Where did the "500" number you used there come from? And the plot generated above does not match the plot from NDSolve when not using FEM as you can see yourself. Also, you did not use Method->{"FiniteElement"} then how does one know that FEM is being used? It is not clear just because one uses DirichletCondition then it will use FEM? I know that FEM is used typically for steady states problems, but it also can be used for time dependent problems. My question is how to do that. May 10 '18 at 22:08
• And thanks to the link. It was useful. But the person there says at the end F.E. method can be used to solve IVP. But more reading is needed. May 10 '18 at 22:45
• @Nasser "……one uses DirichletCondition then it will use FEM?" Yes, this can be verified with NDSolveProcessEquations: ic = DirichletCondition[y[t] == 0, t == 0]; ode = y''[t] + y'[t] + 3 y[t] == Sin[t]; NDSolveProcessEquations[{ode, ic}, y, {t, 0, 20}] As you can see, the output is a …"SteadyState"… object, which is only used for "FiniteElement". May 11 '18 at 10:39
• @Nasser 1) 500 is arbitrary, just to show that if you impose a condition on the other boundary, that works. 2) I had removed `Method -> "FiniteElement" and forgot to put it back, edited. 3) I understand your question, I agree that it can be solved by some FEM, the purpose of my post was to explain that it is not something any FEM can solve, and that very likely the current implementations of FEM in Mathematica can't solve it as it is (which seems confirmed by the two new answers). I know I'm not providing a workaround, sorry, but the other two answers do. May 11 '18 at 15:14