Confusing finite difference results with numpy and scipy sparse for simple ODE

Question

I'm having trouble getting the solution to a simple ODE using finite difference. The equation I'm trying to solve is more complicated, but even for a sin function I'm getting different results with scipy.sparse and numpy.linalg. Neither is particularly close to the real solution, even for 20000 x values, which I would expect to be plenty. Any help appreciated!

import numpy as np
import scipy as sp

nx = 20000
xs = np.linspace(0,2*np.pi,nx)
dx = 2*np.pi/nx

# Solving y'' = v = sin(x)
v = np.sin(xs)
sol_actual = -np.sin(xs)

# Build second-order finite difference matrix with periodic BCs
diagonals = [1,-2,1]
A = sp.sparse.diags_array(diagonals,offsets=[-1,0,1],shape=(nx,nx))
a_l = A.tolil()
a_l[-1,0] = 1
a_l[0,-1] = 1
a_l = a_l/dx**2
A_sparse = a_l.tocsr()
A_numpy = A_sparse.toarray()

sol_sparse = sp.sparse.linalg.spsolve(A_sparse,v)
sol_np = np.linalg.solve(A_numpy,v)


from matplotlib import pyplot as plt
plt.plot(xs,v,label='d2x/dx2')
plt.plot(xs,sol_actual,label='analytical')
plt.plot(xs,sol_sparse,label='scipy sparse')
plt.plot(xs,sol_np,label='numpy')
plt.legend()

Answer 1

You haven't specified any limit condition, as far as I can see.

All 3 solutions are valid solution of y''=sin(x).
-sin(x),
-sin(x)+0.02 (ca)
or -sin(x)+0.6...
You are lucky not to have a -sin(x)+12x-36 [edit. Not really. Order of your matrix is n-1 not n-2. Even tho you have two meanlingless lines, only one of them is redundant. The other is not. Just meanlingless. But it happens to prevent [1,2,3,4,...,n] to be a solution in addition to [1,1,1,1,...,1]. So you were lucky in the sense that what you did wasn't intended to prevent it. But no numerical solver could have returned -sin(x)+12x+36 from this. Only -sin(x)+36 and the like]

Plus, you have a strange loop back in the first and last row. Which is due to your two extra lines. That, I suspect, were precisely intended to add limit condition. But what it does is "close the loop" (like if you were trying to solve the problem of a vibrating cord under a sin(x) excitation, but that vibrating cord is closed (the two ends are connected together, and the cord is under tension by magic)

Maybe you intended to make the first and last line the limit condition line. As people usually do. Any way, they can't be a 2nd derivative line, since those -2 1 0 0 ... 0 0 1 and 1 0 0 ... 0 0 -2 1 lines are not the scheme of a second derivative. Unless, again, you want to "close the loop". But in which case you matrix is singular. Since no random matrix is really singular (and numerical matrix contain random numerical error, so, no numerical matrix is really singular), you have no complain and a solution any way. But that randomness shows, and match the "free degree" of your solutions. Solutions to

[-2  1  0  0  0  0  1 ]
[ 1 -2  1  0  0  0  0 ]
[ 0  1 -2  1  0  0  0 ]
[ 0  0  1 -2  1  0  0 ]
[ 0  0  0  1 -2  1  0 ]
[ 0  0  0  0  1 -2  1 ]
[ 1  0  0  0  0  1 -2 ]

are more than one. Because kernel of that is not {0}. So you can add any combination of the kernel to any solution and still have a solution. Here , [1 1 1 1 1 1 1] is obviously a solution.

Numerically, [1 1 1 1 1 1 1] is only "almost 0". And the randomness of what happens when numpy or scipy is dividing what ever it got by the determinant give the offset.

Another (probably more accurate) way to say it, is that those two lines are redundant. The first is the same as minus the sum of the others. So you can remove it

Which means then that you have n-1 equations for n unknowns. So a 1D variety of solutions. As I said, whatever solution plus any amount as you want of [1,1,1,1,1,1,1], since sum of all columns is 0.

Which is exactly what you see: each method adds it own amount of [1,1,1,1,1,1,1] (that is it's own offset: y=-sin(x)+C)

So, long story short, your matrix is singular, and does not specify limit condition.

What you probably meant to do, starting from the moment you got the initial sparse matrix, is

a_l = A.tolil()
a_l[0,0] = 1
a_l[0,1] = 0
a_l[-1,-1] = 1
a_l[-1, -2] =0

Which then gives this kind of matrix

array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])

(before division by dx² of course)

Which is a legit way to solve y''=sin(x). The inner lines are a derivative scheme for y'' (where it can exist. That is not 1st and last line). And the first and last lines are limit condition. Here, the force y[0] and y[-1] to be v[0] and v[-1], where you should put the limit values. But here, they happen to be already as you obviously want (again, you said nothing about limit. v[0]=15 and v[-1]=-5 would be valid too. But you implicitly said, by saying that any other solution than y=-sin(x) is a bug, that you wanted specifically the solution for which y[0]=y[-1]=0)

All togerther

import numpy as np
import scipy as sp

nx = 2000 # on my machine 20000 crash the interpreter. And any way, I suspect you have increased that nx first expecting it to get "more accurate" result. But since your machine is better than mine (you didn't mention any crash) you can put this back to 20000

xs = np.linspace(0,2*np.pi,nx)
dx = 2*np.pi/nx

# Solving y'' = v = sin(x)
v = np.sin(xs)
sol_actual = -np.sin(xs)

# Build second-order finite difference matrix with periodic BCs
diagonals = [1,-2,1]
A = sp.sparse.diags_array(diagonals,offsets=[-1,0,1],shape=(nx,nx))
a_l = A.tolil()

a_l[0,0] = 1
a_l[0,1] = 0
a_l[-1,-1] = 1
a_l[-1, -2] =0

a_l = a_l/dx**2
A_sparse = a_l.tocsr()
A_numpy = A_sparse.toarray()

sol_sparse = sp.sparse.linalg.spsolve(A_sparse,v)
sol_np = np.linalg.solve(A_numpy,v)


from matplotlib import pyplot as plt
plt.plot(xs,v,label='d2x/dx2')
plt.plot(xs,sol_actual,label='analytical')
plt.plot(xs,sol_sparse,label='scipy sparse')
plt.plot(xs,sol_np,label='numpy')
plt.legend()

Which gives exactly what you wanted, I surmise