Breakdown of Analytical Solution to 4th order ODE

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: 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).

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$$

• 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 – Yuriy S Oct 25 '18 at 20:34
• 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. – Sharat V Chandrasekhar Oct 25 '18 at 20:38
• Yes, most definitely!!! Thanks for the catch. My LaTeX skills are not what they used to be!!! – Sharat V Chandrasekhar Oct 25 '18 at 20:39
• @ 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? – Yuriy S Oct 25 '18 at 20:43
• What solution method are you using? It should be something like multiple shooting with adaptive segmentation. – 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; th_ast = u;
return [ u, u, u, -R*(th-th_ast), 0*th_ast]

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

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.subplot(1,2,1)
plt.plot(x_plot, u_plot, label='R=%.2e'%R)
plt.subplot(1,2,2)
plt.plot(x_plot, u_plot, label='R=%.2e'%R)
vlabels = [" ", "$$\\theta$$", "$$\\theta''$$"]
for k in [1,2]:
plt.subplot(1,2,k)
plt.legend()
plt.xlabel("$$\eta$$")
plt.ylabel(vlabels[k])
plt.grid();
plt.show()


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 • Pretty awesome, I now know whom to ask about numerical methods – Yuriy S Oct 25 '18 at 21:42
• 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$. – Sharat V Chandrasekhar Oct 25 '18 at 22:03
• 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. – Sharat V Chandrasekhar Oct 26 '18 at 6:29
• 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. – LutzL Oct 26 '18 at 9:52
• For the extended system, yes, I agree. – 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

$$D=Be^{-P}$$

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. The temperature and velocity profiles now look like this I found an even simpler solution - balancing the 3$$\times$$3 system prior to inversion for the coefficients $$A,B$$ and $$C$$.