The Problem:

I have the 4th order Ordinary Differential Equation

$$ \frac{\text{d}^4\theta}{\text{d}\eta^4} +R(\theta-\theta_*)=0 $$

in the interval $0\le\eta\le1$, subject to the boundary conditions

$$ \eta=0: \frac{\text{d}\theta}{\text{d}\eta}=-1 ; \frac{\text{d}^2\theta}{\text{d}\eta^2}=0 $$

$$ \eta=1: \frac{\text{d}\theta}{\text{d}\eta}=0 ; \frac{\text{d}^2\theta}{\text{d}\eta^2}=0 $$

and where $\theta_*$ is to be determined such that the clamping constraint

$$\theta(\eta=0)=0$$ is satisfied. Skipping details, it can be shown that the solution to the differential equation is

$$ \theta=\theta_* + e^{P(\eta-1)}(A\cos P\eta+B\sin P\eta) +e^{-P\eta}(C\cos P\eta+D\sin P\eta) $$

where $P=\frac{R^\frac{1}{4}}{\sqrt{2}}$ and (skipping details again) the constants $A,B,C,D$ and $\theta_*$ can be determined from the boundary conditions and constraint. So far, so good.

The issue:

The solution works beautifully, until $R$ approaches $10^7$ whereupon it breaks down due to what I believe is the stiffness of the differential equation - the difference between the largest and smallest roots of the characteristic equation is of the order of $2P$ ~$R^\frac{1}{4}$. This is also apparent from the original differential equation itself, where as $R$ becomes very large $\theta \rightarrow \theta_*$ which tends to violate the Neumann boundary condition $\frac{\text{d}\theta}{\text{d}\eta}(\eta=0)=-1$. What I find very odd however, is that the breakdown in the analytical solution is manifested not at $\eta=0$, where the Neumann BC is actually satisfied very well, but by blowing up in the vicinity of $\eta=1$. This is evident in the graphic below:

enter image description here

My Question

Given that the analytical solution tends to break down at large $R$, how much confidence can I place in the computed values in the vicinity of $\eta=0$. The Neumann condition at $\eta=0$ certainly seems to be honoured for $R=10^7$, but I'm a bit circumspect about the correctness of the peak value in the second derivative (right plot in the graphic above).

Any advice? Thanks in advance.

Note that in practice, I clamp the value of $\eta$ used to compute $\theta$ and its derivatives at $\eta=1.1-0.1\log_{10}R$, for $R\ge 10^6$

  • $\begingroup$ Out of curiousity: is this equation modelling some real life phenomena (perhaps mechanical bending)? Could it be that for such large values of R the system is better described by another model? Again, I don't doubt that you have considered this possibility and maybe this question is motivated by pure mathematical interest $\endgroup$ – Yuriy S Oct 25 '18 at 20:34
  • $\begingroup$ It's a model of fully developed natural convection in a vertical enclosure with a constant heat flux on one boundary with the other boundary insulated. The 4th order ODE arises from combining the energy and momentum equations with the Boussinesq approximation. $\endgroup$ – Sharat V Chandrasekhar Oct 25 '18 at 20:38
  • $\begingroup$ Yes, most definitely!!! Thanks for the catch. My LaTeX skills are not what they used to be!!! $\endgroup$ – Sharat V Chandrasekhar Oct 25 '18 at 20:39
  • $\begingroup$ @ Sharat V Chandrasekhar, thank you for the clarification, then I would ask a more concrete question: are you sure Boussinesq approximation is valid for R ~ 10^7? $\endgroup$ – Yuriy S Oct 25 '18 at 20:43
  • $\begingroup$ What solution method are you using? It should be something like multiple shooting with adaptive segmentation. $\endgroup$ – LutzL Oct 25 '18 at 20:48

With a proper boundary value solver there is no such boundary layer problem. Using the one from python scipy.integrate the code is

import numpy as np
from scipy.integrate import solve_bvp, odeint
import matplotlib.pyplot as plt

def odesys(t,u,R):
    th = u[0]; th_ast = u[4]; 
    return [ u[1], u[2], u[3], -R*(th-th_ast), 0*th_ast]

def boundary(u0, u1):
    return [ u0[0], u0[1]+1, u0[2], u1[1], u1[2] ]

x = np.linspace(0, 1, 3)

for R in [ 1e3, 1e4, 1e5, 1e6, 1e7 ]:
    res = solve_bvp(lambda t,u: odesys(t,u,R), boundary, x, [-x, 0*x+1, 0*x, 0*x, 0*x-1])
    x_plot = np.linspace(0, 1, 500)
    u_plot = res.sol(x_plot)
    plt.plot(x_plot, u_plot[0], label='R=%.2e'%R)
    plt.plot(x_plot, u_plot[2], label='R=%.2e'%R)
vlabels = [" ", "$\\theta$", "$\\theta''$"]
for k in [1,2]:

Note that $\theta_*$ was included in the state as a constant function with derivative $0$ so that the bvp solver automatically also adapts this value. This system with now 5 state dimensions allows to give all 5 boundary conditions at once.

The resulting diagram is

plot of 5 solutions

  • $\begingroup$ Pretty awesome, I now know whom to ask about numerical methods $\endgroup$ – Yuriy S Oct 25 '18 at 21:42
  • $\begingroup$ Very Nice indeed!! I'm pleased to note that even the broken analytical solution for $R=10^7$ matches the numerical (shouldn't it be the other way around?!) solution near $eta=0$. $\endgroup$ – Sharat V Chandrasekhar Oct 25 '18 at 22:03
  • $\begingroup$ Also not to be pedantic, but the problem only has 4 Boundary Conditions - the natural boundary conditions are the Neumann BCs at $\eta=0$ and $\eta=1$ where the first derivatives of $\theta$ are -1 and 0, respectively. The condition $\theta(0)=0$ is not a true BC, but a constraint that clamps the solution which would otherwise be non-unique. $\endgroup$ – Sharat V Chandrasekhar Oct 26 '18 at 6:29
  • $\begingroup$ For the extended system of $(θ,θ_*)$ it is a boundary condition and makes the solution unique. This system of a 4th and a first order equation could also be formulated as 5th order equation $θ^{(5)}+Rθ'=0$ where you can now pose all 5 conditions as boundary conditions in the more restricted sense. $\endgroup$ – LutzL Oct 26 '18 at 9:52
  • $\begingroup$ For the extended system, yes, I agree. $\endgroup$ – Sharat V Chandrasekhar Oct 26 '18 at 20:48

I've figured out how to avoid the issue with the breakdown of the analytical solution as follows:

From the constraint $\theta(0)=0)$ we have

$$\theta_* = -(Ae^{-P} +C)$$

and from the no slip condition $\theta''(0)=0$, we have


whereupon the remaining boundary conditions yield the following 3 $\times$ 3 system for the constants $A,B$ and $C$

$$ \begin{bmatrix} e^{-P} & 2e^{-P} & -1 \\ \cos P - \sin P\ &(1+e^{-P})\cos P +(1-e^{-P})\sin P & -e^{-P}(\cos P + \sin P)\\ -\sin P & (1-e^{-P})\cos P &e^{-P}\sin P \end{bmatrix} \begin{bmatrix} A \\ B \\ C\end{bmatrix} = \begin{bmatrix} -P^{-1} \\ 0 \\ 0\end{bmatrix} $$

For very large values of $R$ (large $P$) , the associated very small values of $e^{-P}$ were causing roundoff errors with the default double precision in the VBA implementation and yielding erroneous evaluations of $B$ and $C$. Assuming a machine roundoff threshold of $10^{-16}$ (typical for double precision), this occurs at a value of $R=7.4\times 10^6$ which is just about where I noticed the analytical solution breaking down! At such low values of $e^{-P}$, the linear system above is functionally equivalent to

$$ \begin{bmatrix} 0 & 0 & -1 \\ \cos P - \sin P\ &\cos P +\sin P & 0\\ -\sin P & \cos P &0 \end{bmatrix} \begin{bmatrix} A \\ B \\ C\end{bmatrix} = \begin{bmatrix} -P^{-1} \\ 0 \\ 0\end{bmatrix} $$

which gives $C=P^{-1}$, $A=B=0$ which yields the correct solution without any restrictions on $R$ and without recourse to the $\eta$- cutoff hack that I was employing earlier.

So in essence, the solution to the problem that I observed is to just set $C=P^{-1}$, $A=B=0$ for $R>7.4\times 10^6$. As seen in the first graphic below, this does result in any perceptibly abrupt changes in $A$, $B$ or $C$ when this fix is implemented.

enter image description here

The temperature and velocity profiles now look like this

enter image description here


I found an even simpler solution - balancing the 3$\times$3 system prior to inversion for the coefficients $A,B$ and $C$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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