cuteBoundaryValueProblem.py
#
# Boundary value problem : u''(x) = x - u(x)
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
from scipy.sparse import spdiags
from scipy.sparse import eye
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['toolbar'] = 'None'
plt.rcParams['figure.facecolor'] = 'silver'
# =========================================================================
def solve(a,b,Ua,Ub,n) :
X = linspace(a,b,n)
h = (b-a)/(n-1)
b = array([0,*ones(n-2),0])
A = spdiags([b,-2*b,b],[-1,0,1],n,n)/(h*h)
B = copy(X); B[0] = Ua; B[-1] = Ub
U = spsolve(A.T+eye(n,n),B)
return X,U
# ============================= mainProgram ===============================
a = 0; b = 10; Ua = 1; Ub = 5
a = 0; b = pi; Ua = 1; Ub = 5
a = 0; b = pi; Ua = 0; Ub = pi
plt.figure('u"(x) = x - u(x)')
for i in range(4):
n = 11*2**i
X,U = solve(a,b,Ua,Ub,n)
plt.plot(X,U,'.-b',linewidth=1)
plt.show()
# =========================================================================