#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 25 12:42:01 2018

@author: chcurry
"""

# Plots numerical against exact solutions for two different methods
# visc*u_xx - u_x = -1
# u(0)=alpha, u(1)=beta

# Note that the central difference approximation has a nasty oscillation
# near the right boundary

# experiment with different values of visc and m. What do you see?

import numpy as np
import scipy as sp
from matplotlib import pyplot as plt

def viscosity_plot_central(visc):
    xx = np.linspace(h,T-h,m)
    u = alpha + xx + (beta-alpha-1)*(np.exp(xx/visc)-1)/(np.exp(1/visc)-1)
    plt.figure()
    plt.plot(np.hstack([0,xx,1]),np.hstack([alpha,u,beta]))
    A = (1/h**2)*sp.sparse.diags([1,-2,1],[-1,0,1],shape=(m,m))
    D = (1/(2*h))*sp.sparse.diags([-1,1],[-1,1],shape=(m,m))
    F = -np.ones(m)
    F[0] = F[0] - alpha*visc/h**2 - alpha/(2*h)
    F[-1] = F[-1] - beta*visc/h**2
    Am = visc*A - D
    U = sp.sparse.linalg.spsolve(Am,F)
    plt.plot(np.hstack([0,xx,1]),np.hstack([alpha,U,beta]))
    plt.show()
    
def viscosity_plot_backward(visc):
    xx = np.linspace(h,T-h,m)
    u = alpha + xx + (beta-alpha-1)*(np.exp(xx/visc)-1)/(np.exp(1/visc)-1)
    plt.figure()
    plt.plot(np.hstack([0,xx,1]),np.hstack([alpha,u,beta]))
    A = (1/h**2)*sp.sparse.diags([1,-2,1],[-1,0,1],shape=(m,m))
    D2 = (1/h)*sp.sparse.diags([-1,1],[-1,0],shape=(m,m))
    F2 = -np.ones(m)
    F2[0] = F2[0] - alpha*visc/h**2 - alpha/h
    F2[-1] = F2[-1] - beta*visc/h**2
    Am2 = visc*A - D2
    V = sp.sparse.linalg.spsolve(Am2,F2)
    plt.plot(np.hstack([0,xx,1]),np.hstack([alpha,V,beta]))
    plt.show()    
    
T=1
alpha=1
beta=3
m=20
h=T/(m+1)    
viscosity_plot_backward(0.01)
viscosity_plot_central(0.01)