# Solve PDE with consraints

I am trying to solve the following problem of the free fall dynamics under gravity of a inextensible horizontal string attached at its end, in a 2D vertical plane. If I'm right, that is the governing equations, for $$x=[x1,x2]\in\mathbb{R}^2$$ and $$s$$ the curvilinear coordinate along the string:

\begin{align} \|x'(s,t)\| &=1 &\text{(inextensibility)}\\ \rho \ddot x(s,t) &= \begin{bmatrix} 0 \\ -\rho g\end{bmatrix} + \lambda(s,t) x'(s,t)^\top x''(s,t) x'(s,t) &\text{(dynamics)}\\ x(0,t)&=[0,0], x(1,t)=[0.5,0] & \text{(boundary cond.)}\\ x(s,0)&=[0,0],\quad \dot x(s,0) = [0,0] & \text{(initial cond.)} \end{align}

Implementation:

x[s_, t_] = {x1[s, t], x2[s, t]};
rho = g = 1;
const = D[x[s, t], s] . D[x[s, t], s] - 1 == 0;
eq1 = rho*D[x1[s, t], {t, 2}] == lambda[s,t]*D[x[s, t], s] . D[x[s, t], {s, 2}] D[x1[s, t], s];
eq2 = rho*D[x2[s, t], {t, 2}] == -rho*g + lambda[s,t]*(D[x[s, t], s] . D[x[s, t], {s, 2}]) D[x2[s, t], s];
bc = {x1[0, t] == 0, x1[1, t] == 0.5, x2[0, t] == 0, x2[1, t] == 0};
ic = {x1[s, 0] == 0, (D[x1[s, t], t] /. t -> 0 ) == 0};


Then,

NDSolve[Flatten@{eq1, eq2, const, bc, ic}, {x1[s, t], x2[s, t],
lambda[s,t]}, {s, 0, 1}, {t, 0, 10}]


returns

The maximum derivative order of the nonlinear PDE coefficients for the Finite Element Method is larger than 1. It may help to rewrite the PDE in inactive form.

So I changed the code to (note that I added some NeumannValue):

eq1 = rho*D[x1[s, t], {t, 2}] == lambda[s, t]*D[x[s, t], s] . Inactive[D][D[x[s, t], {s, 1}], {s, 1}] D[x1[s, t], s] + NeumannValue[0, t == 0];
eq2 = rho*D[x2[s, t], {t, 2}] == -rho*g + lambda[s,t]*(D[x[s, t], s] . Inactive[D][D[x[s, t], {s, 1}], {s, 1}]) D[x2[s, t], s] + NeumannValue[0, t == 0];
bc = {x1[0, t] == 0, x1[1, t] == 0.5, x2[0, t] == 0, x2[1, t] == 0}
ic = {x1[s, 0] == 0, x2[s, 0] == 0, lambda[s, 0] == 0};


and now I have:

Any idea on how to overcome this?

Edit As per xzczd's link, I rewrote the equations using Inactive[Grad] instead of Inactive[D] which is presently not implemented in the FEM package.

eqs = Thread[rho*D[x[s, t], {t, 2}] == {0, -rho*g} + lambda[s,t]
*(D[x[s, t], s] . Inactive[Grad][D[x[s, t], {s, 1}], {s}])
*D[x[s, t], s] + NeumannValue[0, t == 0]];
NDSolveValue[Flatten@{eqs, const, ic}, {x1[s, t], x2[s, t], lambda[s, t]},
{s, 0, 1}, {t, 0, 10}]


I now get the NDSolve:dgsvars error: "The differentiation variables {s} given for Inactive[Grad] should be the spatial independent variables {s,t}". How not to make it believe t is a space variable?

• 1. Inactive@D[...] isn't the correct syntax, should be Inactive[D][...]. 2. Even if you write it right, FiniteElement won't be able to handle it at least for now: mathematica.stackexchange.com/q/217169/1871 Jul 2, 2022 at 10:03
• 3. t is not spatial variable, and you've forgotten 2 i.c.s. Jul 2, 2022 at 11:31
• @xzczd 3. Still no luck with eqs = Thread[rho*D[x[s, t], {t, 2}] == {0, -rho*g} + lambda[s, t]*(D[x[s, t], s] . Inactive[Grad][D[x[s, t], {s, 1}], {s}])* D[x[s, t], s]];, dirichlet = {DirichletCondition[x[s, t] == {0, 0}, s == 0], DirichletCondition[x[s, t] == {0.5, 0}, s == 1]}; and ic = {x1[s, 0] == 0, x2[s, 0] == 0, (D[x1[s, t], t] /. t -> 0) == 0, (D[x2[s, t], t] /. t -> 0) == 0, lambda[s, 0] == 0}: NDSolveValue returns The dependent variable in {x1,x2}=={0,0} in the b.c. DirichletCondition[{x1,x2}=={0,0},s==0] needs to be linear Jul 2, 2022 at 11:49
• 4. It's clear that equations in vector form isn't supported by DirichletCondition, use the traditional one instead. 5. I don't think you can have terms like Inactive[Grad][D[x[s, t], {s, 1}], {s}] in code, see this post for more info: mathematica.stackexchange.com/q/225711/1871 Jul 2, 2022 at 12:13
• Also, a quick test via FDM suggests the convergency of the system is rather bad, according to my (limited) experience, a well-posed IBVP won't be like this. Are you sure the system is correct? Jul 2, 2022 at 13:44

To get the integration started you need to specify sufficient initial conditions:

x[s_, t_] = {x1[s, t], x2[s, t]};
rho = g = 1;
const = D[x[s, t], s] . D[x[s, t], s] - 1 == 0;
eq1 = rho*D[x1[s, t], {t, 2}] ==
lambda[s, t]*D[x[s, t], s] . D[x[s, t], {s, 2}] D[x1[s, t], s];
eq2 = rho*D[x2[s, t], {t, 2}] == -rho*g +
lambda[s, t]*(D[x[s, t], s] . D[x[s, t], {s, 2}]) D[x2[s, t], s];
bc = {x1[0, t] == 0, x1[1, t] == 0.5, x2[0, t] == 0, x2[1, t] == 0};
ic = {x1[s, 0] == 0, (D[x1[s, t], t] /. t -> 0) == 0,
x2[s, 0] == 0, (D[x2[s, t], t] /. t -> 0) == 0,
lambda[s, 0] == 0
};

Monitor[NDSolve[
Flatten@{eq1, eq2, const, bc, ic}, {x1[s, t], x2[s, t],
lambda[s, t]}, {s, 0, 1}, {t, 0, 10},
EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])]
, currentTime]


This still gives messages and I did not run this to the end; but this should be a starting point. Forget FEM for this.