Nonlinear Robin boundary condition involving square root

Question

If you have a nonlinear second-order boundary value problem where the domain of the problem is $x \in [a,b]$, the boundary conditions imposed are the Robins condition at $x=a$ and the Dirichlet condition at $x=b$, we can use the finite difference method to solve this. For example,

\begin{equation} \frac{d^2z}{dx^2} = f(z,z',x) \end{equation}

\begin{equation} c_1 z'(a) + c_2 z(a) = c_3 \end{equation}

In terms of finite differences (central for the interior and forward for the left boundary) where the grid points are $i=1 , \dots, n$ and the step size is $h$.

The governing equation is,

\begin{equation} \frac{z_{i + 1} - 2 z_i + z_{i - 1}}{h^2} = f\left(z_i,\left(\frac{z_{i + 1} - z_{i - 1}}{2 h} \right),x_i\right) \quad i=2,\dots,n-1 \end{equation}

The solution $z(x)$ should satisfy the governing equation but it should also satisfy the boundary conditions. So, we need to write down a difference equation at the boundary $x=a$ that satisfies both the boundary conditions and the governing equation.

Taylor expanding $z_2$ and $z_3$ at the point $z_1$,

\begin{align} & z_2 = z_1 + \frac{dz}{dx}\Bigg|_{i=1} h + \frac{d^2z}{dx^2}\Bigg|_{i=1} \frac{h^2}{2!} + \ldots \\ & z_3 = z_1 + \frac{dz}{dx}\Bigg|_{i=1} (2h) + \frac{d^2z}{dx^2}\Bigg|_{i=1} \frac{4h^2}{2!} + \ldots \end{align}

we have,

\begin{equation} \frac{d^2z}{dx^2}\Bigg|_{i=1} = \frac{-z_3 + 8 z_2 - 7z_1 - 6 h \frac{dz}{dx}\Big|_{i=1}}{2 h^2} \end{equation}

From the Robins condition,

\begin{equation} \frac{dz}{dx}\Big|_{i=1} = \frac{c_3 - c_2 z_1}{c_1} \end{equation}

the governing equation at $i=1$ is,

\begin{align} & \frac{d^2z}{dx^2}\Bigg|_{i=1} = f\left(z_1,\frac{dz}{dx}\Big|_{i=1},x_1\right) \end{align}

Thus,

\begin{equation} \frac{-z_3 + 8 z_2 - 7z_1 - 6 h \frac{c_3 - c_2 z_1}{c_1}}{2 h^2} = f\left(z_1,\frac{c_3 - c_2 z_1}{c_1},x_1\right) \end{equation}

Assuming $f$ is not very complicated, i.e. it could be nonlinear but not anything complicated, say a polynomial, then you could collect like terms and maybe have a few nonlinear polynomial terms.

---

In my situation, the nonlinear second order boundary value differential equation is

\begin{equation} z''(x)-\frac{\frac{1}{100} z(x)^4 \left(2 z'(x)^2+12\right)-600 \left(z'(x)^2+1\right)-\frac{3 z(x)^8}{500000}}{20 z(x) \left(10-\frac{z(x)^4}{1000}\right)}=0 \end{equation}

I have a nonlinear Robins boundary condition at $x=a$. \begin{equation} z'(a) + (d-1) \left(1-\left(\frac{z(a)}{z_h}\right)^{d+1}\right) \sqrt{1+\frac{z'(a)^2}{1-\left(\frac{z(a)}{z_h}\right)^{d+1}}}=0 \end{equation}

where $d=3$, $z_h=10$, and $x\in [10^{-8},10^{-1}]$. Also, the Dirichlet condition at $b=10^{-1}$ is $z(b) = 10^{-3}$. $\textbf{TAKE NOTE:}\; z(a) < z_h = 10$

Discretizing the ODE, we have the governing equation for $i = 2, \ldots, n-1$

\begin{equation} \frac{z_{i + 1} - 2 z_i + z_{i - 1}}{h^2}-\frac{\frac{1}{100} z_i^4 \left(2 \left(\frac{z_{i + 1} - z_{i - 1}}{2 h} \right)^2+12\right)-600 \left(\left(\frac{z_{i + 1} - z_{i - 1}}{2 h} \right)^2+1\right)-\frac{3 z_i^8}{500000}}{20 z_i \left(10-\frac{z_i^4}{1000}\right)}=0 \end{equation}

I have tried discretizing the Robins boundary conditions in the same way as above, i.e. by isolating $z'(a)$, but I'm getting solutions that are doubtful to be correct. So, should I do anything here before doing the finite difference? In particular, I'm thinking the square root is causing some issues. Any guidance on this or suggestion if it does not differ from usual methods?

My code is written in Mathematica below,

(*Setup the equation*)
Needs["VariationalMethods`"]
f = 1 - (z[x]/zh)^(d + 1);
L = (Sqrt[1 + (z'[x]^2/f)]/z[x]^d) + (d - 1) (z'[x]/z[x]^d);(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;
eq = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};(*governing equation displayed in post, s[[1,2]] is like the f(z,z',x) in the post*)
bc = z'[x] + (d - 1) (1 - (z[x]/zh)^(d + 1)) Sqrt[1 + (z'[x]^2/(1 - (z[x]/zh)^(d + 1)))];(*Robins condition displayed in post*)

(*Setting up the finite difference and residuals*)
n = 1000;(*Grid points*)
h = (b - a)/(n - 1);(*Step size*)
a = 10^-8;(*a & b are the domain*)
b = 10^-1;
zf = 10^-3;(*zf is the Dirichlet condition at x=b*)
zp = -Sqrt[((d - 1)^2 (1 - (z[x]/zh)^(d + 1))^2)/(1 - (d - 1)^2 (1 - (z[x]/zh)^(d + 1)))] /. {d -> 3, zh -> 10, z[x] -> z[1]};(*zp is the z'[a] solved from the boundary conditions bc*)
rule = Table[{z''[x] -> ((z[i + 1] - 2 z[i] + z[i - 1])/h^2), z'[x] -> ((z[i + 1] - z[i - 1])/(2 h)), z[x] -> z[i]}, {i, 2, n - 1}];(*finite difference rule*)
eqs = Table[{eq} /. rule[[i]], {i, Length[rule]}];(*substitute the finite difference to the governing equation*)
residual = h^2 eqs // Flatten;(*residual of the governing equation*)
residbound = (-z[3] + 8 z[2] - 7 z[1] - 6 zp h) - 2 h^2 s[[1, 2]] /. {d -> 3, zh -> 10, z[x] -> z[1], z'[x] -> zp};(*residual of the Robins condition, s[[1,2]] is like the f(z,z',x) in the post*)

(*Setup the sparse matrix*)
For[i = 2, i <= n - 1, i++, jac[i] = D[residual[[i - 1]], {{z[i - 1], z[i], z[i + 1]}}]]
DFx = Table[jac[i], {i, 2, n - 1}];
ShiftMatrix[mat_, shift_] := Reverse@PadLeft@MapThread[PadLeft[#1, Length[mat] + #2, 0, #2] &, {Reverse[mat], shift}]
jacbound = D[residbound, {{z[1], z[2], z[3]}}];
sparseresidual = ShiftMatrix[DFx, Table[i, {i, 0, n - 3}]][[All, n - 4 ;;]];
sparse = Join[{Join[jacbound, ConstantArray[0, n - 3]]}, sparseresidual, {Join[ConstantArray[0, n - 1], {1}]}];

m = 90;(*Number of iteration*)
zi = 9.5;(*Initial test point for z[1]*)
z0[0] = Join[{zi}, Reverse[Table[((zi - zf)/(b - a)) (i - a) + zf, {i, a + h, b - h, h}]], {zf}];(*Initial test points for the z[i]'s*)

(*Newton's method*)
For[j = 0, j <= m, j++, residuals = h^2 eqs;
DFxmat = sparse /. z[i_] :> z0[j][[I]];
Residvec = {residbound /. z[i_] :> z0[j][[i]], residuals /. z[i_] :> z0[j][[i]], 0} // Flatten; 
z0[j + 1] = z0[j] + 0.22 LinearSolve[N[DFxmat], N[-Residvec]] // Flatten] // AbsoluteTiming

(*Residual error*)
ResidTol = Total[Table[Abs[residuals[[i]]], {i, 1, n - 2}] /. z[i_] :> z0[j][[i]] // Flatten]/n;
Print["Residual Tolerance = ", ResidTol]

Residual Tolerance = 7.090653622*10^-14

range = Range[a, b, h];
list = MapThread[{#1, #2} &, {range, z0[j]}];
zapprox = Interpolation[list, InterpolationOrder -> 10, Method -> "Spline"];

zapprox[a]
9.999999722

zapprox'[a]
0.01664915436

zapprox''[a]
-498.9821196

(*Check if Robins boundary is satisfied, must be equal to 0*)
bc /. {d -> 3, zh -> 10, z[x] -> zapprox[a], z'[x] -> zapprox'[a]}
0.01666025571

(*Check if governing equation is satisfied, must be equal to 0*)
zapprox''[a] - s[[1, 2]] /. {d -> 3, zh -> 10, z[x] -> zapprox[a], z'[x] -> zapprox'[a]}
-0.003923034814

(*Using zapprox[a] in the Robins condition bc as a consistency check for zapprox'[a]*)
sol = (z' + (d - 1) (1 - (z[x]/zh)^(d + 1)) Sqrt[1 + (z'^2/(1 - (z[x]/zh)^(d + 1)))]) /. {d -> 3, zh -> 10, z[x] -> zapprox[a]};

Solve[sol == 0, z']
z' -> -2.222100815*10^-7

As you see, the governing equation and boundary equation at $x=a$ is not 0. Also, the resulting zapprox'(a) is not consistent with the result given by the Robins condition bc if you plug in zapprox(a).

Answer 1

I would approach the problem in more abstract way:

A linear approximation of the Poisson problem including the Robin boundary conditions can already be written as

\begin{align} \mathbf{A} \vec{z} = \vec{f}. \end{align}

Now in your case the right hand side depends on the solution itself making the problem non-linear

\begin{align} \mathbf{A} \vec{z} = f(\vec{z}), \end{align}

which can be written as minimization problem

\begin{align} \mathbf{A} \vec{z} - f(\vec{z}) =: F(\vec{z}) \stackrel{!}{=} 0. \end{align}

The problem can be tackled using the Newton-Raphson method

\begin{align} \vec{z}^{k+1} =\vec{z}^{k} - \mathbf{J}(\vec{z}^k)^{-1} F(\vec{z}^k). \end{align}

with the Jacobian matrix

\begin{align} \mathbf J = \frac{\partial F_i}{\partial z_j}\longrightarrow\mathbf{A} - \text{diag} \left( \left. \frac{\partial f_i}{\partial z_i}\right|_{i=j}\right). \end{align}

Instead of calculating $\mathbf{J}^{-1}$ you can solve this linear system

\begin{align} \mathbf{J}(\vec{z}^k) \vec{b} = F(\vec{z}^k), \end{align}

and update the solution iteratively with

\begin{align} \vec{z}^{k+1} =\vec{z}^{k} - \vec{b}(\vec{z}^k). \end{align}