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 division
import numpy as np
import 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 ODE
a = 0
b = 1
num_subintervals = 4
num_points = num_subintervals+1
#size of each subinterval
h = (b - a)/num_subintervals
X = np.linspace(a, b, num_points)
def f(x):
return x
vectorized_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 = 0
F[1] = F[1] - (u_0/(h**2))
#We have eliminated U_0 from our equations due to implementation of BCs
F = F[1:]
#eliminated U_0, but U_n is still in the equation, thus num_points-1 to solve for
N = num_points-1
# -2 along main diagonal, 1 along +1 and -1 diagonal
D2 = 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**2
D2[N-1, N-2] = 2/(h**2) #Modifying U_n-1 coefficient (index N-2) in equation for U''_n
D2[N-1, N-1] = (-2 - 2*h)/(h**2) #Modifying U_n (index N-1) in equation for U''_n
I = np.ones((num_points-1, num_points-1))
#LHS
A = -1*D2 + I
U = np.concatenate(([u_0], np.linalg.solve(A, F)))
error = np.linalg.norm(U - u, ord=np.inf)
print "error: ", error
plt.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?
