###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c) 2014 F.J.Gonzales. Portions of the code adopted from the #numericalmooc materials, also under CC-BY.

# The French Connec... eh solution

## Solving the 1-D wave equation

Welcome to this bonus notebook that ties into the second module of ["Practical Numerical Methods with Python"](http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about).

In the first notebook of the second module, we were introduced to the numerical solutions of *partial differential equations (PDEs)*. The equation used in that notebook was the 1-D linear convection equation, AKA *the advection equation*.

However, in this notebook we will explore the 1-D wave equation shown below:

\begin{equation}\frac{\partial^2 u}{\partial t^2}= c^2 \frac{\partial^2 u}{\partial x^2}.
\end{equation}

The wave equation is second order in both spatial and temporal dimensions.

If the initial displacement $\ f(x)$ and a velocity $\ g(x)$ are specified for the wave equation, then the solution of the system is:
\begin{equation}\ u(x,t) = \frac{1}{2}[f(x-ct)+f(x+ct)]+\frac{1}{2c}\int{g(\tau){\partial \tau}}\end{equation}
This solution is known as the d'Alembert solution for the 1-D wave equation.

For this notebook we will focus on the special case where the initial velocity is zero, $g(x)=0$
and the solution is simplified to the equation below:
\begin{equation}\ u(x,t)=\frac{1}{2}[f(x-ct)+f(x+ct)]\end{equation}


## The Guitar pick example

A classical example of the 1-D wave equation is when simulating the motion of a guitar string after it is plucked.



###### Image credit: G. Everstine, 2012 (Fig. 14, page 23).

As the figure above shows, the original triangular displacement splits into two waves traveling in opposite directions.  One can also see that each of the travelling waves is half the height of the original function, which corresponds to d'Alembert's solution shown previously.

## Wave motion for a bar

Now that we have established the solution for the 1-D wave equation, we will use this to simulate the longitudinal wave motion on a finite length bar. 

### Problem Description

Let's say we have a bar of finite length as shown below. The length is $10 \rm{m}$ in the $x$ direction.  As can be seen, the system ends with a dashpot, also known as a damper; this dashpot is an important part of this problem because we will see how this will create a non-reflecting boundary condition. 



###### Image credit: G. Everstine, 2012 (Fig. 22 in page 30).

In addition to the problem description, the initial displacement $f(x)$ is defined below:
\begin{equation}
u(x,0)=\begin{cases}2-2(x-5)^2 & \text{where } 4\leq x \leq 6,\\
0 & \text{everywhere else in } (0, 10)
\end{cases}
\end{equation}

### Discretizing the equation

Since the 1-D wave equation is second order in time and space, we will discretize the equation using central difference for both dimensions.  Below is the discretized equation:

\begin{equation}\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{c^2\Delta t^2} =  \frac{u_{i+1}^n - 2u_i^n+u_{i-1}^n}{\Delta x^2} \end{equation}

Solving for the unknown in the discretized equation we get the following:

\begin{equation}\ u_i^{n+1}=-u_i^{n-1}+2u_i^n+C^2(u_{i+1}^n-2u_i^{n}+u_{i-1}^n) \end{equation}

In this equation $C= c\frac{\Delta t}{\Delta x}$, which equates to the CFL number.

In addition to discretizing the wave equation, we also want to take a look at our boundary conditions.  On the left-hand-side boundary, the bar is fixed, which would correlate with a Dirichlet boundary where the displacement is $u=0$.

The right-hand-side boundary is more tricky and that is where the idea of the nonreflecting boundary condition comes into play.

If we look at a simplified version of d'Alembert's solution, $f(x-ct)+f(x+ct),$ we can see that the first term represents the wave that is moving towards the right boundary while the latter term represents the wave that bounces away from the right boundary.  In the bar we are analyzing, the dashpot at the right-hand-side means that the second return wave term would go to zero.

Differentiating the first term of D'Alembert's solution in terms of space and time would result in the following PDEs:
\begin{equation}\frac{\partial u}{\partial x}= f^\prime,\frac{\partial u}{\partial t}= -cf^\prime\end{equation}

Rearranging the above equations gives the boundary condition we have to enforce on the right-hand side:

\begin{equation}\frac{\partial u}{\partial x}=-\frac{1}{c}\frac{\partial u}{\partial t}\end{equation}

Discretizing our non-reflecting boundary:

\begin{equation}\frac{u_i^n-u_{i-1}^n}{\Delta x} = \frac{-1}{c}\frac{u_i^{n+1}-u_i^n}{\Delta t}\end{equation}

and isolating for the end boundary:

\begin{equation}u_i^n= -\frac{\Delta x}{c\Delta t}(u_i^{n+1}-u_i^n)+u_{i-1}^n\end{equation}

### Lets rev up some code!

The first thing we always do is to import the necessary libraries for this problem.  In this case we will use bothy Numpy and Sympy. We will also use Matplotlib for plotting our displacement.


```
import numpy as np
import sympy as sp
from sympy.utilities.lambdify import lambdify
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
```

We then create the initial condition for the bar.


```
# initial conditions and parameters
nx= 50
x= sp.symbols('x')
func= 2-2*(x-5)**2
init_func= lambdify((x), func)
x_array= np.linspace(0,10,nx+1)
u_initial= np.asarray([init_func(x0) for x0 in x_array])
```

Lets plot the initial function so we can start off on the right foot.


```
plt.plot(x_array, u_initial, color='#003366', ls='--', lw=3)
plt.xlim(0,10)
plt.ylim(0,3)
plt.xlabel('x')
plt.ylabel('Displacement')
```

Now we can define a function to solve for the displacement.


```
def barsolver(u_init, T, C, nx):
    ''' Returns the displacement of the wave equation for a bar.
    
    Parameters
    ----------
    u_init: array of float
        initial displacement of bar
    T: integer
        final time for calculation
    C: integer
        CFL number for stability
    nx: integer
        number of gridsteps
        
    Returns
    -------
    u: array of float
        final displacement of bar wave motion
    '''    
    #initial parameters
    c= 2
    C2= C**2
    dx= 10./nx
    dt = C*dx/c
    nt= int(round(T/dt))
    

    #create arrays
    u= np.zeros(nx+1)    #array holding u for u[n+1]
    u_1= np.zeros(nx+1)  #array holding u for u[n]
    u_2= np.zeros(nx+1)  #array holding u for u[n-1]
    
    u_1[4/dx:6/dx+1]= u_init[4/dx:6/dx+1]
    
    # Loop for first time step
    n = 0
    for i in range(1, nx):
        u[i] = u_1[i] + 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
    # Enforce boundary conditions    
    u[0] = 0;  u[-1] = -dx/(c*dt)*(u[i]-u_1[i])+u[-2]

    # Switch variables before next step
    u_2[:], u_1[:] = u_1, u
    
    # Loops for subsequent time steps
    for n in range(1, nt):
        for i in range(1, nx):
            u[i] = - u_2[i] + 2*u_1[i] + C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])             

        # Enforce boundary conditions
        u[0] = 0;  u[-1] = -dx/(c*dt)*(u[i]-u_1[i])+u[-2]
    
        # Switch variables before next step
        u_2[:], u_1[:] = u_1, u
        
    return u
```

Lets run the solver for $t=0$ and plot the displacement.


```
u_ans= barsolver(u_initial,0,1,50)
```


```
plt.plot(x_array,u_ans)
plt.xlim(0,10)
plt.ylim(-3,3)
```

As we expected, the displacement is equal to the initial parabola.


```
u_ans= barsolver(u_initial,2,1,50)  #solving at t= 2
```


```
plt.plot(x_array,u_ans)
plt.xlim(0,10)
plt.ylim(-2,2)
```


```
u_ans= barsolver(u_initial,4,1,50)  #solving at t= 4 
```


```
plt.plot(x_array,u_ans)
plt.xlim(0,10)
plt.ylim(-2,2)
```

As we can see from the plots, the wave moving in the right direction has been absorbed while the left moving wave has flipped over and is now moving in the right direction.

Our code has expressed correctly the non-reflecting boundary!

##### Now your turn

What do you think will happen as more time passes?

You will need to plot the displacement for t= 6, 8, 10 and check if what you expected was correct.

## Reference

* Barba, Lorena A., et al. "MAE 6286 Practical Numerical Methods with Python," GW Open edX, The George Washingtion University, 2014. http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about .

* Chapra, Steven C and Canale, Raymond P., "Numerical Methods for Engineers," McGraw-Hill Education, Inc., New York, 2010.

* Everstine, Gordon C. "Analytical Solution of Partial Differential Equations," George Washington University, 2012. [PDF](http://gwu.geverstine.com/pdenum.pdf) 

* Fitzpatrick, Richard. "Introduction: Wave Equation," University of Texas, 2006. [Wave Equation](http://farside.ph.utexas.edu/teaching/329/lectures/node89.html).

* Strauss, Walter A., "Partial Differential Equations: An Introduction," Wiley, Inc., New York, 2008.  

---

###### The cell bellow loads the style of the notebook


```
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
```




<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>
<style>

@font-face {
    font-family: "Computer Modern";
    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

#notebook_panel { /* main background */
    background: rgb(245,245,245);
}

div.cell { /* set cell width */
    width: 750px;
}

div #notebook { /* centre the content */
    background: #fff; /* white background for content */
    width: 1000px;
    margin: auto;
    padding-left: 0em;
}

#notebook li { /* More space between bullet points */
margin-top:0.8em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running { 
    border: 1px solid #111;
}

/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
    background-color: rgb(256,256,256); 
    border-radius: 0px; 
    padding: 0.5em;
    margin-left:1em;
    margin-top: 1em;
}

div.text_cell_render{
    font-family: 'Alegreya Sans' sans-serif;
    line-height: 140%;
    font-size: 125%;
    font-weight: 400;
    width:600px;
    margin-left:auto;
    margin-right:auto;
}


/* Formatting for header cells */
.text_cell_render h1 {
    font-family: 'Nixie One', serif;
    font-style:regular;
    font-weight: 400;    
    font-size: 45pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.5em;
    margin-top: 0.5em;
    display: block;
}	
.text_cell_render h2 {
    font-family: 'Nixie One', serif;
    font-weight: 400;
    font-size: 30pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    display: block;
}	

.text_cell_render h3 {
    font-family: 'Nixie One', serif;
    margin-top:16px;
	font-size: 22pt;
    font-weight: 600;
    margin-bottom: 3px;
    font-style: regular;
    color: rgb(102,102,0);
}

.text_cell_render h4 {    /*Use this for captions*/
    font-family: 'Nixie One', serif;
    font-size: 14pt;
    text-align: center;
    margin-top: 0em;
    margin-bottom: 2em;
    font-style: regular;
}

.text_cell_render h5 {  /*Use this for small titles*/
    font-family: 'Nixie One', sans-serif;
    font-weight: 400;
    font-size: 16pt;
    color: rgb(163,0,0);
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.8em;
    display: block;
}

.text_cell_render h6 { /*use this for copyright note*/
    font-family: 'PT Mono', sans-serif;
    font-weight: 300;
    font-size: 9pt;
    line-height: 100%;
    color: grey;
    margin-bottom: 1px;
    margin-top: 1px;
}

.CodeMirror{
        font-family: "PT Mono";
        font-size: 90%;
}

</style>





a3998d5ac27acb91942d852a4aa7f3ae55e42cf7
ipynb
Jupyter Notebook
Lessons.and.Assignments/1D.Wave.Bar/1D.WaveBar.ipynb
udaypandit/black_scholes
705c3e47e5fd441c30a38c1ab17a80a75441e7d5
MIT
2015-03-30T13:29:15.000Z
2018-06-23T03:39:11.000Z
Lessons.and.Assignments/1D.Wave.Bar/1D.WaveBar.ipynb
udaypandit/black_scholes
705c3e47e5fd441c30a38c1ab17a80a75441e7d5
MIT
2015-01-01T13:55:56.000Z
2015-01-16T19:46:58.000Z
Lessons.and.Assignments/1D.Wave.Bar/1D.WaveBar.ipynb
udaypandit/black_scholes
705c3e47e5fd441c30a38c1ab17a80a75441e7d5
MIT
2015-01-10T08:13:13.000Z
2021-08-10T15:18:26.000Z
Qwen/Qwen-72B
1. YES
2. YES