Exercise 1: Solitary wave

Figure 1: Solitary wave.

This exercise focuses on a solitary wave propagating on water. (PDF-version)

The differential equation for a solitary wave can be written as: $$ \begin{align} \frac{d^2 Y}{d X^2} = \frac{3Y}{D^2}\left(\frac{A}{D}-\frac{3}{2D}Y\right) \label{eq:wave} \end{align} $$

where $ D $ is the middle depth, $ Y(X) $ is the wave height above middle depth and A is the wave height at $ X=0 $. The wave is symmetric with respect to $ X=0 $. See Figure 1. The coordinate system follows the wave.

By using dimensionless variables: $ x=\frac{X}{D} $, $ a=\frac{A}{D} $, $ y=\frac{Y}{A} $, Eq. \eqref{eq:wave} can be written as: $$ \begin{align} y''(x)=a \, 3 \, y(x) \, \left(1 - \frac{3}{2}y(x) \right) \label{eq:wave2} \end{align} $$

initial conditions: $ y(0)=1 $, $ y'(0)=0 $. Use $ a=\frac{2}{3} $

Pen and paper

The following problems should be done using pen and paper:
a) Calculate y(0.6), and y'(0.6) using euler's method with $ \Delta x = 0.2 $
b) Solve a) using Heuns's method.
c) Perform a taylor expansion series around $ x=0 $ (include 3 parts) on Eq. \eqref{eq:wave2}, and compare results with a) and b).
d) Solve Eq. \eqref{eq:wave2} analytically for $ a=\frac{2}{3} $, given: $$ \begin{align} \int \frac{dy}{y \, \sqrt{1-y}} = -2 \, arctanh \left(\sqrt{1-y}\right) \nonumber \end{align} $$

Compare with solutions in a), b) and c).

Programing

Write a program that solve a), b) and c) numerically, and compare with the analytical solution found in d). Solve first with $ \Delta x = 0.2 $, and experiment with different values.

Hints for both Pen and paper and Programming problems:

Hint 1.

Solutions:
a) $ y(0.6)=0.88 $, $ y'(0.6)=-0.569 $.

b) $ y(0.6)=0.8337 $, $ y'(0.6)=-0.4858 $.

c) $ y \approx 1 - \frac{x^2}{2} + \frac{x^4}{6} $, $ y' \approx -x +\frac{2}{3}x^3 $

d $ y = \frac{1}{cosh^2\left(x \, /\sqrt{2}\right)} = \frac{1}{1 + cosh\left(\sqrt{2} \cdot x\right)} $

Hint 2.

Figure 2: Plot should look something like this.

If you want you can use this template and fill in the lines where it's indicated.

# src-ch1/solitaryWave.py
#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()
import numpy as np

#### set default plot values: ####
LNWDT=3; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT

""" This script solves the problem with the solitary wave:

        y'' = a*3*y*(1-y*3/2)
        
        y(0) = 1, y'(0) = 0
        
    or as a system of first order differential equations (y0 = y, y1 = y'):
        
        y0' = y'
        y1' = a*3*y0*(1-y0*3/2)
        
        y0(0) = 1, y1(0) = 0
        
"""
a = 2./3
h = 0.2 # steplength dx
x_0, x_end = 0, 0.6

x = np.arange(x_0, x_end + h, h) # allocate x values

#### solution vectors: ####
Y0_euler = np.zeros_like(x) # array to store y values
Y1_euler = np.zeros_like(x) # array to store y' values

Y0_heun = np.zeros_like(x)
Y1_heun = np.zeros_like(x)

#### initial conditions: ####
Y0_euler[0] = 1 # y(0) = 1
Y1_euler[0] = 0 # y'(0) = 0

Y0_heun[0] = 1 
Y1_heun[0] = 0 


#### solve with euler's method ####

for n in range(len(x) - 1):
    y0_n = Y0_euler[n] # y at this timestep
    y1_n = Y1_euler[n] # y' at this timestep
    
    "Fill in lines below"
    f0 = 
    f1 = 
    "Fill in lines above"
    
    Y0_euler[n + 1] = y0_n + h*f0
    Y1_euler[n + 1] = y1_n + h*f1

#### solve with heun's method: ####

for n in range(len(x) - 1):
    y0_n = Y0_heun[n] # y0 at this timestep (y_n)
    y1_n = Y1_heun[n] # y1 at this timestep (y'_n)
    
    "Fill in lines below"
    f0 = 
    f1 = 
    
    y0_p = 
    y1_p = 
    
    f0_p = 
    f1_p = 
    "Fill in lines above"
    
    Y0_heun[n + 1] = y0_n + 0.5*h*(f0 + f0_p)
    Y1_heun[n + 1] = y1_n + 0.5*h*(f1 + f1_p)
    

Y0_taylor = 1 - x**2/2 + x**4/6
Y1_taylor = -x + (2./3)*x**3

Y0_analytic = 1./(np.cosh(x/np.sqrt(2))**2)


#### Print and plot solutions: ####

print "a) euler's method: y({0})={1}, y'({2})={3}".format(x_end, round(Y0_euler[-1], 4), x_end, round(Y1_euler[-1], 4))
print "b) heun's method: y({0})={1}, y'({2})={3}".format(x_end, round(Y0_heun[-1], 4), x_end, round(Y1_heun[-1], 4))
print "c) Taylor series: y({0})={1}, y'({2})={3}".format(x_end, round(Y0_taylor[-1], 4), x_end, round(Y1_taylor[-1], 4))
print "d) Analytical solution: y({0})={1}".format(x_end, round(Y0_analytic[-1], 4))

plt.figure()
plt.plot(x, Y0_euler, 'r-o')
plt.plot(x, Y0_heun, 'b-^')
plt.plot(x, Y0_taylor, 'g-*')
plt.plot(x, Y0_analytic, 'k--')

eulerLegend = 'euler, y({0})={1}'.format(x_end, round(Y0_euler[-1], 4))
heunLegend = 'heun, y({0})={1}'.format(x_end, round(Y0_heun[-1], 4))
taylorLegend = 'taylor, y({0})={1}'.format(x_end, round(Y0_taylor[-1], 4))
analyticLegend = 'analytic, y({0})={1}'.format(x_end, round(Y0_analytic[-1], 4))

plt.legend([eulerLegend, heunLegend, taylorLegend, analyticLegend], loc='best', frameon=False)
plt.show()