# Are there any more time-efficient ways to solve a matrix equation than with SparseArray and LinearSolve?

I am trying to solve a matrix equation, with the matrix being tridiagonal. The problem I am having is that a message keeps appearing that I don't have enough memory. Here is the code:

n = 10109;(*matrix dimension*)

a = SparseArray[{Band[{1, 1}] -> -2, Band[{1, 2}] -> 1,Band[{2, 1}] -> 1}, {n, n}](*the tridiagonal matrix*)

f[x_] := Piecewise[{{x + 1, x <= 0}, {-x + 1, x > 0}}](*initial condition, just for constructing the right-hand side vector*)

g = Table[f[i], {i, -1 + 2/(n + 1), 1 - 2/(n + 1), 2/(n + 1)}];(*the vector*)

u = LinearSolve[a, g];(*the solution I need*)


Actualy what I need is only

Norm[u,1]


and

u[5055].


-
It might make sense to write an explicit Gaussian elimination tailored to banded matrices. If LinearSolve is running out of memory it may be doing something less efficient than that. (I don't know why it would do otherwise, I'm just guessing based on the reported symptom.) –  Daniel Lichtblau May 3 at 19:34
@DanielLichtblau Interestingly, the exact solution can be obtained in (perhaps) acceptable 30 seconds by breaking it up into ff=LinearSolve[a];ff[g]. –  Jens May 3 at 20:39

Update

The speed issue in this case was actually a simple matter: the matrices had only exact numbers, and that slowed the solution down to the extent that it became impractical.

But when I converted the matrix a and the vector g to machine precision, the solution was found almost instantaneously:

n = 10109;(*matrix dimension*)a =
SparseArray[{Band[{1, 1}] -> -2, Band[{1, 2}] -> 1,
Band[{2, 1}] -> 1}, {n, n}];(*the tridiagonal matrix*)

f[x_] := Piecewise[{{x + 1, x <= 0}, {-x + 1,
x > 0}}](*initial condition,just for constructing the right-hand \
side vector*)

g = Table[
f[i], {i, -1 + 2/(n + 1), 1 - 2/(n + 1), 2/(n + 1)}];(*the vector*)


Here is the simple solution based on converting to machine precision:

vv = LinearSolve[N[a], N[g]];

ListPlot[vv]


The absolute error is now much smaller than with my initial approach based on a differential-equation approximation.

diff2 = a.vv - g;

Max[Abs[diff2]]

(* ==> 5.82634*10^-9 *)


If exact numbers are needed

In case you don't want only machine-precision results, you will have to accept a more time-consuming calculation. But it can still be done:

ff = LinearSolve[a];

Timing[ww = ff[SparseArray[g]];]

(* ==>  {27.6843, Null} *)


Leaving out the SparseArray in the last command doesn't make a big difference, but it also doesn't hurt:

diff3 = a.ww - g;

Max[Abs[diff3]]

(* ==> 0 *)


This is therefore an exact solution, requiring less than half a minute of patience. It works by breaking the solution process into two steps. First get a LinearSolveFunction object that is independent of the right-hand side g, and then apply that function to the specific g. This sequence is still kind of slow, but much faster than the single LinearSolve[a,g] (which I never ran to completion).

I'll keep the previous approach below, because it's perhaps useful to observe the special nature of the matrix in this question.

Previous approximate solution

The matrix that you're dealing with is so large that it's best treated in a continuum approximation. This means you go from the discrete solution vector to a continuos function which you can find by solving a differential equation that is the continuous form of the matrix equation.

In your case, that's really easy to do because the matrix a is the finite-difference approximation to a second derivative operator. So here is how the problem can be solved.

Define a continuous right-hand side for our equation in terms of the function f, and call it h. It has the same domain 1,...n as represented by the vector g which it is supposed to replace. Then solve the differential equation that is equivalent to the matrix equation:

h[x_] := f[-1 + 2/(n + 1) + 2 (x - 1)/(n + 1)]

s[x_] = y[x] /.
First@DSolve[D[y[x], x, x] == h[x], y[x], x] /. {C[1] -> 0,
C[2] -> 0}


$$\begin{cases} \frac{x^3}{30330} & x \le 5055 \\ -\frac{x^3}{30330}+x^2-5055 x+8517675 & \text{True} \\ \end{cases}$$

Here, two integration constants C[1], C[2] arise, corresponding to homogeneous solutions. I choose them to be zero.

Now test the solution by sampling the solution of DSolve at discrete points and applying the matrix to it:

u = Table[s[x], {x, 1, n}];

ListPlot[Most[a.u]]


Compare this to the ListPlot of the original right-hand side vector:

ListPlot[g]


The output looks the same, except that I had to throw out the right-most element of the vector a.u from the differential equation because boundary effects cause it to deviate too much. To be more quantitative, here are the differences:

diff = Most[a.u - g];

N[Max[Abs[diff]]]

(* ==> 0.0000659413 *)


Within this error estimate, we can be confident that the following results answer your questions:

N[Norm[u, 1]]

(* ==> 7.53367*10^10 *)

u[[5055]]

(* ==> 8517675/2 *)


Here is what the solution vector looks like:

ListPlot[u]


This solution looks very different from the exact vector obtained earlier. But this is only one of the possible solutions to the differential equations, because a second-order differential equation always leads to two integration constants which I chose arbitrarily. To make the continuos equation yield the same solution plot as the discrete equation, I have to choose C[1] -> 0 and C[2] -> -(774869930101/306605970).

-