I would like to inquire about a certain approach to numerically solve the equations that result from a finite difference method (FDM) discretization of boundary value problems (BVPs). For concreteness, consider the following problem $$1) \ \frac{d^2\mathbf{y}}{dx^2}=\mathbf{f}(x,\mathbf{y}), x \in [a,b], \ \mathrm{subject \ to} \ \mathbf{y}(a)=\mathbf{\alpha}, \ \mathbf{y}(b)=\mathbf{\beta},$$ where $\mathbf{y}$ is a function $\mathbf{y}:[a,b] \to \mathbb{R}^n$ and $\mathbf{f}$ is some function $\mathbf{f}:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n$.
Let us discretize 1) by introducing a grid $\{x_i=a+ih|i=0,...,N+1\}$, where the step size $h$ is given by $h=(b-a)/(N+1)$, and corresponding samples $\{\mathbf{y}_i=\mathbf{y}(x_i)|i=0,...,N+1\}$. We have $$ \frac{d^2\mathbf{y}}{dx^2}(x_i)=\frac{1}{h^2}(\mathbf{y}_{i+1}+\mathbf{y}_{i-1}-2\mathbf{y}_{i})+\mathcal{O}(h^2).$$ Thus, a 2:nd order accurate discretization of 1) is given by $$2) \ \left\{\begin{aligned}& \mathbf{y}_2+\alpha-2\mathbf{y}_1=h^2\mathbf{f}(x_1,\mathbf{y}_1) \\ & \mathbf{y}_{i+1}+\mathbf{y}_{i-1}-2\mathbf{y}_i=h^2\mathbf{f}(x_i,\mathbf{y}_i), \ i = 2,...,N-1 \\ & \beta+\mathbf{y}_{N-1}-2\mathbf{y}_N=h^2\mathbf{f}(x_N,\mathbf{y}_N)\end{aligned}. \right. $$ Equation 2) is a non-linear system of equations whose solution gives an approximate solution to 1). One way of solving the system 2) would be as follows: Introduce the $n\cdot N$ vector $\mathbf{Y} = [\mathbf{y}_1^\mathrm{T} \dots \mathbf{y}_N^\mathrm{T}]^\mathrm{T}$, where $\mathrm{T}$ denotes transposition. Further define the function $\mathbf{F}$ by $\mathbf{F}:\mathbb{R}^{n\cdot N}\to \mathbb{R}^{n\cdot N}$ with $$\mathbf{F}(\mathbf{Y}) = \left[\begin{aligned}& \hspace{1cm}\mathbf{y}_2+\alpha-2\mathbf{y}_1-h^2\mathbf{f}(x_1,\mathbf{y}_1) \\ & \hspace{1cm} \mathbf{y}_{3}+\mathbf{y}_{1}-2\mathbf{y}_2-h^2\mathbf{f}(x_2,\mathbf{y}_2) \\ &\hspace{3.5cm}\vdots \\ & \mathbf{y}_{N}+\mathbf{y}_{N-2}-2\mathbf{y}_{N-1}-h^2\mathbf{f}(x_{N-1},\mathbf{y}_{N-1}) \\ & \hspace{1cm}\beta+\mathbf{y}_{N-1}-2\mathbf{y}_N-h^2\mathbf{f}(x_N,\mathbf{y}_N) \end{aligned} \right]. $$ Then 2) is equivalent to the system $$3) \ \mathbf{F}(\mathbf{Y})=\mathbf{0},$$ which can be solved numerically using e.g. Newton's method.
My actual question, however, concerns a different approach. For simplicity, take $N=3$. Then 2) takes the explicit form $$4) \ \left\{\begin{aligned}& \mathbf{y}_2+\alpha-2\mathbf{y}_1=h^2\mathbf{f}(x_1,\mathbf{y}_1) \\ & \mathbf{y}_{3}+\mathbf{y}_{1}-2\mathbf{y}_2=h^2\mathbf{f}(x_2,\mathbf{y}_2)\\& \beta+\mathbf{y}_{2}-2\mathbf{y}_3=h^2\mathbf{f}(x_3,\mathbf{y}_3)\end{aligned}. \right. $$ We see that on row $i$, the right hand side only depends on $\mathbf{y}_i$ (and $x_i$) but not on $\mathbf{y}_j$ for $j\neq i$). This allows us to eliminate $\mathbf{y}_2$ in terms of $\mathbf{y}_1$ and $\mathbf{y}_3$ in terms of $\mathbf{y}_2$. Thus we find $$5) \ \mathbf{y}_2=2\mathbf{y}_1-\alpha+h^2\mathbf{f}(x_1,\mathbf{y}_1) $$ and $$5') \ \mathbf{y}_3=-\mathbf{y}_1+2\mathbf{y}_2+h^2\mathbf{f}(x_2,\mathbf{y}_2) =-2\alpha+3\mathbf{y}_1+2h^2\mathbf{f}(x_1,\mathbf{y}_1)+h^2\mathbf{f}(x_2,2\mathbf{y}_1-\alpha+h^2\mathbf{f}(x_1,\mathbf{y}_1))$$ Substituting 5) and 5') into the last equation in 4) yields a single equation in $\mathbf{y}_1$ of the form $$6) \ \mathbf{g}(\mathbf{y}_1)=\mathbf{0}, $$ for a function $\mathbf{g}$ defined by $$\begin{align}7) \ \mathbf{g}(\mathbf{y})=&\beta+3\alpha-4\mathbf{y}-2h^2\mathbf{f}(x_1,\mathbf{y})-2h^2\mathbf{f}(x_2,\mathbf{y}-\alpha+h^2\mathbf{f}(x_1,\mathbf{y}))\\&-h^2\mathbf{f}(x_3,-2\alpha+3\mathbf{y}+2h^2\mathbf{f}(x_1,\mathbf{y})+h^2\mathbf{f}(x_2,2\mathbf{y}-\alpha+h^2\mathbf{f}(x_1,\mathbf{y}))) \end{align}.$$ In principle, solving 6) for $\mathbf{y}_1$ and the using 5) and 5') to solve for $\mathbf{y}_2$ and $\mathbf{y}_3$, respectively, solves the system 4).
QUESTIONS:
- Is this a sensible approach? The obvious advantage of 6) vs. 3) is that 6) has a lower dimensionality ($n$ vs $n\cdot N$). The obvious disadvantage of 6) is the complicated nature of the nested expressions of $\mathbf{f}$.
- Does this approach have a particular name? If so, are there any good references where I can read more about it?
- What would be the best way to solve equation 6) ? Newton's method, or maybe fixed point iteration?
Yes, this is possible, this would make the method a single-shooting method based on the weakly stable central Euler method. If you go that way, you could also do better by using a more stable and higher order method. All the disadvantages of single-shooting methods apply, a possible extreme sensitivity to $x_1$, unknown ranges of inadmissible values for $x_1$ where the forward integration diverges or drastically leaves the domain of the expected solution values, a tendency to move far away from the expected solution if multiple solutions are possible.
An implementation as a multiple-shooting method, that is, solving the whole system as a non-linear system via some approximation to the Newton method, ameliorates or avoids these problems.