3
$\begingroup$

I have a second order differential equation that I'd like to solve numerically, and then integrate its solution twice to get the parametric equations of a curve. The ODE is:

$2 \kappa''(s)+\kappa^3(s)-\kappa(s)=0$

and the parametric equations are $x(s)=\int \cos (\theta(s))$ and $y(s)=1+\int \sin (\theta(s))$, where: $\theta(s)=\int \kappa(s) \ ds$.

I accomplish this by doing:

    curve = NDSolve[{2 \[Kappa]''[s] + \[Kappa][s]^3 - \[Kappa][s] == 0, \[Kappa][0] == -1, \[Kappa][1] == 1}, \[Kappa], {s, 0, 1} ];]
    \[Theta][t_?NumericQ] := NIntegrate[\[Kappa][s] /. First@curve, {s, 0, t}]
    x[u_?NumericQ] := NIntegrate[Cos[\[Theta][t]], {t, 0, u}]
    y[u_?NumericQ] := 1 + NIntegrate[Sin[\[Theta][t]], {t, 0, u}]
    ParametricPlot[{x[u], y[u]}, {u, 0, 1}]

Question: How do I either fix the arclength of the curve or fix the end points of the curve at $x(1)=0$ and $y(0)=1$?

My first attempt was along the lines of this similar question: Solve differential equation using a integral form boundary condition with a solution by Jens. This eliminates one numerical integration, getting me closer to putting the boundary conditions into the initial NDSolve. (Note: I've cleaned up the code with easier symbols, and a better way to incorporate the boundary conditions from MathX).

    curve = ParametricNDSolve[{k'''[s] + 1/6 k'[s] k[s]^2 - 1/2 k'[s] == 0, k[0] == 0, k'[0] == bc0, k'[1] == bc1}, k, {s, 0, 1}, {bc0, bc1} ];
    x[u_?NumericQ, a_?NumericQ, b_?NumericQ] := NIntegrate[Cos[k[a, b][t] /. curve], {t, 0, u}] 
    y[u_?NumericQ, a_?NumericQ, b_?NumericQ] := 1 + NIntegrate[Sin[k[a, b][t] /. curve], {t, 0, u}]
    ParametricPlot[{x[u, -1, 1], y[u, -1, 1]}, {u, 0, 1}]

EDIT: I'm closer now by continuing the antiderivative procedure out to fourth order. Now I am able to fix $x(1)=0$, but not both $x(1)=0$ and $y(0)=1$:

    curve = ParametricNDSolve[{k''''[s] + 1/12 k''[s] k[s] - 1/2 k''[s] == 0, k'[0] == 0, k[1] == 0, k''[0] == bc0, k''[1] == bc1}, k, {s, 0, 1}, {bc0, bc1} ];
    Manipulate[Plot[Evaluate[k[-1, n][s] /. curve], {s, 0, 1}, PlotRange -> {{0, 1}, {-0.1, 1}}], {n, -1, 1.4, 0.05, Appearance -> "Labeled"}]

EDIT - Question Clarification: Imagine an elastic beam (say, a thin sheet of metal) bent into an arch. One end of the arch is fixed, and the other end of the arch is being twisted (about an axis orthogonal to the planar curve). The ODE above describes the beam's curvature along its length. In addition to boundary conditions that prescribe the beam's curvature at the fixed end and the twisting end, I am looking for a way to ensure that the beam's end points in $x$ and $y$ don't move, and/or that the beam doesn't change length.

Improve this question
$\endgroup$
5
  • $\begingroup$ Where exactly in this code do you want to apply the boundary condition of x[1]==0 and y[0]==1? Everything in your current code (after fixing the syntax error) seems to be working. Do you mean you want to change the range of integration in definition of x and y? $\endgroup$
    – MathX
    Apr 13, 2016 at 16:44
  • $\begingroup$ I don't have working code with those boundary conditions. I would like to study this equation as I vary \[Kappa][1] (which, in the code above is \[Kappa][1] == 1). If you run the code for different values of \[Kappa][1], you'll see that x[1] varies, though I'd like it to remain fixed at x[1]==0. $\endgroup$
    – dpholmes
    Apr 13, 2016 at 17:03
  • $\begingroup$ After @J.M.'s wonderful answer and also because you changed the ODE from what you posted initially, which i solved for, I'm gonna delete my answer. I think you have received the answers you were looking for, right? $\endgroup$
    – MathX
    Apr 13, 2016 at 21:25
  • $\begingroup$ I'm still baffled by what you're doing differentiating the original ODE to allow for additional boundary conditions, but I'm pretty sure that in doing so, you're modifying the original ODE. One derivative adds a constant while two derivatives add a linear term. There will be some boundary conditions which result in the constant or linear term being 0, but in general they won't be. $\endgroup$
    – obsolesced
    Jul 18, 2016 at 8:37
  • $\begingroup$ Being 2nd order, the original ODE admits 2 boundary conditions. If you've chosen to fix the curvature at the end points, \[Kappa][0] and \[Kappa][1], you can't specify any more boundary conditions. If, instead, you want to fix the end points x[1] and y[0], you can't specify \[Kappa][0] and \[Kappa][1]. By the way, I think you mean y[1], since y[0] is fixed at 1 by definition. Also, the arc length is fixed to be 1 by construction: Integrate[Sqrt[x'[s]^2+y'[s]^2],{s,0,1}] == 1, so did you mean fixing the distance between the end points Sqrt[(x[1]-x[0])^2+(y[1]-y[0])^2]? $\endgroup$
    – obsolesced
    Jul 18, 2016 at 8:56

1 Answer 1

Reset to default
3
$\begingroup$

I should confess that I am still slightly perplexed by what you actually want to do, but I'll proceed with the interpretation that you're starting from a Cesàro equation implicitly defined by a second-order ODE (which, BTW, can be solved in terms of elliptic functions, but let's pretend we don't know that), and you want to play with the boundary conditions.

Let me propose the alternative of forming the requisite ODEs from the Frenet-Serret equations augmented with your ODE for the curvature. That is,

k0 = 0; k1 = -3; (* boundary conditions for curvature *)
v1 = {1, 0}; (* unit tangent vector *)
{xx, kk} = {x, κ} /. First @ NDSolve[
                     {x'[s] == t[s], x[0] == {0, 1},
                      t'[s] == κ[s] n[s], t[0] == v1,
                      n'[s] == - κ[s] t[s], n[0] == Cross[v1],
                      2 κ''[s] + κ[s]^3 - κ[s] == 0, κ[0] == k0, κ[1] == k1},
                     {x, κ}, {s, 0, 1}, Method -> "StiffnessSwitching"]

and then

ParametricPlot[xx[s], {s, 0, 1}]

my modest attempt

You will have to figure out how to adjust k0, k1, and v1 to suit your needs. ParametricNDSolve[] might be one possible way to help you in these adjustments.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks @J.M., I really like this framework. I am familiar with this equation as Euler's elastica, and I'm aware of its analytical solution by way of elliptic integrals. Part of the focus behind this question is to have a numerical toolset to go along with the analytical one. $\endgroup$
    – dpholmes
    Apr 14, 2016 at 1:39
  • $\begingroup$ I have edited the original question to clear up the confusion you mention in your answer. In short, the curve's lengths/end-points shouldn't change, since in my case I'm modeling an inextensible elastica. $\endgroup$
    – dpholmes
    Apr 14, 2016 at 1:48
  • $\begingroup$ Silly me, I didn't recognize the elastica at once! Well, the arclength shouldn't be a problem; the thing is set up such that the resulting curve is unit-speed (in this case, the arclength is 1 since that's the upper limit integrated to). I'll think about how to maintain the other endpoint. $\endgroup$
    – J. M.'s eventual burnout ♦
    Apr 14, 2016 at 1:51

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.