# What can I do to solve the numerical solution of this partial differential equation?

I want to solve the following differential equation. It is about modeling the vibration of a linearized elastic rod.

$$\frac{\partial^2} {\partial x^2} ( E I \frac{\partial^2 w} {\partial x^2}) + \rho S \frac{\partial^2 w} {\partial t^2} =0$$

In the above formula, $$E I = 1, \rho S = 1$$, $$w(x,t)$$ is a binary function of $$x$$ and $$t$$.

The boundary and initial conditions are as follows:

$$w(x,t) \Big| _{t=0}=\frac{x^2} {6} (3 - x)$$ $$\frac{\partial w} {\partial t}\Big| _{x=0}=0$$

$$\frac{\partial^2 w} {\partial t^2}\Big| _{x=1}=0$$ $$\frac{\partial^3 w} {\partial t^3}\Big| _{x=1}=0$$

I wrote the following code according to the above conditions:

ClearAll["Global*"]
tau = 10;
L = 1;
Elastic = 1;
Imoment = 1;
ρ = 1;
S = 1;
sol = NDSolveValue[{D[Elastic*Imoment*D[w[x, t], {x, 2}], {x, 2}] +
S*ρ*D[w[x, t], {t, 2}] == 0, w[x, 0] == x^2/6 (3 - x),
D[w[0, t], {t, 1}] == 0,
D[w[L, t], {t, 2}] == 0 D[w[L, t], {t, 3}] == 0},
w[x, t], {x, 0, L}, {t, 0, tau},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100},
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 100, "MinPoints" -> 100,
"DifferenceOrder" -> 2}}, MaxSteps -> 10^6]


But I can't get the numerical solution of $$w(x,t)$$, so I can't draw the vibration image of the first 10 seconds.

What can I do to solve this partial differential equation?

• As I said already in a comment your other post, the boundary condition on the third derivative in time does not make much sense. Actually all three boundary conditions in time do not make sense to me. Also, you need to impose initial conditions for the time derivatives. Commented Apr 6, 2020 at 8:15
• Also, for people new to this, you should give at least a bit of context (it is about modelling the vibration of a linearized elastic rod) and link to you previous post. Commented Apr 6, 2020 at 8:16
• hi. It is simple. You have PDE with order 4 in space and order 2 in time. Hence you need 4 boundary conditions and 2 initial conditions. You list above only 4 "conditions". You need a total of "6 conditions". two for time and 4 for space. You can not have a "condition" with equal or higher derivative than the order of derivative in the PDE. So for time, IC can not have higher than first order derivative, since the time in the PDE has order 2 derivative. Similar for space. Commented Apr 6, 2020 at 8:43
• What I said is: The partial differential equation that you wrote does not make sense. As we do not know which problem you actually want to solve, we cannot tell how to modify the boundary conditions. As for the initial conditions: If you mean the cantilever to be at rest in the beginning, then use D[w[x, t], {t, 1}] == 0 /. t -> 0 is what you should. Commented Apr 6, 2020 at 8:46
• I guess from your other post that you want the left and of the cantilever to be clamped. That would be {w[0, t] == 0, D[w[x, t], {x, 1}] == 0 /. x -> 0}. I also guess that the right end should be free. I am not 100% sure, but I think the correct boundary conditions for that would be {D[w[x, t], {x, 2}] == 0 /. x -> L, D[w[x, t], {x, 3}] == 0 /. x -> L}. Other users will know that better than me. Commented Apr 6, 2020 at 8:47

To extend an answer by Henrik Schumacher and my answer here we conseder numericaly finished code

tau = 3;
L = 1;
Elastic = 1;
Imoment = 1;
\[Rho] = 1;
S = 1; A = 1/10; Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
sol = NDSolveValue[{S \[Rho] D[w[x, t], {t, 2}] +
D[Elastic Imoment D[w[x, t], {x, 2}], {x, 2}] == 0,
w[x, t] == A x^2/6 (3 - x) /. t -> 0,
D[w[x, t], {t, 1}] == 0 /. t -> 0, w[x, t] == 0 /. x -> 0,
D[w[x, t], {x, 1}] == 0 /. x -> 0,
D[w[x, t], {x, 2}] == 0 /. x -> L,
D[w[x, t], {x, 3}] == If[t <= 10^-4, -A, -A Exp[-10 t]] /.
x -> L}, w, {x, 0, L}, {t, 0, tau},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> False,
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 80, "MinPoints" -> 80,
"DifferenceOrder" -> 4}},
EvaluationMonitor :> (currentTime = t;)];]

Plot3D[sol[x, t], {x, 0, L}, {t, 0, tau}, ColorFunction -> "Rainbow",
AxesLabel -> Automatic, Mesh -> None, PlotPoints -> 50]


• But I ran this code in Mathematica version 12.1, but I get a lot of error information, and I can't get the image shown in your answer. Commented Apr 7, 2020 at 9:56
• This code is debugged for version "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)". What messages did you receive in version 12.1? Commented Apr 7, 2020 at 10:38
• I received the following information in v12.1: NDSolveValue::deqn: Equation or list of equations expected instead of If[t<=1/10000,-A,-A exp(-10 t)] in the first argument {w^(0,2)(x,t)+w^(4,0)(x,t)==0,w(x,0)==1/60 (3-x) x^2,w^(0,1)(x,0)==0,w(0,t)==0,w^(1,0)(0,t)==0,w^(2,0)(1,t)==0,If[t<=1/10000,-A,-A exp(-10 t)]}. NDSolveValue::dsvar: 0.0000204285714 cannot be used as a variable. General::stop: Further output of NDSolveValue::dsvar will be suppressed during this calculation. Commented Apr 7, 2020 at 22:54
• @PleaseCorrectGrammarMistakes Ah, sorry, there is a missing part of code. Now should be good. Commented Apr 7, 2020 at 23:14
• Thank you very much for your great answer. Commented Apr 7, 2020 at 23:41

This might work.

sol = NDSolveValue[{
S \[Rho] D[w[x, t], {t, 2}] + D[Elastic Imoment D[w[x, t], {x, 2}], {x, 2}] == 0,
w[x, t] == x^2/6 (3 - x) /. t -> 0,
D[w[x, t], {t, 1}] == 0 /. t -> 0,
w[x, t] == 0 /. x -> 0,
D[w[x, t], {x, 1}] == 0 /. x -> 0,
D[w[x, t], {x, 2}] == 0 /. x -> L,
D[w[x, t], {x, 3}] == 0 /. x -> L
},
w,
{x, 0, L}, {t, 0, 10},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100},
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 100, "MinPoints" -> 100, "DifferenceOrder" -> 2}
},
MaxSteps -> 10^6
]
`
• Thank you very much for your answer, but the image of w at t = 0 is obviously not a parabola w=x^2/6 (3 - x). When t = 10, the amplitude of the cantilever beam is obviously too large, which is far greater than the initial amplitude L^2/6 (3 - L) at the first right end. Commented Apr 6, 2020 at 9:06