I'm trying to solve numerically the following IVP nonlinear second-order differential equation: $$f''(x)=A\left(1+f'(x)^2\right)^{3/2}+\frac{\left(1+f'(x)^2\right)}{f(x)}$$

where $A$ is a constant and the problem has the initial conditions (IVP) $f(0)=R$ and $f'(0)=0$, where $R$ is also a constant.

I'm used to applying the RK4 method for first-order differential equations, but is it possible to solve it by the RK4 method as well? If so could someone provide me a good reference to understand this specific application for method? I also open to any other useful numerical method that can solve this IVP. Thanks in advance.


2 Answers 2


I had worked on this problem before and am glad to share it with you.

Recognize the slope $\phi=y'$ and main curvature $\kappa_1:$

$$ \dfrac{y''}{(1+y'^2)^{\frac32}}=A +\dfrac{1}{y\sqrt{1+y'^2}} \tag1 $$

Let $y'= \tan \phi,\text{so that }\cos \phi =\dfrac{ 1}{\sqrt{1+y'^2}} =u \;$

$$\kappa_1= A+ \dfrac{u}{y} \tag2 $$

With proper sign

$$ \kappa_1+\kappa_2 = A = 2 H \tag3 $$

where $\kappa_2$ is second curvature perpendicular to the principal meridional direction.

It is an important surface of constant mean curvature CMC when internal pressure acts on a soap bubble as a minimal surface area of revolution for a given enclosed volume.

Btw It can be formulated using calculus of variations nicely.

Having a constant mean curvature=H, they are also known as CMC Delaunay surfaces.

You can directly integrate to first order ODE. Its alternate form first degree:

$$\dfrac { \sin \phi\; d \phi}{dy}-\dfrac{u}{y} \to -\dfrac { du}{dy}-\dfrac{u}{y} = 2H \tag4 $$


$$ u= - H y+\dfrac{c}{y}=\dfrac{1}{\sqrt{1+y'^2}} \tag 5 $$

$$ \boxed{\cos \phi = \dfrac{c}{y}- H y} \tag 6 $$ is the differential equation of DeLaunay surfaces. Vertical tangent at $ y=\sqrt{c/H}.$

Characteristics of shape:

When $ H= 0, $ it is a Catenary; when $ H=1 , $ it is a Circle.

The meridians $(x,y)$ are calculated numerically from above ODEs and plotted below using RK4 available with Mathematica. Initial values $ R=1, y_0'= 0$

enter image description here

An interesting fact about Unduloids: If a conic rolls on a straight line then the Unduloids form as loci of their foci.

The second order RK4 is handled similar to first order, one more initial condition should be supplied.

  • $\begingroup$ @ Mate0us Forcelini: Let me know after you obtain such unduloid plots.. $\endgroup$
    – Narasimham
    Commented Jan 15, 2021 at 20:23
  • $\begingroup$ Thank you very much for your reply! Unfortunately, I just saw your solution now, and I used the RK4 directly to solve the system. This problem came for me dealing with the uplift capacity of plate anchors. I'll review what I did and see how your explanation fits on it. Best regards. $\endgroup$ Commented Jan 25, 2021 at 18:23
  • $\begingroup$ Glad it was helpful. $\endgroup$
    – Narasimham
    Commented Jan 25, 2021 at 19:06
  • $\begingroup$ But am not clear how plate anchor deformation can be significant when plates are not thin compared to span between anchor bolts. $\endgroup$
    – Narasimham
    Commented Jan 25, 2021 at 19:07
  • $\begingroup$ Actually, the present ODE describes the slip surface (the surface where the discontinuous velocity field occurs) that minimizes a functional form for the anchor pullout force. When I say plate anchors I refer to offshore shallow foundations. $\endgroup$ Commented Jan 25, 2021 at 19:17

You can just write your equation as a system of two first order equations and apply RK to that system.

$$ \begin{cases} f' = &g \\ g' = &A(1+g^2)^{3/2} + \dfrac{1+g^2}{f}\\ \end{cases} $$

with initial conditions $f(0) = R, g(0)=0$. If the method is appropriate to this particular problem is another matter. There seems to be a singularity near $x=1$.

  • $\begingroup$ Thank you Pierre. I just noticed that singularity, where it happens depends on the scalar $A$. Do you think is a problem related to the method employed or is the differential equation itself? $\endgroup$ Commented Jan 15, 2021 at 15:24
  • 1
    $\begingroup$ @MateusForcelini Yes, it depends on $A$ and $R$, I just took $A=1$ and $R=10$ in my experiment. The singularity is not related to the method, it prevails over different time steps and numerical methods. $\endgroup$ Commented Jan 15, 2021 at 16:26

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .