%matplotlib inline
%precision 16
import numpy
import matplotlib.pyplot as plt
import scipy.sparse as sparse
import scipy.linalg
import scipy.sparse.linalg as linalg

Before you turn this problem in, make sure everything runs as expected. First, restart the kernel (in the menubar, select Kernel $\rightarrow$ Restart) and then run all cells (in the menubar, select Cell $\rightarrow$ Run All).

Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE", as well as your name and collaborators below:

Enter your name and UNI here (just in case).

NAME = Karim Abdallah

UNI = ka2566

HW 4: Parabolic Problems¶

Question 1 - Separation of Variables¶

Consider the heat equation with zero boundary conditions:

$$ \frac{\partial u}{\partial t} = \kappa \frac{\partial^2 u}{\partial x^2} \qquad \mbox{on } [0,L]\times \mathbb R_+\\ u(x,0) = u_0(x) ~~~~~~ u(0,t) = 0 ~~~~~ u(L,t) = 0. $$

(a) [20] Implement the backward Euler and Crank-Nicholson methods to solve this equation. You are responsible for picking a $\Delta t$ so that the numerical method has enough accuracy but hits the end time $t_f$ exactly (or as much as this is possible with floating point arithmetic). Note that making $\Delta t$ too small may harm your accuracy as much as picking one that is too big.

def solve_heat_BE(m, kappa, L, t_0, t_final, U_0, g_0, g_1):
    """Solve the heat equation using backward Euler and given boundaries
    
    Function will take the necessary number of time steps to reach t_final
    from t_0.  Note that this may mean taking a slightly smaller time step 
    size to hit t_f exactly.
    
    :Input:
     - *m* (int) Number of points use to discretize the domain.  Note that
       the total number of points is *m+1*.
     - *kappa* (float) Diffusion coefficient
     - *L* (float) Length of half of the domain
     - *t_0* (float) Starting time
     - *t_final* (float) Time to integrate to
     - *U_0* (numpy.ndarray) Initial condition at time t_0, should be m+1
     - *g_0* (func) Boundary condition at x=0
     - *g_1* (func) Boundary condition at x=L
     
    :Output:
     - (numpy.ndarray) Solution at time t_final.  Note that this vector should m+1
    """
    
    # Spatial discretization
    # x = numpy.linspace(0.0, L, m+2)
    #delta_x = 2.0/(m-1.0)
    delta_x = L/float(m + 1)
    # Time discretization
    delta_t = 0.41*delta_x
    N = 101
    #t = numpy.arange(0.0, N * delta_t, delta_t)
    
    c = kappa*delta_t/(delta_x**2)
    
    # U
    U = numpy.zeros(m+2)
    U_1 = numpy.zeros(m+2)    
    
    A = numpy.zeros((m+2,m+2))
    b = numpy.zeros(m+2)

    A[0,0] = 1.0
    A[-1,-1] = 1.0

    for i in range(1,m+1):
        A[i,i-1] = -c
        A[i,i] = 1.0+2.0*c
        A[i,i+1] = -c

    
    for i in range(0, m+2):
        U_1[i] = U_0[i]
        
    for n in range(1, N):
        # Compute b and solve linear system
        for i in range(1, m):
            b[i] = -U_1[i]
            
        b[0] = g_0(n-1) + c*U_1[0]
        b[-1] = g_1(n-1) + c*U_1[-1]
        U[:] = scipy.linalg.solve(A, b)

        # Update u_1 before next step
        U_1[:] = U    
    
    
    #raise NotImplementedError()
    
    return U

beta = 150.0
kappa = 0.02
L = 1.0
u_true = lambda x, t: numpy.exp(-(x - 0.4)**2. / (4.0 * kappa * t + 1.0 / beta)) \
                                            / numpy.sqrt(4.0 * kappa * beta * t + 1.0)
g_0 = lambda t: u_true(0.0, t)
g_1 = lambda t: u_true(L, t)

# Discretization and output times
m = 80
x = numpy.linspace(0.0, L, m + 2)
x_fine = numpy.linspace(0.0, L, 100)
delta_x = L / float(m + 1)
output_times = (0.0, 0.5, 1.0, 1.5, 2.0)

# Solve
U_BE = numpy.empty((len(output_times), m + 2))
U_BE[0, :] = u_true(x, 0.0)
for (n, t) in enumerate(output_times[1:]):
    U_BE[n + 1, :] = solve_heat_BE(m, kappa, L, output_times[n], t, U_BE[n, :], g_0, g_1)

error = numpy.linalg.norm(delta_x * (U_BE[-1, :] - u_true(x, output_times[-1])), ord=1)
print "Error BE = %s" % error
#assert error < 1e-3

# Plot some of the results
colors = ['k', 'r', 'b', 'g', 'c']
fig = plt.figure()
axes_BE = fig.add_subplot(1, 1, 1)
for (n, t) in enumerate(output_times):
    axes_BE.plot(x_fine, u_true(x_fine, t), 'r')
    axes_BE.plot(x, U_BE[n, :], "o%s" % colors[n], label='t=%s' % numpy.round(t, 4))
    axes_BE.set_xlabel("x")
    axes_BE.set_ylabel("$U_(x,t)$")
    axes_BE.set_title("Solution to Heat Equation")
    axes_BE.set_xlim([0, L])
    axes_BE.legend()
    
plt.tight_layout()
plt.show()
print t

def solve_heat_CN(m, kappa, L, t_0, t_final, U_0, g_0, g_1):
    """Solve the heat equation using Crank-Nicholson and given boundaries
    
    Function will take the necessary number of time steps to reach t_final
    from t_0.  Note that this may mean taking a slightly smaller time step 
    size to hit t_f exactly.
    
    :Input:
     - *m* (int) Number of points use to discretize the domain.  Note that
       the total number of points is *m+1*.
     - *kappa* (float) Diffusion coefficient
     - *L* (float) Length of half of the domain
     - *t_0* (float) Starting time
     - *t_final* (float) Time to integrate to
     - *U_0* (numpy.ndarray) Initial condition at time t_0, should be m+1
     - *g_0* (func) Boundary condition at x=0
     - *g_1* (func) Boundary condition at x=L
     
    :Output:
     - (numpy.ndarray) Solution at time t_final.  Note that this vector should m+1
    """
    
    # YOUR CODE HERE
    delta_x = L / float(m+1)
    
    C = 10
    
    delta_t = C * delta_x
    
    r = numpy.ones(m+2) * delta_t * kappa / (2.0 * delta_x**2)
    A = sparse.spdiags([-r, 1.0 + 2.0 * r, -r], [-1, 0, 1], m+2, m+2).tocsr()
    
    B = sparse.spdiags([r, 1.0 - 2.0 * r, r], [-1, 0, 1],  m+2, m+2).tocsr()
    
    b = B.dot(U_0)
    b[0] += delta_t * kappa/ (2.0 * delta_x**2) * (g_0(t_0) + g_0(t_final))
    b[-1] += delta_t * kappa / (2.0 * delta_x**2) * (g_1(t_0) + g_1(t_final))
    
        # Solve system
        
    
    U = linalg.spsolve(A, b)
    
    #raise NotImplementedError()
    
    return U

beta = 150.0
kappa = 0.02
L = 1.0
u_true = lambda x, t: numpy.exp(-(x - 0.4)**2. / (4.0 * kappa * t + 1.0 / beta)) \
                                            / numpy.sqrt(4.0 * kappa * beta * t + 1.0)
g_0 = lambda t: u_true(0.0, t)
g_1 = lambda t: u_true(L, t)

# Discretization and output times
m = 20
x = numpy.linspace(0.0, L, m + 2)
x_fine = numpy.linspace(0.0, L, 100)
delta_x = L / float(m + 1)
output_times = (0.0, 0.5, 1.0, 1.5, 2.0)

# Solve
U_CN = numpy.empty((len(output_times), m + 2))
U_CN[0, :] = u_true(x, 0.0)
for (n, t) in enumerate(output_times[1:]):
    U_CN[n + 1, :] = solve_heat_CN(m, kappa, L, output_times[n], t, U_CN[n, :], g_0, g_1)

error = numpy.linalg.norm(delta_x * (U_CN[-1, :] - u_true(x, output_times[-1])), ord=1)
print "Error CN = %s" % error
assert error < 1e-3
print "Success!"

# Plot some of the results
colors = ['k', 'r', 'b', 'g', 'c']
fig = plt.figure()
axes_CN = fig.add_subplot(1, 1, 1)
for (n, t) in enumerate(output_times):
    axes_CN.plot(x_fine, u_true(x_fine, t), 'k-')
    axes_CN.plot(x, U_CN[n, :], "o%s" % colors[n], label='t=%s' % numpy.round(t, 4))
    axes_CN.set_xlabel("x")
    axes_CN.set_ylabel("$U_(x,t)$")
    axes_CN.set_title("Solution to Heat Equation")
    axes_CN.set_xlim([0, L])
    axes_CN.legend()

plt.tight_layout()
plt.show()

Error CN = 0.00774874523861

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-5-ca2094b702f3> in <module>()
     21 error = numpy.linalg.norm(delta_x * (U_CN[-1, :] - u_true(x, output_times[-1])), ord=1)
     22 print "Error CN = %s" % error
---> 23 assert error < 1e-3
     24 print "Success!"
     25 

AssertionError:

(b) [5] Perform a convergence analysis on the methods above (i.e. plot their convergence).

# YOUR CODE HERE
raise NotImplementedError()

Question 2 - Periodic boundary conditions¶

Consider the temperature distribution in a thin circular ring with perimeter $2L$. Assume that the lateral surfaces are perfectly insulated, such that the temparature is distributed evenly througout the ring (see sketch) and that the ends at $x=-L$ and $x=L$ are in perfect thermal contact. Thin ring, may not show up if you are not online.

(a) [5] Use the heat equation to formulate the corresponding mathematical problem.

Since we consider that the ring is thin, we can interpret the problem as heat distribution in a circular wire. We begin with the general form of the heat equation: $$ u_t = \kappa u_{xx} $$ Where $\kappa$ is the conductivity of the ring material. We therefore have a problem of heat diffusion, where we can apply the following boundary conditions: $$ u(-L,t) = u(L,t) \space and \space \frac {du} {dx} (-L,t) = \frac {du} {dx} (L,t) $$ And the following initial conditiosn: $$ u(x,0) = f(x)$$ for a given f(x)

(b) [10] Solve the problem analytically.

We begin with a seperation of variables to solve the 2nd order PDE: $$ u(x,t) = \phi(x)G(x) \\ \frac{\phi''(x)}{\phi(x)} = \frac {G\dot(t)}{\kappa G(t)} = - \lambda^2 \\ $$ IVP: $$ \frac {G\dot(t)}{\kappa G(t)} = - \lambda^2 \\ G(x) = ce^{-\kappa\lambda^2t} $$ BVP: $$ \frac{\phi''(x)}{\phi(x)} = - \lambda^2 \\ \phi''(x) + \lambda^2 \phi(x) = 0 \\ \phi(x) = Acos(\lambda x) + Bsin(\lambda x) \\ $$ We determine $\lambda$ for a non-trivial solution as being $\lambda = \frac{n\pi}{L}$ We therefore obtain the solution: $$ u(x,t) = e^{-\kappa(\frac{n\pi}{L})^2t}[\sum_{n=0}^{\infty} A_ncos(\frac{n\pi}{L}x) + \sum_{n=1}^{\infty} B_nsin(\frac{n\pi}{L}x)] \\$$ We then use the orthogonality of cosine and sine to determine the coefficients $A_0$, $A_n$ and $B_n$: $$ A_0 = \frac{1}{2L} \int_{-L}^L f(x)dx \\ A_n = \frac{1}{2L} \int_{-L}^L f(x)cos(\frac{n\pi x}{L})dx \space for \space n = 1,2,3... \\ B_n = \frac{1}{2L} \int_{-L}^L f(x)sin(\frac{n\pi x}{L})dx \space for \space n = 1,2,3... \\ $$

(c) [15] Write a function to solve this problem using the Crank-Nicholson method.

def solve_heat_periodic_CN(m, kappa, L, t_0, t_final, U_0):
    """Solve the heat equation on a periodic domain using Crank-Nicholson
    
    :Input:
     - *m* (int) Number of points use to discretize the domain.  Note that
       the total number of points is *m+1*.
     - *kappa* (float) Diffusion coefficient
     - *L* (float) Length of half of the domain
     - *t_0* (float) Starting time
     - *t_final* (float) Time to integrate to
     - *U_0* (numpy.ndarray) Initial condition at time t_0, should be m+1
     
    :Output:
     - (numpy.ndarray) Solution at time t_final.  Note that this vector should m+1
    """
    
    # YOUR CODE HERE
    delta_x = L / float(m+1)
    
    C = 0.2
    
    delta_t = C * delta_x
    print U_0[0]
    
    r = numpy.ones(m+2) * delta_t * kappa / (2.0 * delta_x**2)
    A = sparse.spdiags([-r, 1.0 + 2.0 * r, -r], [-1, 0, 1], m+1, m+1).tocsr()
    C = numpy.zeros((m+1,m+1))
    C[0][0] += 1.0
    C[-1][-1] += 1.0
    
    A = A + C
    
    B = sparse.spdiags([r, 1.0 - 2.0 * r, r], [-1, 0, 1],  m+1, m+1).tocsr()
    
    b = B.dot(U_0)
    b[0] += delta_t * kappa/ (2.0 * delta_x**2) * (U_0(0) + g_0(t_final))
    b[-1] += delta_t * kappa / (2.0 * delta_x**2) * (g_1(t_0) + g_1(t_final))
    
        # Solve system
        
    
    U = linalg.spsolve(A, b)
    #raise NotImplementedError()
    
    return U

L = numpy.pi
kappa = 0.075
u_true = lambda x, t: 6.0 * numpy.sin(2.0 * numpy.pi / L * x) * numpy.exp(-4.0 * kappa * (numpy.pi / L)**2 * t) + 3.0

# Discretization and output times
# Note that there are m unknowns including one of the points either at 0 or L
m = 200
delta_x = (2.0 * L) / float(m)
x = numpy.linspace(-L, L, m + 2)
x_fine = numpy.linspace(-L, L, 100)
output_times = (0.0, 2.5, 5.0, 7.5, 10.0)

# Solve
U = numpy.empty((len(output_times), m + 1))
U[0, :] = u_true(x[:-1], 0.0)
for (n, t) in enumerate(output_times[1:]):
    U[n + 1, :] = solve_heat_periodic_CN(m, kappa, L, output_times[n], t, U[n, :])

error = numpy.linalg.norm(delta_x * (U[-1, :] - u_true(x[:-1], output_times[-1])), ord=1)
print "Error = %s" % error
assert error < 1e-3
print "Success!"

# Plot some of the results
colors = ['k', 'r', 'b', 'g', 'c']
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
for (n, t) in enumerate(output_times):
    axes.plot(x_fine, u_true(x_fine, t), 'k-')
    axes.plot(x[:-1], U[n, :], "o%s" % colors[n], label='t=%s' % numpy.round(t, 4))
    axes.set_xlabel("x")
    axes.set_ylabel("$U_(x,t)$")
    axes.set_title("Solution to Heat Equation")
    axes.set_xlim([-L, L])
    axes.legend()
plt.show()

Question 3 - Diffusion and Decay¶

Consider the partial differential equation

$$ \frac{\partial u}{\partial t} = \kappa \frac{\partial^2 u}{\partial x^2} -\gamma u, $$

which models a diffusion with parameter $\kappa \in \mathbb R_{>0}$ with decay specified by $\gamma \in \mathbb R_{>0}$.

Consider a $\theta$-method of the form

$$ U_j^{n+1} = U_j^n + \frac {\Delta t} {2\Delta x^2} \left( U_{j-1}^n - 2 U_j^n + U_{j+1}^n + U_{j-1}^{n+1} - 2U_j^{n+1}+U_{j+1}^{n+1} \right)-\Delta t\gamma\left( (1-\theta) U_j +\theta U_j^{n+1} \right) $$

where $\theta$ is a parameter.

(a) [15] Assume $\kappa = 1 $, by computing the local truncation error, show that this method is $p$th order accurate in time and second order accurate in space, where $p=2$ for $\theta = 0.5$ and $p=1$ otherwise.

YOUR ANSWER HERE

(b) [15] Using von Neumann analysis, show that this method is unconditionally stable if $\theta \geq 0.5$.

We begin by substituting the function $U^n_j = e^{ij\Delta x \xi}$ into the above method, which leads to: $$ U^{n+1}_j = U^n_j + \frac{\Delta t}{2 \Delta x^2} \left[ \left(e^{i \Delta x \xi} - 2 + e^{-i\Delta x \xi} \right) U^n_{j}+ \left(e^{i \Delta x \xi} - 2 + e^{-i\Delta x \xi} \right) U^{n+1}_{j} \right ]-\Delta t\gamma\left( (1-\theta) U_j +\theta U_j^{n+1} \right) $$ We note the substitution $U^{n+1}_j = g(\xi) U^n_j$ which leads to: $$ g = 1 + \frac{\Delta t}{2 \Delta x^2} \left( e^{i \Delta x \xi} - 2 + e^{-i\Delta x \xi} \right ) (1 + g)-\Delta t\gamma\left( (1-\theta)+\theta g \right) $$

By solving for the amplification factor g, we obtain that:

$$ g = \frac{1+ z -\Delta t \gamma (1-\theta)}{1 - z + \Delta t \gamma \theta} $$

where $$ z = \frac{\Delta t}{\Delta x^2} \left (e^{i \Delta x \xi} - 2 + e^{-i\Delta x \xi} \right ) = \frac{2 \Delta t}{\Delta x^2} (\cos(\xi \Delta x) - 1). $$ We therefore would like to ensure that $ -1 \leq g \leq 1$ for stability: $$ \frac{1+ z -\Delta t \gamma (1-\theta)}{1 - z + \Delta t \gamma \theta} \leq 1 \\ 2z \leq \Delta t \gamma \\ $$ Which holds so long as $ \Delta t \geq \frac {2z}{\gamma} $. Similarly $$ \frac{1+ z -\Delta t \gamma (1-\theta)}{1 - z + \Delta t \gamma \theta} \geq -1 \\ 2z \leq \Delta t \gamma \\ $$

(c) [15] Show that, if $\theta = 0$, then the method is stable provided $\Delta t \leq \frac 2 \gamma$, independent of $\Delta x$.

If $ \theta = 0$, we have the following expression:

$$ g = \frac{1+ z -\Delta t \gamma}{1 - z} $$

We must ensure that $ -1 \leq g \leq 1$ for stability, therefore we have: $$ \frac{1+ z -\Delta t \gamma}{1 - z} \leq 1 \\ \Delta t \gamma \geq 2z \rightarrow \Delta t \geq 0 \\ $$

or

$$ \frac{1+ z -\Delta t \gamma}{1 - z} \geq -1 \\ \Delta t \leq \frac 2 \gamma $$