DUE DATE: The PROJECT is do at the end of the term.

NOTE: (5-4-97) 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 pseudo-arclength 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: (5-29-97) 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 pseudo-arclength continuation in FEMBIF, just use the simplest predictor-corrector iteration with an Euler predicter, and a Newton corrector. As discussed in class Wednesday (5-28-97), 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: (5-30-97) In your interpretation of the inner-products <.,.> for implementing the augmented system, keep in mind that while the correct inner-product for two finite element functions is multiplication and then integration over the domain, the correct inner-product 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 inner-product 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: (6-9-97) 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 Euler-predictor, Newton-corrector, with pseudo-arc-length continuation, as discussed in class. I *DO* use the inner-product we described the last day we met; namely the big L2 inner-product for finite element functions (rather than the H1 semi-inner-product as in Glowinski et al.), and simply multiplication as an inner-product 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 inner-product. Look at my code above to see how I deal with it.

Note that the sign convention used in the write-up 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 S-shaped 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!