12 Week 12 Participation

Name: [Type Your Name Here]

To begin all assignments (whether participation or homework), please save a copy of this notebook to your Google Drive by clicking File -> Save a copy in Drive

import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from ipywidgets import interactive

12.1 2D Heat Equation - Numerical Solution

After writing out your steps, replace the ??? with the appropriate indices.

\[ U_{i,j}^{n+1} = U_{???,???}^{???} + \frac{D \cdot ???}{???} \left( ??? + ??? - ??? + ??? + ??? \right) \]

Practice

Now, we need to implement the finite difference scheme that you developed in the previous problem. As a model problem, consider the 2D heat equation

\[ u_t = D(u_{xx} + u_{yy}) \]

on the domain \((x,y) \in [0,1] \times [0,1]\) with the initial condition \(u(0,x,y) = \sin(\pi x)\sin(\pi y)\), homogeneous Dirichlet boundary conditions, and
\(D=1\). Fill in the empty spaces in the following code chunks.

First we import the proper libraries and set up the domains for \(x, y,\) and \(t\).

np.linspace(0,1,11)

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm # this allows for color maps
from ipywidgets import interactive

# Write code to build a linearly spaced array of x values
# starting at 0 and ending at exactly 1
x = np.linspace(0,1,26)# your code here
y = x # This is a step that allows for us to have y = x
# The consequence of the previous line is that dy = dx.
#dx = 1/100 # Extract dx from your array of x values. #
dx = x[1] - x[0]
# Now write code to build a linearly spaced array of time values
# starting at 0 and ending at 0.25.
# You will want to use many more values for time than for space
# (think about the stability conditions from the 1D heat equation).
t = np.linspace(0,0.25,1001) # your code here
dt = t[1] - t[0] #Extract dt from your array of t values

# Next we will use the np.meshgrid() command to turn the arrays of
# x and y values into 2D grids of x and y values.
# If you match the corresponding entries of X and Y then you get
# every ordered pair in the domain.
X, Y = np.meshgrid(x,y)

# Next we set up a 3 dimensional array of zeros to store all of
# the time steps of the solutions.
U = np.zeros( (len(t), len(x), len(y)))

U[0,:,:] = np.sin(np.pi*x) * np.sin(np.pi*y) # initial condition depending on X and Y
U[:,0,:] = 0# boundary condition for x=0
U[:,-1,:] = 0# boundary condition for x=1 (NOTE: -1 represents the last element of the array)
U[:,:,0] = 0# boundary condition for y=0
U[:,:,-1] = 0# boundary condition for y=1

D = 1
a = D*dt/dx**2
print(a)

0.15625

for n in range(len(t)-1):
  U[n+1,1:-1,1:-1] = U[n,1:-1,1:-1] + \
    a*(U[n, 2:, 1:-1] + \
       U[n, :-2, 1:-1] - \
       4*U[n, 1:-1, 1:-1] + \
       U[n, 1:-1, 2:] + \
       U[n, 1:-1, :-2])

def plotter(Frame):
  fig = plt.figure(figsize=(12,10))
  ax = plt.axes(projection='3d')
  ax.plot_surface(X,Y,U[Frame,:,:], cmap=cm.coolwarm)
  ax.set_zlim(0,1)
  plt.show()

interactive_plot = interactive(plotter, Frame=(0,len(t)-1))
interactive_plot

fig = go.Figure(data=go.Heatmap(
                   z=U[100,:,:],
                   x=x,
                   y=y,
                   hoverongaps = False))
fig.show()

import plotly.graph_objects as go
import pandas as pd
import numpy as np
TEMP = 100

# Create figure

z = U[TEMP,:,:]
fig = go.Figure(data=[go.Surface(z=z, x=x, y=y)])
fig.update_layout(title='Heat Equation Solution', autosize=False,
                  width=500, height=500,
                  margin=dict(l=65, r=50, b=65, t=90))
fig.show()

# this is not working :( --- UPDATE, it's working now :]
fig = go.Figure(
    data=[go.Surface(x=[], y=[], z=[],cmin=0, cmax=1)])


fig.update_layout(

         scene = dict(

        xaxis=dict(range=[min(x), max(x)], autorange=False),
        yaxis=dict(range=[min(y), max(y)], autorange=False),
        zaxis=dict(range=[0, 1], autorange=False),
        )),

frames = [go.Frame(data= [go.Surface(
                                       x=x,
                                       y=y,
                                       z=U[k+1,:,:])],

                   traces= [0],
                   name=f'frame{k}'
                  )for k  in  range(len(U[:,:,:])-1)]

fig.update(frames=frames),


fig.update_layout(updatemenus=[dict(type="buttons",
                          buttons=[dict(label="Play",
                                        method="animate",
                                        args=[None, dict(frame=dict(redraw=True,fromcurrent=True, mode='immediate'))      ])])])


fig.layout.updatemenus[0].buttons[0].args[1]['frame']['duration'] = 1

fig.show()

# this is not working :( --- UPDATE, it's working now :]
fig = go.Figure(
    data=[go.Heatmap(x=[], y=[], z=[], hoverongaps = False, zmin=0,zmax=1)])

fig.update_layout(

         scene = dict(

        xaxis=dict(range=[0, 1], autorange=False),
        yaxis=dict(range=[0, 1], autorange=False),
        zaxis=dict(range=[0, 1], autorange=False),
        )),

frames = [go.Frame(data= [go.Heatmap(
                                       x=x,
                                       y=y,
                                       z=U[k+1,:,:])],
                   traces= [0],
                   name=f'frame{k}'
                  )for k  in  range(len(U[:,:,:])-1)]

fig.update(frames=frames),


fig.update_layout(updatemenus=[dict(type="buttons",
                          buttons=[dict(label="Play",
                                        method="animate",
                                        args=[None, dict(frame=dict(redraw=True,fromcurrent=True, mode='immediate'))      ])])])


fig.layout.updatemenus[0].buttons[0].args[1]['frame']['duration'] = 10

fig.show()

Theorem. In order for the finite difference solution to the 2D heat equation on a square domain to be stable then we need

\[ D\Delta t / \Delta x^2 < ??? \]

Experiment with several parameters to imperically determine the bound.

12.2 The Wave Equation

The problems that we’ve dealt with thus far all model natural diffusion processes: heat transport, molecular diffusion, etc. Another interesting physical phenomenon is that of wave propagation. Previously it was given that the 1D wave equation is

\[ u_{tt} = c \; u_{xx} \]

where \(c\) is a parameter modeling the stiffness of the medium the wave is traveling through.

Let’s take a few minutes to derive this in the guitar string example.

Practice

As before, we use the notation \(U^n_i\) to represent the approximate solution \(u(t,x)\) at the point \(t=t_n\) and \(x=x_i\).

Based on our previous work, give an approximation of \(u_{tt}\)

\[ u_{tt}(t_n, x_i) \approx \frac{??? - ??? + ???}{???}. \]

Give an approximation of \(u_{xx}\)

\[ u_{xx}(t_n, x_i) \approx \frac{??? - ??? + ???}{???}. \]

Put your answers from parts (a) and (b) together using the 1D wave equation

\[ \frac{??? - ??? + ???}{\Delta t^2} = c \left( \frac{??? - ??? + ???}{\Delta x^2} \right). \]

Be sure that your indexing is correct: the superscript \(n\) is the index for time and the subscript \(i\) is the index for space.

ANS:

Rearrange your result from part (c) to solve for \(U_i^{n+1}\)

\[ U_i^{n+1} = \]

You should notice that the finite difference scheme for the wave equation references two different times:

\[ U^n_i \text{ and } U^{n-1}_i. \]

Based on this observation, what information do we need to in order to actually start our numerical solution?

ANS:

Practice

Now we want to implement your answers of the previous exercise to approximate the solution to the following problem:

\[ \text{Solve: } u_{tt} = 2 u_{xx} \]

with

\[ x \in (0,1), \, u(0,x) = 4x (1-x), \, u_t(0,x) = 0, \, \text{and} \, u(t,0) = u(t,1) = 0. \]

12.3 The Wave Equation - 2D

Now consider the 2D wave equation

\[ u_{tt} = c\left( u_{xx} + u_{yy} \right). \]

We want to build a numerical solution to this new PDE. Just like with the 2D heat equation we propose the notation \(U_{i,j}^n\) for the approximation of the function \(u(t,x,y)\) at the point \(t=t_n, x=x_i, y=y_i\).

Practice

Give discretizations of the derivatives \(u_{tt}, u_{xx},\) and \(u_{yy}\)

ANS:

Substitute your discretizations into the 2D wave equation, make the simplifying assumption that \(\Delta x = \Delta y\), and solve for \(U^{n+1}_{i,j}\). This is the finite difference scheme for the 2D wave equation.

ANS:

Write code to implement the finite difference scheme from part (b) on the domain \((x,y) \in (0,1)\times (0,1)\) with homogeneous Dirichlet boundary conditions, initial condition \(u(0,x,y) = \sin(2\pi(x-0.5))\sin(2\pi(y-0.5))\) and zero initial velocity.

ANS:

Draw the finite difference stencil for the 2D WAVE equation.

ANS: