|
AMa 220C (Numerical Continuation Methods)
Course Topics: Nonlinear PDEs, bifurcation theory,
numerical methods for continuation
Instructors: Dr. M. Holst and Prof. H. Keller
Term: Spring 1997
Outline:
Bifurcation Theory and Numerical Methods
Software:
MATLAB Finite Element Bifurcation Package
for 2D Nonlinear Elliptic Equations
A writeup of our class PROJECT can be found
here
(5-4-97).
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!
|