<jupyter_start><jupyter_text>Conjugate-Gradient Method[AMath 585, Winter Quarter 2020](http://staff.washington.edu/rjl/classes/am585w2020/) at the University of Washington. Developed by <NAME> and distributed under the [BSD license](https://github.com/rjleveque/amath585w2020/blob/master/LICENSE).  You are free to modify and use as you please, with attribution.These notebooks are all [available on Github](https://github.com/rjleveque/amath585w2020/).-----This notebook illustrates an implmentation of the Conjugate-Gradient method on a simple 1-dimensional boundary value problem $u''(x) = f(x)$ with Dirichlet boundary conditions.  This illustrates the basic algorithm, although for this particular problem it would of course be better to solve the tridiagonal system directly.<jupyter_code>%matplotlib inline
from pylab import *<jupyter_output><empty_output><jupyter_text>Animating the iterationHere we allow two possible ways to animate figures.  The version defined in the module [jsanimate_figs.py](jsanimate_figs.py) creates an animation that works well even when the notebook is rendered as html (by embedding javascript into the notebook, hence the js). The code to do this is a bit long so it has been put in a separate file, but it is based on the matplotlib [animation tools](https://matplotlib.org/3.1.3/api/animation_api.html), in particular [FuncAnimation](https://matplotlib.org/api/_as_gen/matplotlib.animation.FuncAnimation.html) to loop over figures and the `to_jshtml` function to convert to javascript.The version using [ipywidgets](https://ipywidgets.readthedocs.io/en/latest/examples/Widget%20Basics.html) does not render in the html version but you may prefer using it if you are running the notebook live, since it executes more quickly.To use widgets, set `use_widgets = True`.<jupyter_code>use_widgets = True

if use_widgets:
    from ipywidgets import interact
    import ipywidgets as widgets

    def animate_figs(figs):
        show_frame = lambda frameno: display(figs[frameno])
        interact(show_frame, frameno=widgets.IntSlider(min=0,max=len(figs)-1, value=0))
else:
    from jsanimate_figs import animate_figs<jupyter_output><empty_output><jupyter_text>Define the BVPAgain we solve $u''(x) = f(x)$ with Dirichlet boundary conditions and choose a problem where the truncation error is zero so that the exact solution of the linear system is the solution of the ODE evaluated at the grid points.  Discrete matrixIn the notebook [IterativeMethods.ipynb](IterativeMethods.ipynb) we solved the system\begin{align*}\frac{1}{h^2} \left[\begin{array}{ccccccccccccccc}1&0\\1&-2&1\\&1&-2&1\\&&&\ddots\\&&&&1&-2&1\\&&&&&0&1\end{array}\right]~\left[  \begin{array}{ccccccccccccccc}U_0 \\ U_1 \\ U_2 \\ \vdots \\ U_m \\ U_{m+1}\end{array} \right]= \left[  \begin{array}{ccccccccccccccc}\alpha \\ f(x_1) \\ f(x_2) \\ \vdots \\ f(x_m) \\ \beta\end{array} \right].\end{align*}But note that this matrix is not symmetric and hence the conjugate gradient method cannot be used, or at least it will not work for general initial data.  It turns out if the initial data satisfies the Dirichlet boundary conditions then it does in fact work, because the first and last equations are always satisfied, but to avoid concerns that this is affecting the behavior we will instead solve the tridiagonal system for the interior unknowns,\begin{align*}\frac{1}{h^2} \left[\begin{array}{ccccccccccccccc}-2&1\\1&-2&1\\&&&\ddots\\&&&&1&-2&1\\&&&&&1&-2\end{array}\right]~\left[  \begin{array}{ccccccccccccccc}U_1 \\ U_2 \\ \vdots \\ U_{m-1} \\ U_m\end{array} \right]= \left[  \begin{array}{ccccccccccccccc}f(x_1) - \alpha/h^2 \\ f(x_2) \\ \vdots \\ f(x_{m-1}) \\ f(x_m) - \beta/h^2\end{array} \right].\end{align*} This complicates things a bit since we generally want to plot the full solution so we will use `U_full` to be the full set of grid points, including boundary points, and `U` to be the solution on the interior points, the actual solution to the linear system, with similar notation for other variables. **Note:** This matrix $A$ is symmetric *negative* definite rather than positive definite.  We could negate both the matrix and right hand side to make it positive definite, but this isn't necessary:  the C-G algorithm as written in the text still works fine.  Also note that the first matrix above, in addition to being non-symmetric, has two positive eigenvalues equal to 1 and the rest of the eigenvalues are negative.  So even if it were symmetric, C-G would not work without negating some of the equations.  Define the matrix-vector multiply functionTo implement the C-G algorithm we need to be able to multiply $Av$ for various vectors $v$.  The following `matvec` routine computes and returns this product.  This is the only place the matrix $A$ is needed, so not that in general we do not need to form or store the large sparse matrix.<jupyter_code>def matvec(v):
    """
    Given v of length m 
    Return b = A*v
    """
    
    m = len(v)
    h = 1./(m+1)
    
    # initialize vector for product:
    b = empty(m, dtype=float)
    
    # set each value of b, noting that first and last rows are special:
    for i in range(m):
        if i>0:
            vim = v[i-1]
        else:
            vim = 0.
        if i<m-1:
            vip = v[i+1]
        else:
            vip = 0.
            
        b[i] = (vim - 2*v[i] + vip) / h**2
        
    return b<jupyter_output><empty_output><jupyter_text>Different version: extend v to include boundary pointsHere's a simpler way to define the `matvec` function that is more like what is done in the notebook [ConjugateGradient2D.ipynb](ConjugateGradient2D.ipynb).<jupyter_code>def matvec(v):
    """
    Given v of length m 
    Return b = A*v
    """
    
    # pad v with zeros at each boundary, needed for computing centered differences:

    m = len(v)
    h = 1./(m+1)
    v_full = zeros(m+2)
    v_full[1:-1] = v
    
    # compute b of length m at interior points only. 
    # Note that slicing with [1:-1] corresponds to interior points,
    # [:-2] correspond to points to the left and 
    # [2:] are points to the right.
    
    b = (v_full[:-2] - 2*v_full[1:-1] + v_full[2:]) / h**2
    
    return b<jupyter_output><empty_output><jupyter_text>The Conjugate-Gradient algorithmThe next function implements the C-G algorithm.<jupyter_code>def solve_bvp_CG(f_fcn, utrue_fcn, m, maxiter, kplot, verbose=False):
    
    h = 1./(m+1)
    x_full = linspace(0,1,m+2)
    x = x_full[1:m+1] # interior points

    utrue_full = utrue_fcn(x_full)
    utrue = utrue_full[1:m+1]  # at interior points
    
    # Dirichlet boundary values from true solution:
    alpha = utrue_fcn(0.)
    beta = utrue_fcn(1.)

    f_full = f_fcn(x_full)

    # right-hand side:
    f = f_full[1:m+1]  # at interior points

    # adjust for Dirichlet BCs:
    f[0] = f[0] - alpha / h**2
    f[m-1] = f[m-1] - beta / h**2

    # initial guess:
    U0_full = linspace(alpha, beta, m+2)  # linear
    U0 = U0_full[1:m+1]  # interior points

    U = U0.copy() # current iterate
    r = f - matvec(U)  # initial residual
    p = r.copy()  # initial direction

    tol = 1e-8  # stop if the residual 

    enorm = abs(U-utrue).max()
    errors = [enorm]
    figs = []  # for the list of figures we generate

    rTr_km = dot(r,r)  # r^T * r at iteration k-1

    for k in range(1,maxiter+1):
        w = matvec(p)   # the only matrix-vector multiply
        a = rTr_km / dot(p,w) # alpha_{k-1} in CG algorithm
        U = U + a*p
        r = r - a*w

        enorm = abs(U-utrue).max()
        errors.append(enorm)

        if mod(k,kplot)==0 or k==maxiter:
            # every kplot iterations create a plot:
            fig = figure(figsize=(12,5))
            plot(x_full,U0_full,'r-o', label='initial guess')
            plot(x_full,utrue_full,'k-o', label='true solution')
            U_full = hstack([alpha,U,beta])
            plot(x_full,U_full,'bo-', label= 'iteration k = %i' % k)
            legend()
            grid(True)
            xlim(0,1)
            ylim(0,3)
            title('After %i iterations, norm(error) = %.2e' \
                  % (k, enorm))
            figs.append(fig)
            close(fig)

        rTr_k = dot(r,r)
        rnorm = sqrt(rTr_k)
        
        if verbose:
            print('iteration %3i:   2-norm(r) = %.2e,   max-norm(E) = %.2e' \
                  % (k,rnorm,enorm))
        
        if rnorm < tol:
            if verbose: print('Stopping after %i iterations' % k)
            break

        # determine next search direction:
        b = rTr_k / rTr_km   # beta_{k-1} in CG algorithm
        rTr_km = rTr_k       # for next iteration
        p = r + b*p           # next search direction
         
    return errors, figs<jupyter_output><empty_output><jupyter_text>Test this on a cubic functionAgain we choose an example for which the exact solution of the linear system is just the ODE solution evaluated at the grid points...<jupyter_code>f_fcn = lambda x: 6*x + 2
utrue_fcn = lambda x: x**3 + x**2 - x + 1
errors,figs = solve_bvp_CG(f_fcn, utrue_fcn, m=19, maxiter=25, kplot=1, verbose=True)<jupyter_output><empty_output><jupyter_text>Note that this converged after 19 iterations, as expected.  Note that neither the residual nor the error were very small in iteration 18, and dropped dramatically in the final iteration!<jupyter_code>semilogy(range(1,len(errors)+1), errors, 'b-x')
grid(True)<jupyter_output><empty_output><jupyter_text>Note the behavior of the approximate solution, now plotted every iteration:<jupyter_code>animate_figs(figs)<jupyter_output><empty_output><jupyter_text>A different ODEHere's an even simpler ODE where the true solution is quadratic.  In this case we obtain convergence in 10 iterations, again with a big drop in the final iteration.<jupyter_code>f_fcn = lambda x: -6*ones(x.shape)
utrue_fcn = lambda x: 2-3*(x-0.2)*(x-0.7)
errors,figs = solve_bvp_CG(f_fcn, utrue_fcn, m=19, maxiter=25, kplot=1, verbose=True)
semilogy(range(1,len(errors)+1), errors, 'b-x')
grid(True)<jupyter_output><empty_output><jupyter_text>The behavior of the iterates is particularly simple in this case:<jupyter_code>animate_figs(figs)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1105660747571863, 'fraction_numerical': 0.018149710585696067, 'mean_word_length': 3.3048986486486487, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 43, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '635034', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3489, 'n_tokens_neox': 3283, 'n_words': 1099, 'path': 'notebooks/ConjugateGradient.ipynb', 'repo_name': 'rjleveque/amath585w2020', 'size': '10199'}