Discrepancy between analytic solution and solution by relaxation method

Question

So I am trying to solve the differential equation $\frac{d^2y}{dx^2} = -y(x)$ subject to boundary conditions y(0) = 0 and y(1) = 1 ,the analytic solution is y(x) = sin(x)/sin(1).

I am using three point stencil to approximate the double derivative. The curves obtained through these ways should match at least at the boundaries ,but my solutions have small differences even at the boundaries. I am attaching the code, Please tell me what is wrong.

import numpy as np
import scipy.linalg as lg
from scipy.sparse.linalg import eigs
from scipy.sparse.linalg import inv
from scipy import sparse
import matplotlib.pyplot as plt

a = 0
b = 1
N = 1000
h = (b-a)/N
r = np.arange(a,b+h,h)

y_a = 0
y_b = 1

def lap_three(r):
    h = r[1]-r[0]
    n = len(r)
   
    
    M_d = -2*np.ones(n)
   #M_d = M_d + B_d
    O_d = np.ones(n-1)
    mat = sparse.diags([M_d,O_d,O_d],offsets=(0,+1,-1))
   #print(mat)
    return mat


def f(r):
    h = r[1]-r[0]
    n = len(r)
    return -1*np.ones(len(r))*(h**2)

def R_mat(f,r):
    r_d = f(r)
    R_mat = sparse.diags([r_d],offsets=[0])
   #print(R_mat)
    return R_mat

#def R_mat(r):
 #   M_d = -1*np.ones(len(r))
    
    
def make_mat(r):
    main = lap_three(r) - R_mat(f,r)
    return main

main = make_mat(r)
main_mat = main.toarray()
print(main_mat)
'''
eig_val , eig_vec = eigs(main, k = 20,which = 'SM')

#print(eig_val)
Val = eig_vec.T
plt.plot(r,Val[0])
'''
main_inv = inv(main)
inv_mat = main_inv.toarray()
#print(inv_mat)
#print(np.dot(main_mat,inv_mat))
n = len(r)
B_d = np.zeros(n)
B_d[0] = 0
B_d[-1] = 1
#print(B_d)
#from scipy.sparse.linalg import spsolve
A = np.abs(np.dot(inv_mat,B_d))
plt.plot(r[0:10],A[0:10],label='calculated solution')

real = np.sin(r)/np.sin(1)

plt.plot(r[0:10],real[0:10],label='analytic solution')
plt.legend()
#plt.plot(r,real)
#plt.plot(r,A)

'''diff = A-real
plt.plot(r,diff)'''

Answer 1

There is no guarantee of what the last point in arange(a,b+h,h) will be, it will mostly be b, but could in some cases also be b+h. Better use

r,h = np.linspace(a,b,N+1,retstep=True)

The linear system consists of the equations for the middle positions r[1],...,r[N-1]. These are N-1 equations, thus your matrix size is one too large.

You could keep the matrix construction shorter by including the h^2 term already in M_d.

If you use sparse matrices, you can also use the sparse solver A = spsolve(main, B_d).

The equations that make up the system are

A[k-1] + (-2+h^2)*A[k] + A[k+1] = 0

The vector on the right side thus needs to contain the values -A[0] and -A[N]. This should clear up the sign problem, no need to cheat with the absolute value.

The solution vector A corresponds, as constructed from the start, to r[1:-1]. As there are no values for postitions 0 and N inside, there can also be no difference.

---

PS: There is no relaxation involved here, foremost because this is no iterative method. Perhaps you meant a finite difference method.