# 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,

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

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

$$$$\textbf{A} \cdot \Delta \textbf{z} = - \textbf{e}^{\rm{old}}$$$$

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}];
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$$.

• 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. Commented 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. Commented 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. Commented 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? Commented 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? Commented 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[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} *)


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? Commented Dec 22, 2022 at 3:35