import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
from ipywidgets import interactive
xe = np.linspace(0,2, 6);
A = 1; B = -np.sinh(2)/np.cosh(2);
ye = A*np.cosh(xe) + B*np.sinh(xe);
plt.plot(xe, ye, label="Exact soln"); plt.legend();
N = 10;
x = np.linspace(0,2,N); dx = x[1] - x[0];
A1 = np.zeros((N, N)); b = np.zeros((N,1));
A2 = np.zeros((N, N));
alpha = -2 - dx**2;
dv = alpha*np.ones((N,));
dv1 = np.ones((N-1,));
A1 = np.diag(dv) + np.diag(dv1,1) + np.diag(dv1,-1);
A1[0,:] = 0; A1[N-1,:] = 0;
A2 = np.diag(dv) + np.diag(dv1,1) + np.diag(dv1,-1);
A2[0,:] = 0; A2[N-1,:] = 0;
A1[0,0] = 1; A2[0,0] = 1;
A1[-1,-2] = -1; A1[-1,-1] = 1;
A2[-1,-1] = alpha; A2[-1,-2] = 2;
b[0] = 1;
f = np.dot(np.linalg.inv(A1), b);
g = np.dot(np.linalg.inv(A2), b);
plt.plot(x, f, label="First order BC");
plt.plot(x, g, label="Second order")
plt.plot(xe, ye, 'ok', label="Exact");
plt.legend(); plt.xlabel("$x$"); plt.ylabel("$f(x)$")
Text(0, 0.5, '$f(x)$')
N = 10;
x = np.linspace(0,2,N); dx = x[1] - x[0];
A1 = np.zeros((N, N)); b = np.zeros((N,1));
alpha = -2 - dx**2;
dv = alpha*np.ones((N,));
dv1 = np.ones((N-1,));
A1 = np.diag(dv) + np.diag(dv1,1) + np.diag(dv1,-1);
A1[0,:] = 0; A1[N-1,:] = 0;
A1[0,0] = 1; A2[0,0] = 1;
A1[-1,-2] = -1; A1[-1,-1] = 1;
b[0] = 1;
from scipy import sparse
import scipy.sparse.linalg as sl
As = sparse.csc_matrix(A1);
f = As.dot(b)
%%time
#f = np.dot(np.linalg.inv(A1), b);
As = sparse.csc_matrix(A1);
f = sl.spsolve(As, b)
print(f)
[1. 0.81117727 0.66241267 0.54635982 0.45728769 0.39079767 0.34360631 0.31338315 0.29863571 0.29863571] Wall time: 2.53 ms
plt.plot(x, f, label="First order BC");
plt.plot(xe, ye, 'ok', label="Exact");
plt.legend(); plt.xlabel("$x$"); plt.ylabel("$f(x)$")
Text(0, 0.5, '$f(x)$')