5
$\begingroup$

I am trying to reproduce results from this paper, which involves solving a fourth-order differential equation.

After setting the various parameters:

l = 5.;
L = 25.;
rad = 0.4;
splay = Tan[17. Pi/180];

and defining the so-called "pressure" function:

press[x_] := Max[ 0.85 (1 - x/6.), 0.0]

I used NDSolve to solve the differential equation, imposing 4 boundary conditions:

sol = First[
  NDSolve[ { 
    l^3  D[ R[s], {s, 4}] - (L - s) D[R[s], {s, 2}] + D[R[s], s] - 
      press[R[s]] == 0, 
    R[0.0] == rad, 
    Derivative[1][R][0.0] == splay , Derivative[2][R][L] == 0.0, 
    Derivative[3][R][L] == 0.0}, {R}, {s, 0, L}]]

The resulting solution:

rap[s_] = R[s] /. sol

looks absolutely fine to the eye. The problem occurs when I look at one of the properties obtained from the solution, which the authors term "elasticity":

elastic[s_] = l^3 D[rap[s], {s, 4}]

which is proportional to the 4th derivative of the solution. It has more or less the same behavior as the result shown in the paper, but has jagged oscillations superimposed on it. Below I compare the authors' result on the left with the Mathematica result on the right - in both cases they are the curves shown in black.

Comparison of the paper's results, with those I calculated using Mathematica

Clearly taking the 4th derivative numerically will amplify any irregularities in the solution. Is there any way to obtain a cleaner result for the elasticity?

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ Use exact values. Like 4/10 instead of 0.4 and add WorkingPrecision -> 30 as an option of NDSolve. $\endgroup$ Commented Oct 19, 2024 at 16:20

2 Answers 2

Reset to default
6
$\begingroup$

Using only exact values like 4/10 instead of 0.4 or 17 instead of 17. and NDSolve[..., WorkingPrecision -> 30].

LogPlot[elastic[s], {s, 0, 25}, PlotRange -> {10^-4, 10}]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Ahhh, so it's as simple as just increasing the WorkingPrecision... Should I delete the question for triviality? $\endgroup$ Commented Oct 19, 2024 at 16:33
  • $\begingroup$ @ClaraDíazSanchez I am fine with such posts from new contributors. $\endgroup$ Commented Oct 19, 2024 at 16:42
6
$\begingroup$

There are machine-precision solutions. The problem is taking the fourth derivative of an NDSolve solution. The Automatic interpolation order is low, and differentiating the interpolating function increases the noise.

Setting InterpolationOrder -> All gives an interpolating function that is smooth up to the order of the interpolation method. That order, in this case is 5 or 6 (it varies over the course of the integration. Or use Method -> "ExplicitRungeKutta" with InterpolationOrder -> All. It works nicely on this problem and gives an order 9 solution.

Another approach is to solve for the 3rd derivative. The first derivative of a variable in an NDSolve solution is usually smooth.

Similar problem is address here: questions about InterpolationOrder

This shows the difference in noise between two InterpolationOrder settings:

ClearAll[solmeth];
mem : solmeth[iorder_, meth_] := mem =
   First[
    NDSolve[{l^3 D[R[s], {s, 4}] - (L - s) D[R[s], {s, 2}] + 
        D[R[s], s] - press[R[s]] == 0, R[0.0] == rad, 
      Derivative[1][R][0.0] == splay, Derivative[2][R][L] == 0.0, 
      Derivative[3][R][L] == 0.0}, {R}, {s, 0, L}
     , InterpolationOrder -> iorder, Method -> meth]];

Plot[Evaluate@{
   R''''[s] /. solmeth[Automatic, Automatic],
   R''''[s] /. solmeth[All, Automatic]},
 {s, 10, 13},
 Mesh -> (R["Coordinates"] /. solmeth[Automatic, Automatic]), 
 MeshStyle -> Magenta, 
 PlotLegends -> 
  HoldForm /@ {InterpolationOrder -> Automatic, 
    InterpolationOrder -> All}]

Noisy plot with InterpolationOrder -> Automatic, smooth plot with InterpolationOrder ->All

Here's a refactoring of the code to make the 3rd derivative of R a separate variable R3. It also makes a smooth plot of elastic[].

{sol, sol3} = First[NDSolve[{
     R3[s] == D[R[s], {s, 3}],
     l^3 R3'[s] - (L - s) D[R[s], {s, 2}] + D[R[s], s] - press[R[s]] == 0,
     R[0.0] == rad, Derivative[1][R][0.0] == splay, 
     Derivative[2][R][L] == 0.0, R3[L] == 0.0},
    {R, R3}, {s, 0, L}
    ]];
rap[s_] = R[s] /. sol;
elastic[s_] = l^3 R3'[s] /. sol3;

LogPlot[elastic[s], {s, 0, 25}, PlotRange -> {10^-4, 10}]
smooth decreasing log plot of elastic[s]

In part, the high-order derivative in this problem is a bit easier because the differential order of the ODE is also high.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for the nice explanation. I thought that the numerical differentiation was amplifying noise, but I didn't know about raising the order of the interpolation function. $\endgroup$ Commented Oct 20, 2024 at 8:48

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.