0
$\begingroup$

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?

Improve this question
$\endgroup$
7
  • 4
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Apr 6, 2020 at 8:15
  • 2
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Apr 6, 2020 at 8:16
  • 4
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Apr 6, 2020 at 8:43
  • 2
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Apr 6, 2020 at 8:46
  • 1
    $\begingroup$ 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. $\endgroup$
    – Henrik Schumacher
    Apr 6, 2020 at 8:47

2 Answers 2

Reset to default
2
$\begingroup$

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]

Figure 1

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – A little mouse on the pampas
    Apr 7, 2020 at 9:56
  • 1
    $\begingroup$ 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? $\endgroup$
    – Alex Trounev
    Apr 7, 2020 at 10:38
  • $\begingroup$ 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. $\endgroup$
    – A little mouse on the pampas
    Apr 7, 2020 at 22:54
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes Ah, sorry, there is a missing part of code. Now should be good. $\endgroup$
    – Alex Trounev
    Apr 7, 2020 at 23:14
  • $\begingroup$ Thank you very much for your great answer. $\endgroup$
    – A little mouse on the pampas
    Apr 7, 2020 at 23:41
1
$\begingroup$

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
  ]
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – A little mouse on the pampas
    Apr 6, 2020 at 9:06

Your Answer

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

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