538 Day 14

imports

In [6]:
import numpy as np
import sympy as sp
sp.init_printing()
from matplotlib import rcdefaults
rcdefaults()  
%config InlineBackend.figure_format='retina'   
import seaborn as sns
import matplotlib.pyplot as plt
sns.set()
#%matplotlib notebook
In [7]:
x0 = 0; x1 = 1
alpha = 1; beta = 10
def p(x): return 50*x
def q(x): return -30*x**2
def r(x): return 20*np.sin(x)

plt.figure(figsize=(10,10))

# by 2nd order finite differences
for N in [200]:
    x = np.linspace(x0,x1,N+2)
    h = x[1]-x[0]
    a = -1 + h**2/2*q(x)
    b = 1/2*(1-h/2*p(x))
    c = 1/2*(1+h/2*p(x))
    A = np.zeros((N,N))
    i = np.arange(N)  # a range of indices
    A[i,i] = a[1:-1]            # main diagonal
    A[i[1:],i[:-1]] = b[2:N+1]  # subdiagonal
    A[i[:-1],i[1:]] = c[1:N]    # superdiagonal
    R = h**2/2*r(x)
    R[1] -= b[1]*alpha
    R[N] -= c[N]*beta
    y = np.empty(N+2)
    y[ 0] = alpha
    y[-1] = beta
    y[1:-1] = np.linalg.solve(A,R[1:-1])
    plt.plot(x,y,'.',label='FD N='+str(N))