# LinearSolve used in Iteration is extremely slow

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

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]][] // 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}];
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 = 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. 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$$. • Haven't had a chance to scour your code, but anytime I see linear algebra going slowly, I wonder if it isn't trying to solve symbolically. Otherwise, it should be lickety-split. Dec 20, 2022 at 17:53
• OK, ran it. Is your rationalize statement absolutely needed? Get rid of it and your codes takes as long as it takes to press the key. Dec 20, 2022 at 17:56
• Last comment, if you get rid of the rationalize, the answer you get is within $10^{-9}$ of the answer with rationalize. Dec 20, 2022 at 18:01
• @MikeY Ah! Rationalize, but it was already stored in alpha beforehand, so what's causing LinearSolve to slow down? Dec 20, 2022 at 18:07
• @MikeY Ah, I see what you mean, there's a Rationalize also in my For loop. In any case, I changed to n=1000 and m=5. It's taking a bit of time, I don't have a good grasp of what is fast and slow in this case, does my code have a normal speed for this new setting? Dec 20, 2022 at 18:11

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 = 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 = 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} *)


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.
• Could you run code with n=800, m=5 and plot ListLinePlot[Table[zi[i], {i, 1, m}]] to show how it converges? Dec 22, 2022 at 3:35