# The Implicit Backward Time Centered Space (BTCS) Difference Equation for the Heat Equation
## The Heat Equation
The Heat Equation is the first order in time ($t$) and second order in space ($x$) Partial Differential Equation: 
$$  \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2},$$
The equation describes heat transfer on a domain
$$ \Omega = \{ t \geq, 0\leq x \leq 1\}. $$
with an initial condition at time $t=0$ for all $x$ and boundary condition on the left ($x=0$) and right side ($x=1$).

## Backward Time Centered Space (BTCS) Difference method
This note book will illustrate the Backward Time Centered Space (BTCS) Difference method for the Heat Equation with the __initial conditions__  
$$ u(x,0)=2x, \ \ 0 \leq x \leq \frac{1}{2}, $$
$$ u(x,0)=2(1-x), \ \ \frac{1}{2}  \leq x \leq 1, $$
and __boundary condition__
$$ u(0,t)=0,  u(1,t)=0. $$



```python
# LIBRARY
# vector manipulation
import numpy as np
# math functions
import math 

# THIS IS FOR PLOTTING
%matplotlib inline
import matplotlib.pyplot as plt # side-stepping mpl backend
import warnings
warnings.filterwarnings("ignore")
from IPython.display import HTML

HTML('''
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')

```





<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>



## Discete Grid
The region $\Omega$ is discretised into a uniform mesh $\Omega_h$. In the space $x$ direction into $N$ steps giving a stepsize of
$$ h=\frac{1-0}{N},$$
resulting in 
$$x[i]=0+ih,$$
and into $N_t$ steps in the time $t$ direction giving a stepsize of 
$$ k=\frac{1-0}{N}$$
resulting in 
$$t[i]=0+ik.$$
The Figure below shows the discrete grid points for $N=10$ and $N_t=15$, the red dots are the unknown values, the green dots are the known boundary conditions and the blue dots are the known initial conditions of the Heat Equation.


```python
N=10
Nt=200
h=1/N
k=1/Nt
r=k/(h*h)
time_steps=15
time=np.arange(0,(time_steps+.5)*k,k)
x=np.arange(0,1.0001,h)
X, Y = np.meshgrid(x, time)
fig = plt.figure()
plt.plot(X,Y,'ro');
plt.plot(x,0*x,'bo',label='Initial Condition');
plt.plot(np.ones(time_steps+1),time,'go',label='Boundary Condition');
plt.plot(x,0*x,'bo');
plt.plot(0*time,time,'go');
plt.xlim((-0.02,1.02))
plt.xlabel('x')
plt.ylabel('time (ms)')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.title(r'Discrete Grid $\Omega_h,$ h= %s, k=%s'%(h,k),fontsize=24,y=1.08)
plt.show();
```

## Discrete Initial and Boundary Conditions

The discrete initial conditions are 
$$ w[i,0]=2x[i], \ \ 0 \leq x[i] \leq \frac{1}{2} $$
$$ w[i,0]=2(1-x[i]), \ \ \frac{1}{2}  \leq x[i] \leq 1 $$
and the discete boundary conditions are 
$$ w[0,j]=0,  w[10,j]=0, $$
where $w[i,j]$ is the numerical approximation of $U(x[i],t[j])$.

The Figure below plots the inital and boundary conditions for $t[0]=0.$


```python
w=np.zeros((N+1,time_steps+1))
b=np.zeros(N-1)
# Initial Condition
for i in range (1,N):
    w[i,0]=2*x[i]
    if x[i]>0.5:
        w[i,0]=2*(1-x[i])
    

# Boundary Condition
for k in range (0,time_steps):
    w[0,k]=0
    w[N,k]=0

fig = plt.figure(figsize=(8,4))
plt.plot(x,w[:,0],'o:',label='Initial Condition')
plt.plot(x[[0,N]],w[[0,N],0],'go',label='Boundary Condition t[0]=0')
#plt.plot(x[N],w[N,0],'go')
plt.xlim([-0.1,1.1])
plt.ylim([-0.1,1.1])
plt.title('Intitial and Boundary Condition',fontsize=24)
plt.xlabel('x')
plt.ylabel('w')
plt.legend(loc='best')
plt.show()
```

## The Implicit Backward Time Centered Space (BTCS) Difference Equation
The implicit Backward Time Centered Space (BTCS) difference equation of the Heat Equation is
\begin{equation*}
\frac{w_{ij+1}-w_{ij}}{k}=\frac{w_{i+1j+1}-2w_{ij+1}+w_{i-1j+1}}{h^2}
\end{equation*}
Rearranging the equation we get
\begin{equation}
-rw_{i-1j+1}+(1+2r)w_{ij+1}-rw_{i+1j+1}=w_{ij}
\end{equation}
for $i=1,...9$ where $r=\frac{k}{h^2}$.

This gives the formula for the unknown term $w_{ij+1}$ at the $(ij+1)$ mesh points
in terms of $x[i]$ along the jth time row.

Hence we can calculate the unknown pivotal values of $w$ along the first row of $j=1$ in terms of the known boundary conditions.

This can be written in matrix form 
$$ A\mathbf{w}_{j+1}=\mathbf{w}_{j} +\mathbf{b}_{j+1} $$
for which $A$ is a $9\times9$ matrix:
$$
\left(\begin{array}{cccc cccc}
1+2r&-r& 0&0&0 &0&0&0\\
-r&1+2r&-r&0&0&0 &0&0&0\\
0&-r&1+2r &-r&0&0& 0&0&0\\
0&0&-r&1+2r &-r&0&0& 0&0\\
0&0&0&-r&1+2r &-r&0&0& 0\\
0&0&0&0&-r&1+2r &-r&0&0\\
0&0&0&0&0&-r&1+2r &-r&0\\
0&0&0&0&0&0&-r&1+2r&-r\\
0&0&0&0&0&0&0&-r&1+2r\\
\end{array}\right)
\left(\begin{array}{c}
w_{1j+1}\\
w_{2j+1}\\
w_{3j+1}\\
w_{4j+1}\\
w_{5j+1}\\
w_{6j+1}\\
w_{7j+1}\\
w_{8j+1}\\
w_{9j+1}\\
\end{array}\right)=
\left(\begin{array}{c}
w_{1j}\\
w_{2j}\\
w_{3j}\\
w_{4j}\\
w_{5j}\\
w_{6j}\\
w_{7j}\\
w_{8j}\\
w_{9j}\\
\end{array}\right)+
\left(\begin{array}{c}
rw_{0j+1}\\
0\\
0\\
0\\
0\\
0\\
0\\
0\\
rw_{10j+1}\\
\end{array}\right).
$$
It is assumed that the boundary values $w_{0j}$ and $w_{10j}$ are known for $j=1,2,...$, and $w_{i0}$ for $i=0,...,10$ is the initial condition.

The Figure below shows the matrix $A$ and its inverse $A^{-1}$ in color plot form for $r$.


```python
A=np.zeros((N-1,N-1))
for i in range (0,N-1):
    A[i,i]=1+2*r

for i in range (0,N-2):           
    A[i+1,i]=-r
    A[i,i+1]=-r
    
Ainv=np.linalg.inv(A)   
fig = plt.figure(figsize=(12,4));
plt.subplot(121)
plt.imshow(A,interpolation='none');
plt.xticks(np.arange(N-1), np.arange(1,N-0.9,1));
plt.yticks(np.arange(N-1), np.arange(1,N-0.9,1));
clb=plt.colorbar();
clb.set_label('Matrix elements values');
clb.set_clim((-1,1));
plt.title('Matrix A r=%s'%(np.round(r,3)),fontsize=24)

plt.subplot(122)
plt.imshow(Ainv,interpolation='none');
plt.xticks(np.arange(N-1), np.arange(1,N-0.9,1));
plt.yticks(np.arange(N-1), np.arange(1,N-0.9,1));
clb=plt.colorbar();
clb.set_label('Matrix elements values');
clb.set_clim((-1,1));
plt.title(r'Matrix $A^{-1}$ r=%s'%(np.round(r,3)),fontsize=24)

fig.tight_layout()
plt.show();
```

## Results
To numerically approximate the solution at $t[1]$ the matrix equation becomes 
$$ \mathbf{w}_{1}=A\mathbf{w}_{0} +\mathbf{b}_{0} $$
where all the right hand side is known. 
To approximate solution at time $t[2]$ we use the matrix equation
$$ \mathbf{w}_{2}=A\mathbf{w}_{1} +\mathbf{b}_{1}. $$
Each set of numerical solutions $w[i,j]$ for all $i$ at the previous time step is used to approximate the solution $w[i,j+1]$. 
The left and right plot below show the numerical approximation $w[i,j]$ of the Heat Equation using the BTCS method for $x[i]$ for $i=0,...,10$ and time steps $t[j]$ for $j=1,...,15$.


```python
fig = plt.figure(figsize=(12,6))

plt.subplot(121)
for j in range (1,time_steps+1):
    w[1:(N),j]=np.dot(Ainv,w[1:(N),j-1])
    plt.plot(x,w[:,j],'o:',label='t[%s]=%s'%(j,time[j]))
plt.xlabel('x')
plt.ylabel('w')
#plt.legend(loc='bottom', bbox_to_anchor=(0.5, -0.1))
plt.legend(bbox_to_anchor=(-.4, 1), loc=2, borderaxespad=0.)

plt.subplot(122)
plt.imshow(w.transpose())
plt.xticks(np.arange(len(x)), x)
plt.yticks(np.arange(len(time)), time)
plt.xlabel('x')
plt.ylabel('time')
clb=plt.colorbar()
clb.set_label('Temperature (w)')
plt.suptitle('Numerical Solution of the  Heat Equation r=%s'%(np.round(r,3)),fontsize=24,y=1.08)
fig.tight_layout()
plt.show()
```


```python

```

81d562c5da2b031fedbc5960caf2115ce9be8a77
ipynb
Jupyter Notebook
Chapter 08 - Heat Equations/.ipynb_checkpoints/03_Heat Equation-Crank_Nicolson-checkpoint.ipynb
jjcrofts77/Numerical-Analysis-Python
97e4b9274397f969810581ff95f4026f361a56a2
MIT
2019-09-05T21:39:12.000Z
2022-03-26T14:00:25.000Z
Chapter 08 - Heat Equations/.ipynb_checkpoints/03_Heat Equation-Crank_Nicolson-checkpoint.ipynb
jjcrofts77/Numerical-Analysis-Python
97e4b9274397f969810581ff95f4026f361a56a2
MIT
Chapter 08 - Heat Equations/.ipynb_checkpoints/03_Heat Equation-Crank_Nicolson-checkpoint.ipynb
jjcrofts77/Numerical-Analysis-Python
97e4b9274397f969810581ff95f4026f361a56a2
MIT
2021-06-17T15:34:04.000Z
2022-01-14T14:53:43.000Z
Qwen/Qwen-72B
1. YES
2. YES