<jupyter_start><jupyter_text>Transient Partial Differential Equationwriten by <NAME>(<EMAIL>) IntroductionPartial differential equations ([PDEs](https://en.wikipedia.org/wiki/Partial_differential_equation))are equations that involve rates of change with respect to continuous variables. Classic domains where PDEs are used include acoustics, fluid dynamics, electrodynamics, and heat transfer.Now, we generalize to the instationary heat equation $T_t = T_{xx} + T_{yy}$ on the unit square $]0; 1[^2$ with the temperature T, the two-dimensional coordinates x and y, the time t, homogeneous Dirichlet boundary conditions $T(x, y, t) = 0$  forall $(x, y)$ in $ ∂]0; 1[^2$; $t$ in $]0; ∞[$and$T(x, y, 0) = 1.0$ forall $(x, y)$ in $]0; 1[^2$as initial condition.For the spatial discretization, we again use the grid points$ \{(x_i, y_j ) = (i ·\frac{1}{h_x}, j ·\frac{1}{h_y}) 	\forall i = 0, 1, . . . , N_x, N_x + 1; j = 0, 1, . . . , Ny, Ny + 1\}$and the finite difference approximation of the second derivatives$T_{xx|i,j}\approx \frac{T_{i−1,j} − 2T_{i,j} + T_{i+1,j}}{h_x^2}$and$T_{yy|i,j}\approx \frac{T_{i,j-1} − 2T_{i,j} + T_{i,j+1}}{h_y^2}$for $i = 1, . . . , N_x$, $j = 1, . . . , N_y$ with $h_x =\frac{1}{N_x+1}$ and $h_y =\frac{1}{N_y+1}$<jupyter_code>%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Boundary of the PDE 
left=0
right= 1

down=0
up= 1
N =[3,7,31] # number of division for x
dt_all=[1/64, 1/128, 1/2048, 1/4096]# dt 

#initial Condition
T=np.zeros((N[-1]+2,N[-1]+2))
T[1:N[-1]+1,1:N[-1]+1]=1
#T
fig = plt.figure()
ax = Axes3D(fig)
dx = (right-left) / (N[-1]+2)
x = np.arange(left, right, dx)
y = np.arange(down, up, dx)
X, Y = np.meshgrid(x, y)


plt.title('Problem at time, t=0')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Temperature')

surf=ax.plot_surface(X,Y,T,rstride=1, cstride=1, antialiased=True)

plt.show()<jupyter_output><empty_output><jupyter_text>Von Nuemann Stability AnalysisIn numerical analysis, [von Neumann stability analysis](https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis) is a procedure used to check the stability of finite difference schemes as applied to linear partial differential equations.<jupyter_code>for Nx in N:
    for dt in dt_all:
        dx = (right-left) / (Nx+1)
    
        if (dt > 0.5*dx*dx):
            print("STABLE  : Explicit Method for Nx="+str(Nx)+" & dt="+str(dt))
        else:
            print("UNSTABLE: Explicit Method for Nx="+str(Nx)+" & dt="+str(dt))
#########################Variables to Store Values only################
def SetupValues():
    explicit={}
    implicit={}  
    for Nx in N:
        explicit[Nx]={}
        implicit[Nx]={}
        for dt in dt_all:
            explicit[Nx][dt]=list()
            implicit[Nx][dt]=list()
    return explicit, implicit
#####################################################################<jupyter_output><empty_output><jupyter_text>Explicit Time SteppingThis function is used to solve the finite difference discreditised matrix using an explicit euler formulation.<jupyter_code>def explicitEuler(Nx,Tin,dt,dx):
    
    tau =dt/(dx*dx)
    Tnew=Tin
    Told=Tin
    explicit_=[]
    
    for i in range(int(0.5/dt)):
        explicit_.append(Told.copy())
        
        for x in range(1,Nx+1):
            for y in range(1,Nx+1):
                Tnew[x,y]= Told[x,y]+  tau*( Told[x-1,y]+Told[x+1,y]+Told[x,y+1]+Told[x,y-1]-4*Told[x,y] )

        Told=Tnew
    return explicit_.copy()

explicit,implicit=SetupValues()
for Nx in N:
    T=np.zeros((Nx+2,Nx+2))
    T[1:Nx+1,1:Nx+1]=1
    dx = (right-left) / (Nx+1)
    
    for dt in dt_all:
        explicit[Nx][dt]=explicitEuler(Nx,T.copy(),dt,dx)
        
#Set Values to be visuaised here
NValue=31
step=0
dtValue=0.015625
#########Visualisation############
figimp = plt.figure()
ax = Axes3D(figimp)
dx = (right-left) / (NValue+2)
x = np.arange(left, right, dx)
y = np.arange(down, up, dx)
X, Y = np.meshgrid(x, y)


plt.title('Problem at timestep='+str(step)+',  for dt='+str(dtValue)+' N='+str(Nx))
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Temperature')
surf=ax.plot_surface(X,Y,explicit[NValue][dtValue][step],rstride=1, cstride=1, antialiased=True)
plt.show()<jupyter_output><empty_output><jupyter_text>Implicit Time Stepping<jupyter_code>def implicitEuler():
    tau =dt/(dx*dx)
    Tnew=Tin
    Told=Tin
    implicit_=[]
    ########### ToDo: Implement the implicit time stepping scheme
    
    #######################################################
    return implicit_.copy()
explicit,implicit=SetupValues()
for Nx in N:
    T=np.zeros((Nx+2,Nx+2))
    T[1:Nx+1,1:Nx+1]=1
    dx = (right-left) / (Nx+1)
    
    for dt in dt_all:
        implicit[Nx][dt]=implicitEuler(Nx,T.copy(),dt,dx)
#Set Values to be visuaised here
NValue=3
step=2
dtValue=0.015625
#########Visualisation############
figimp = plt.figure()
ax = Axes3D(figimp)
dx = (right-left) / (NValue+2)
x = np.arange(left, right, dx)
y = np.arange(down, up, dx)
X, Y = np.meshgrid(x, y)


plt.title('Problem at timestep='+str(step)+',  for dt='+str(dtValue)+' N='+str(Nx))
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Temperature')
surf=ax.plot_surface(X,Y,implicit[NValue][dtValue][step],rstride=1, cstride=1, antialiased=True)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20365535248041775, 'fraction_numerical': 0.025550167847817978, 'mean_word_length': 4.37374749498998, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '215539', 'license': 'MIT', 'n_tokens_mistral': 2186, 'n_tokens_neox': 2118, 'n_words': 472, 'path': 'Notebooks/6.Transient PDEs.ipynb', 'repo_name': 'AbinavRavi/NumericalMethods', 'size': '5384'}