6
$\begingroup$

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:

approximate vs. exact

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?

$\endgroup$
  • $\begingroup$ You are correct that you have a Robin condition at $x = 1$. $\endgroup$ – Geoff Oxberry Oct 14 '14 at 8:16
  • $\begingroup$ @GeoffOxberry Thanks for confirming that. $\endgroup$ – user89 Oct 14 '14 at 15:30
5
$\begingroup$

This is actually a very frustrating kind of bug. Change the line

I = 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 is

array([[ 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.

$\endgroup$
  • $\begingroup$ Thank you so much! This was exactly it. Checking the residual is a smart idea. $\endgroup$ – user89 Oct 14 '14 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.