1
$\begingroup$

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

$\endgroup$
8
  • $\begingroup$ You may simplify your problem in a first step and check if the checks are fullfilled pointwise, e.g. using two Dirichlet BC's and a linear right hand side. If this does not hold, then something else is wrong. $\endgroup$
    – ConvexHull
    Jan 7, 2023 at 14:39
  • $\begingroup$ @ConvexHull Take note that $z(a) < z_h = 10$. Taking $z(a) \in (8, 9.999)$ would be reasonable. $\endgroup$
    – mathemania
    Jan 7, 2023 at 14:40
  • 1
    $\begingroup$ As a total aside, it should be "Robin", not "Robins" -- en.wikipedia.org/wiki/Victor_Gustave_Robin $\endgroup$ Jan 8, 2023 at 16:37
  • $\begingroup$ A "quick and dirty" approach is solving it by shooting, integrating the ODE from the end x=b where Dirichlet BC is set to the end x=a where the Robin BC is set. Probably something better than that can be proposed, if we know the strategy for integration in x. How would you solve this equation if Dirichlet BC were set on both ends? $\endgroup$ Jan 10, 2023 at 0:41
  • $\begingroup$ It is (also) a simple nonlinear 1D Poisson equation. Or am I missing something? $\endgroup$
    – ConvexHull
    Jan 10, 2023 at 2:57

1 Answer 1

1
$\begingroup$

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}

$\endgroup$
6
  • $\begingroup$ That is exactly how I solved the problem (I added a Mathematica the code with comments), however, I'm not sure if writing the Robins condition in the same way is correct, i.e. like the simpler setup I wrote at the start of my post. $\endgroup$
    – mathemania
    Jan 7, 2023 at 13:43
  • $\begingroup$ Of course, your problem simply boils down to a linear discretization of the Poisson problem using Robin boundary conditions. $\endgroup$
    – ConvexHull
    Jan 7, 2023 at 13:43
  • $\begingroup$ You mean it is correct to use the same way of discretizing as in the simpler setup that I wrote? But, in the last part of my code you can see that the solution I got doesn't really satisfy the governing equation and the boundary condition, i.e. it's not equal to 0. $\endgroup$
    – mathemania
    Jan 7, 2023 at 13:46
  • $\begingroup$ I didn't say anything that $z(x=a) = 0$, what I'm saying is that looking at the Robins condition where all terms are on one side so that it is an equation (something = 0). I have written my code in such a way too so that if you check everything, the result should give 0 (of course machine precision zero), i.e. it satisfies the Robins condition. $\endgroup$
    – mathemania
    Jan 7, 2023 at 14:08
  • $\begingroup$ Ok, I see the problem. In general, the consistence checks should be fullfilled, pointwise (using a FD method). However, I am no Mathematica user, so I am not able to debug the code. I must think about it a bit. $\endgroup$
    – ConvexHull
    Jan 7, 2023 at 14:13

Your Answer

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

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