In [1]:
import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
from ipywidgets import interactive
In [2]:
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();
2021-02-02T20:04:30.532310 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/
In [3]:
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)$")
Out[3]:
Text(0, 0.5, '$f(x)$')
2021-02-02T20:04:33.653797 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/
In [4]:
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;
In [5]:
from scipy import sparse
import scipy.sparse.linalg as sl
As = sparse.csc_matrix(A1); 
f = As.dot(b)
In [6]:
%%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
In [7]:
plt.plot(x, f, label="First order BC");
plt.plot(xe, ye, 'ok', label="Exact");
plt.legend(); plt.xlabel("$x$"); plt.ylabel("$f(x)$")
Out[7]:
Text(0, 0.5, '$f(x)$')
2021-02-02T20:04:34.016712 image/svg+xml Matplotlib v3.3.3, https://matplotlib.org/
In [ ]: