I am attempting to solve the following ODE problem:$$-u''+ u = x$$$$u(0) = 0$$$$u'(1) = -u(1)$$The exact solution is: $u(x) = e^{-x-1} - e^{x-1} + x$I have a Dirichlet at $x = 0$ and a Robin condition at $x = 1$. I refer to the exact solution with little $u$, and the approximate solution with big $U_i$, where the subscript $i$ refers to point number: if I have $m$ sub-intervals, then I will have $m+1$ points indexed from $0$ to $m$.In order to implement the Dirichlet condition, I eliminated $u_0$ from the finite difference matrix, and modified the right hand side, but in this case I didn't really modify the RHS since $u(0) = 0$, thus $U_0 = 0$.$$u''(0+h) \approx \frac{U_0 - 2U_1 + U_2}{h^2} = \frac{- 2U_1 + U_2}{h^2}$$In order to implement the Robin condition, I decided to use the ghost point method along with a centred difference approximation of the first derivative ($U_m = u(1)$, if there are $m$ points).$$u'(1) = -u(1),\;\therefore u'(1) \approx \frac{U_{m+1} - U_{m-1}}{2h} = -U_m$$$$U_{m+1} = U_{m-1} - 2hU_{m}$$$$\therefore u''(1) \approx \frac{U_{m-1} - 2U_m + U_{m+1}}{h^2} = \frac{2U_{m-1} + (-2 - 2h)U_m}{h^2}$$So, if I have 4 sub-intervals (5 points), where $U_0$ and $U_4$ are the correspond to $u(0)$ and $u(1)$ respectively, I should have the following matrix equation:$$\left(\frac{-1}{h^2}\begin{bmatrix}-2 & 1 & 0 & 0\\ 1 & -2 & 1 & 0\\ 0 & 1 & -2 & 1\\ 0 & 0 & 2 & -2 - 2h \end{bmatrix} + \mathbf{I}\right)\begin{bmatrix}U_1 \\ U_2 \\ U_3 \\ U_4\end{bmatrix} = \begin{bmatrix}F_1 \\ F_2 \\ F_3 \\ F_4\end{bmatrix}$$I then implemented the above approximation using some simple Python code (you should be able to run it too), and I have included comments liberally to aid in readability:from __future__ import divisionimport numpy as npimport matplotlib.pyplot as plt#Want to solve ODE: -u'' + u = x, u(0) = 0, u'(1) = -u(1)#Actual solution: exp(-X-1)# a and b mark the interval over which we want to solve our ODEa = 0b = 1num_subintervals = 4num_points = num_subintervals+1#size of each subintervalh = (b - a)/num_subintervalsX = np.linspace(a, b, num_points)def f(x):    return xvectorized_f = np.vectorize(f)# actual solution: -u'' + u = x, u(0) = 0, u'(1) = -u(1)u = np.exp(-X-1) - np.exp(X-1) + X#Implementing Dirichlet condition at u(0)F = vectorized_f(X)    u_0 = 0F[1] = F[1] - (u_0/(h**2))#We have eliminated U_0 from our equations due to implementation of BCsF = F[1:]#eliminated U_0, but U_n is still in the equation, thus num_points-1 to solve forN = num_points-1# -2 along main diagonal, 1 along +1 and -1 diagonalD2 = np.diag((-2/(h**2))*np.ones(N), 0) + np.diag((1/(h**2))*np.ones(N-1), 1) + np.diag((1/(h**2))*np.ones(N-1), -1)#D1(U_n) = (U_n+1 - U_n-1)/2h = -U_1# -2h*U_n + U_n-1 = U_n+1# D2(U_n) = (U_n-1 - 2*U_n + U_n+1)/h**2 = (U_n-1 - 2*U_n - 2*h*U_n + U_n-1)/h**2# D2(U_n) = (-2*U_n - 2*h*U_n + 2*U_n-1)/h**2D2[N-1, N-2] = 2/(h**2) #Modifying U_n-1 coefficient (index N-2) in equation for U''_nD2[N-1, N-1] = (-2 - 2*h)/(h**2) #Modifying U_n (index N-1) in equation for U''_nI = np.ones((num_points-1, num_points-1))#LHSA = -1*D2 + IU = np.concatenate(([u_0], np.linalg.solve(A, F)))error = np.linalg.norm(U - u, ord=np.inf)print error: , errorplt.plot(X, U, label=approximate)plt.plot(X, u, label=exact)plt.legend(loc=best)plt.show()I get the following image: And the error does not decrease as I increase number of sub-intervals. So, clearly, I have an error, but I am not sure where. Could you help?
Error implementing Robin boundary conditions in toy ODE problem
finite difference;homework
This is actually a very frustrating kind of bug. Change the lineI = np.ones((num_points-1, num_points-1))to I = np.eye(num_points-1)If you look at your matrix $A=-D_2+I$, it isarray([[ 129.,  -63.,    1.,    1.,    1.,    1.,    1.,    1.],       [ -63.,  129.,  -63.,    1.,    1.,    1.,    1.,    1.],       [   1.,  -63.,  129.,  -63.,    1.,    1.,    1.,    1.],       [   1.,    1.,  -63.,  129.,  -63.,    1.,    1.,    1.],       [   1.,    1.,    1.,  -63.,  129.,  -63.,    1.,    1.],       [   1.,    1.,    1.,    1.,  -63.,  129.,  -63.,    1.],       [   1.,    1.,    1.,    1.,    1.,  -63.,  129.,  -63.],       [   1.,    1.,    1.,    1.,    1.,    1., -127.,  145.]])which is not what you want. You can catch these kinds of bugs by computing $AU-F$ and $Au-F$ and checking the residual.