LinearSolve used in Iteration is extremely slow

Question

I'm trying to solve the nonlinear second order boundary ODE where the method I'm using is Relaxation method for ODEs,

$$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 \label{1} \tag{1}$$

---

Brief info about relaxation method:

Given a boundary ODE,

$$\frac{d^2z}{dx^2} = g(x,z,z'), \quad z(a) = \alpha, \quad z(b) = \beta, \quad x \in [a,b]$$

We replace the dependent variables by finite differences,

$$z''(x) \rightarrow \frac{z_{i+1} - 2 z_i + z_{i-1}}{h^2}, \quad z'(x) \rightarrow \frac{z_{i+1} - z_{i-1}}{2h}, \quad z(x) \rightarrow z_i, \\ h \; \textrm{is the step size}$$

I'll show a simpler ODE to illustrate as an example,

$$z''(x) = -\frac{1}{(1+z)^2} \rightarrow \frac{z_{i+1} - 2 z_i + z_{i-1}}{h^2} = -\frac{1}{(1+z_i)^2}$$

where the BCs are $z(0) = z(1) = 0$

Here is the key step, we rewrite this equation to a new equation called the residual $e_i$ (measure of error),

$$e_i = z_{i+1} - 2 z_i + z_{i-1} + \frac{h^2}{(1+z_i)^2}$$

Taylor expand it,

\begin{equation} e^{\rm{new}}_i(z_{i-1} + \Delta z_{i-1}, z_{i} + \Delta z_{i}, z_{i+1} + \Delta z_{i+1}) = e^{\rm{old}}_i + \frac{\partial e_i}{\partial z_{i-1}} \Delta z_{i-1} + \frac{\partial e_i}{\partial z_{i}} \Delta z_{i} + \frac{\partial e_i}{\partial z_{i+1}} \Delta z_{i+1}\\ = e^{\rm{old}}_i + \alpha_i \Delta z_{i-1} + \beta_i \Delta z_{i} + \gamma_i \Delta z_{i+1} \; (i = 1, \dots n) \end{equation}

We want the newest error to be zero $\textbf{e}^{\rm{new}} = 0$, thus we form the linear equation,

\begin{equation} \textbf{A} \cdot \Delta \textbf{z} = - \textbf{e}^{\rm{old}} \end{equation}

We solve this to find the correction vector $\Delta \textbf{z}$, $\textbf{z}_{old}$ is the initial test value, and $\textbf{z}_{new}$ is the new test value, where $\textbf{z}_{new} = \textbf{z}_{old} + \Delta \textbf{z}$. They call $A$ sparse matrix since it contains many zeros (actually it is tridiagonal) and could be inverted quickly. Last note, only row $i=2$ to $i=n-1$ of the partial derivatives I've really calculated in the code because rows $i=1$ and $i=n$ are $(1,0,\dots,0)$ and $(0,0,\dots,1)$ because they represent the boundary conditions (I don't go into details but that is the case).

---

Going back to eq.\eqref{1}. The code that I've written is divided into three parts: (1) Equation setup (2) Building the matrix $A$ (3) Solving the linear equation then do iteration. For the iteration I've used the Multidimensional Newton's method.

My question is, step (1) & (2) are calculated very quickly, but step (3) is slow, for $n=10$ and iteration $m=2$ it takes about a minute. However, since this is a finite difference approach, I really need a big $n$ as well as several iterations. Can anyone help me reassess the code that I've written? I'm just a beginner in numerical analysis as well as not a Mathematica expert.

(*Equation setup*)
ClearAll["Global`*"]
Needs["VariationalMethods`"]
f = 1 - (z[x]/zh)^(d + 1);(*Blackening factor*)
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;(*2nd order EOM*)
eqsample = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};

(*Building the matrix A*)
a = 0;
b = 1;
n = 10;
h = (b - a)/(n - 1);
alpha = Rationalize[9.306854, 10^-6];
beta = 10^-3;
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}];
eqs = Table[{eqsample} /. rule[[i]], {i, Length[rule]}];
residual = h^2 eqs // Flatten;
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}]
sparseresidual = ShiftMatrix[DFx, Table[i, {i, 0, n - 3}]][[All, n - 4 ;;]];
sparse = Join[{Join[{1}, ConstantArray[0, n - 1]]}, sparseresidual, {Join[ConstantArray[0, n - 1], {1}]}];

(*Solve linear equation then iterate*)
m = 2;
For[j = 0, j <= m, j++, residuals = h^2 eqs;
zi[0] = Join[{alpha}, Reverse[Table[(Rationalize[9.306854, 10^-6] - 10^-3) h i , {i, 1, n - 2}]], {beta}];
zr[j] = MapThread[#1 -> #2 &, {Array[z, Length[Table[i, {i, 1, n}]]], zi[j]}];
DFxmat = sparse /. zr[j];
Residvec = Join[{{0}}, Table[residuals[[i]], {i, 1, n - 2}] /. zr[j], {{0}}]; 
zi[j + 1] = zi[j] + LinearSolve[DFxmat, -Residvec] // Flatten]
zi[m] // N

The initial test value is zi[0]. The sparse matrix $A$ is shown below, the first and last row is discussed in the last sentence of the relaxation method I've written, i.e. $\beta_1 = \beta_M = 1$ and $\gamma_1 = \alpha_M = 0$. Also, the image used $M$ as opposed to my $n$.

Answer 1

Your problem is with defining your sparse matrix in terms of z[i] and then using rules to set their values. This actually takes a really long time. With $n$ rules and $O(n^{2})$ terms in the matrix, you are looking at $O(n^3)$ operations to assign values. Here's a way to work around that without completely rewriting all of your code.

Setting n = 500, the linear solving of your code after getting rid of the Rationalize statement...the timing is 17 seconds. Most of that time is spent setting up your matrix.

(m = 2;
 zi[0] = Flatten@{{alpha}, Reverse[Table[(9.306854 - 10^-3) h i, {i, 1, n - 2}]], {beta}} // N;
 For[j = 0, j <= m, j++, 
    residuals = h^2 eqs;
    zr = MapThread[#1 -> #2 &, {Array[z, Length[Table[i, {i, 1, n}]]], zi[j]}];
    DFxmat = sparse /. zr;
    Residvec = Flatten@{{{0}}, Table[residuals[[i]], {i, 1, n - 2}] /. zr, {{0}}};
    zi[j + 1] = zi[j] + LinearSolve[DFxmat, -Residvec] // Flatten
    ]); // Timing

(* {17.5469, Null} *)

Now the modification to how DFxmat and Residvec get assigned...a third of a second.

(m = 2;
 zi[0] = Flatten@{{alpha}, Reverse[Table[(9.306854 - 10^-3) h i, {i, 1, n - 2}]], {beta}} // N;
 For[j = 0, j <= m, j++, 
    residuals = h^2 eqs;
    DFxmat = sparse /. z[i_] :> zi[j][[i]];
    Residvec = Flatten@{0, (residuals /. z[i_] :> zi[j][[i]]), 0};
    zi[j + 1] = zi[j] + LinearSolve[DFxmat, -Residvec] // Flatten
    ]); // Timing

(* {0.296875, Null} *)

Both answers are the same.

What is going on here, is that I am scanning the sparse term and wherever I find a v[i] term I replace it immediately with zi[j][[i]]. This takes a single pass to do the matrix. Much, much faster than creating a long set of rules and applying them. Note the :> instead of -> which is due to the need to hold off evaluating the right side of the rule until it is used.

with n=1000 and m=5 it now takes 1.8 seconds.