Michael Holst  
http://ccom.ucsd.edu/~mholst/ 
Professor of Mathematics and Physics UC San Diego 


AMa 220C (Numerical Continuation Methods)
Course Topics: Nonlinear PDEs, bifurcation theory,
numerical methods for continuation
DUE DATE: The PROJECT is do at the end of the term. NOTE: (5497) A couple of crises have delayed my getting a formal writeup of this posted; I'll have a LaTeX file up here in the next few days. Meanwhile, get familiar with the MATLAB FEM package above. I reworked it this weekend to put the abstract nonlinear/bilinear form interface on it, so it will be really easy to work with. The Bratu problem on the unit square is currently setup as the problem (meaning that I provide F(u)(v), DF(u)(w,v), and the initial simplicial mesh). The parameter lambda is currently chosen so that the nonlinear system is uniquely solvable. The package uses the form F(u)(v) to build nonlinear algebraic equations, and then does a Newton iteration using the bilinear form DF(u)(w,v), as described in class. Our task in the PROJECT will simply be to form the augmented (bordered) system to handle the fold that occurs on the solution path. We will first look at continuation in lambda (natural continuation), notice that it fails, and then implement pseudoarclength continuation. This will pretty much boil down to sticking an extra row and column (the right choice of row and column is of course the whole game) onto the Jacobian; an easy task in MATLAB. NOTE: (52997) Tobias found a mistake in my definition of DF(u)(w,v) (the routine DFu_wv.m) for the Bratu problem in the posted version of FEMBIF. I corrected the problem and posted the new version at the link above. In implementing pseudoarclength continuation in FEMBIF, just use the simplest predictorcorrector iteration with an Euler predicter, and a Newton corrector. As discussed in class Wednesday (52897), you should run into the turning point around lambda=7, with some variation depending on how accurately you are discretizing the problem (how many mesh refinements you are doing before solving the discrete problem). NOTE: (53097) In your interpretation of the innerproducts <.,.> for implementing the augmented system, keep in mind that while the correct innerproduct for two finite element functions is multiplication and then integration over the domain, the correct innerproduct of two real numbers (i.e., the constraint) can be simply multiplication. Also, keep in mind that I keep the dirichlet points in the list of "unknowns" in FEMBIF; the equations for these points (which are really known values) that appear in the discrete system are then just identity equations. This is commonly done in finite element implementations to avoid dealing with two different types of points. When you add the constraints (the new column and row), keep in mind that any entry in the new column that is at a dirichlet point is simply zero, and any entry in the row that will multiply a dirichlet point in an innerproduct will also need to simply be set to zero. If you are confused about this, come to my office and I'll make it clear on the board. NOTE: (6997) Believe it or not, it is actually possible to get the FEMBIF program working; my solution is now posted as the MATLAB program at the link above. I do a simple Eulerpredictor, Newtoncorrector, with pseudoarclength continuation, as discussed in class. I *DO* use the innerproduct we described the last day we met; namely the big L2 innerproduct for finite element functions (rather than the H1 semiinnerproduct as in Glowinski et al.), and simply multiplication as an innerproduct for scalars. I now believe that the problems people were encountering had to do with how the Dirichlet conditions were handled in the Jacobian, rather than a problem with the innerproduct. Look at my code above to see how I deal with it. Note that the sign convention used in the writeup is slightly different than that used in the FEMBIF code; make sure you take this into account, or you will get weird results (like one unique solution, no turning point). Using this implementation, it is possible to recover the Sshaped curve (lambda versus u), and the turning point approaches the value in Herb's papers as the mesh is refined. (I believe now that the slightly different value some of you were seeing was due to the coarsness of the mesh you were using.) I realize this was a very difficult problem to do, but everyone at least made some progress toward a solution (a couple of you actually got the thing working!) In any event, have a look at my solution above; I tried to document the code reasonably well, so it should be possible to understand what I am doing. Thanks for attending my lectures this quarter; things were a bit disorganized sometimes due to several trips I had to take, my car accident, etc, but I learned some things and I hope you guys did too. Have a great summer!
