0
$\begingroup$

I am trying to perform numerical integration using the following code:

lagrangian = 
  1/(6 \[Pi]) Sec[r]^3 Sqrt[4/(1 + Cos[r]^2 + Cos[r]^4) + (x1'[r])^2];

Needs["VariationalMethods`"]

eqofmotion = 
 Last[EulerEquations[lagrangian, x1[r], r][[1]] // 
     Numerator] /. {x1'[r] -> 1/y'[x], x1''[r] -> -y''[x]/y'[x]^3} /. 
  r -> y[x]

eq = eqofmotion y'[x]^3 // Simplify

R[rc_] := 
 NDSolveValue[{eq == 0, y[0] == rc, y'[0] == 0}, y, {x, -8, 8}]

Rp[rc_] := D[R[rc][x], x]

xp[rc_] := 1/Rp[rc]

rc[lq_] := 
 rc1 /. Quiet@
   FindRoot[R[rc1][lq/2] == Pi/2 - 10^-3, {rc1, .0001}, 
    PrecisionGoal -> 4, AccuracyGoal -> 3, Evaluated -> False]

rc[1]

En[lq_?NumericQ] := -2 NIntegrate[(
   Sec[\[Chi]]^3 Sqrt[
    4/(1 + Cos[\[Chi]]^2 + Cos[\[Chi]]^4) + (xp[\[Chi]])^2])/(
   6 \[Pi]), {\[Chi], \[Pi]/2, rc[lq]}]

If I want to get the output of, say, En[1], I get errors. How to resolve this issue? Am I missing something here?

Improve this question
$\endgroup$
6
  • $\begingroup$ You should really test each step on its own before adding more steps that uses the last step. This way you know where the problem is. I just did R[1] which is your first very call to NDSolve and got an error NDSolveValue::ndsz: At x == 0.46930048260328944`, step size is effectively zero; singularity or stiff system suspected. $\endgroup$
    – Nasser
    Commented Oct 20 at 15:06
  • $\begingroup$ @Nasser I have just modified this answer for my specific problem of finding the integration. I thought it was just a warning message and, hence, ignored it. $\endgroup$
    – codebpr
    Commented Oct 21 at 1:06
  • $\begingroup$ It is not clear what do you try to solve? $\endgroup$
    – Alex Trounev
    Commented Oct 21 at 3:50
  • $\begingroup$ @AlexTrounev Basically, I am trying to reproduce Figure 2 of this paper, for which I need to solve their integral number 9 and 10. $\endgroup$
    – codebpr
    Commented Oct 21 at 4:16
  • $\begingroup$ @codebpr Could you post the figure you computed? $\endgroup$
    – Alex Trounev
    Commented Oct 22 at 7:41

1 Answer 1

Reset to default
1
$\begingroup$

I tried a different approach, as given in this answer by Alex, and now the integral evaluates fine.

lagrangian = 
 1/(6 \[Pi]) Sec[r]^3 Sqrt[4/(1 + Cos[r]^2 + Cos[r]^4) + (x1'[r])^2]

Needs["VariationalMethods`"]

eqofmotion = 
 Last[EulerEquations[lagrangian, x1[r], r][[1]] // 
     Numerator] /. {x1'[r] -> 1/y'[x], x1''[r] -> -y''[x]/y'[x]^3} /. 
  r -> y[x]

eq = eqofmotion y'[x]^3 // Simplify

lag = lagrangian /. {x1'[r] -> 1/y'[x], x1''[r] -> -y''[x]/y'[x]^3} /. 
  r -> y[x]

ClearAll[R1, R1p, lag, lq, rc]

R1[rc_, lq_] := 
 NDSolveValue[{eq == 0, y[0] == rc, y'[0] == 0}, y, {x, -lq/2, lq/2}]

R1p[rc_, lq_] := D[R1[rc, lq][x], x]

rc[x1_] := 
 rc1 /. Quiet@
   FindRoot[R1[rc1, x1][x1/2] == Pi/2 - 10^-3, {rc1, .0001}, 
    PrecisionGoal -> 4, AccuracyGoal -> 3, Evaluated -> False]

ClearAll[integral]

integral[lq_] := 
 lag /. y[x] -> R1[rc[lq], lq][x] /. y'[x] -> R1p[rc[lq], lq]

En[lq_] := -2 NIntegrate[integral[lq], {x, -lq/2, lq/2}]
Improve this answer
$\endgroup$
2
  • $\begingroup$ I don't understand what they have shown in figure 2, since both integrals $E, E_0$ diverge. $\endgroup$
    – Alex Trounev
    Commented Oct 22 at 7:51
  • $\begingroup$ @AlexTrounev maybe they have used a trick like taking a very small value like $\epsilon$ added in the expressions or something similar. $\endgroup$
    – codebpr
    Commented Oct 22 at 11:44

Your Answer

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

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