# How to solve a matrix dominated by zeros?

I am trying to solve a matrix of this form:

Is there a known algorithm or a method to solve this kind of matrices more efficiently than a normal Gauß elimination method?

I input the diagonals as vectors, so the main diagonal (green) in A, B is the 1st lower diagonal (lower yellow), C = B[::-1] is the 1st upper diagonal (upper yellow), D is the 2nd lower diagonal (lower blue), E = D[::-1] is the 2nd upper diagonal (upper blue), F is the 3rd lower diagonal (lower orange) and G = F[::-1] is the 3rd upper diagonal (upper orange).

The result vector is H = [0,value,value,value,value,value,value,value,0]

How get the solution vector x?

Here is my code:

def fd_8(x, v, w, t, a, b):
"""Implements the shooting method to solve linear second order BVPs

Compute finite difference solution to the BVP but with the 8-order formula!

t should be passed in as an n element array.   x, v, and w should be
either n element arrays corresponding to x(t), v(t) and w(t) or
scalars, in which case an n element array with the given value is
generated for each of them.

USAGE:
u'' = x(t) + v(t) u + w(t) u'
u = fd(x, v, w, t, a, b)

INPUT:
x,v,w - arrays containing x(t), v(t), and w(t) values.  May be
specified as Python lists, NumPy arrays, or scalars.  In
each case they are converted to NumPy arrays.
t     - array of n time values to determine u at
a     - solution value at the left boundary: a = u(t[0])
b     - solution value at the right boundary: b = u(t[n-1])

OUTPUT:
u     - array of solution function values corresponding to the
values in the supplied array t.
"""

# Get the dimension of t and make sure that t is an n-element vector

if type(t) != np.ndarray:
if type(t) == list:
t = np.array(t)
else:
t = np.array([ float( t ) ])

n = len(t)

# Make sure that x, v, and w are either scalars or n-element vectors.
# If they are scalars then we create vectors with the scalar value in
# each position.

if type(x) == int or type(x) == float:
x = np.array([float(x)] * n)

if type(v) == int or type(v) == float:
v = np.array([float(v)] * n)

if type(w) == int or type(w) == float:
w = np.array([float(w)] * n)

# Compute the stepsize.  It is assumed that all elements in t are
# equally spaced.

h = t[1] - t[0];

if (n > 7):
# Construct tridiagonal system; boundary conditions appear as first and
# last equations in system.

#A = -( 1.0 + w[1:n] * h / 2.0 )

A = - (490/(180*(h**2))) - v

A[1:3] = -(2/h**2) - v[1:3]
A[-3:-1] = -(2/h**2) - v[-3:-1]
A[0] = A[-1] = 1.0

B = np.zeros(n-1)
for i in range(n-1):
B[i] = 270/(180*(h**2))
B[0:2] = B[-3:-1] = B[-1] = 1/h**2
B[-1] = 0.0

C = B[::-1].copy()

D = np.zeros(n-2)
for i in range(n-2):
D[i] = - (27/(180*(h**2)))
D[0] = 0.0
D[-3:-1] = D[-1] = 0.0

E = D[::-1].copy()

F = np.zeros(n-3)
for i in range(n-3):
F[i] = 2/(180*(h**2))
F[-3:-1] = F[-1] = 0.0

G = F[::-1].copy()

H = x
H[0] = a
H[-1] = b
'''
# Solve tridiagonal system?? How?

• I didn't read your post. What does "solve a matrix" mean? – Ted Shifrin May 19 at 1:24
• It seems that the dimension of the matrix could be reduced to $7 \times 7$. – zongxiang yi May 19 at 1:50
• solve is have the solution vector.. – ZelelB May 19 at 14:13
• Iterative solvers are usually good for sparse systems, especially those from PDEs and ODEs. You can probably use CG since most finite difference matrice are positive definite. – tch May 27 at 17:06