The following content is provided under a Creative
Commons license.
Your support will help MIT OpenCourseWare
continue to offer high quality educational resources for free.
To make a donation or to view additional materials
from hundreds of MIT courses, visit MIT OpenCourseWare
at ocw.mit.edu.
PROFESSOR: OK, let's go ahead and get started again.
So we had a nice review on Monday
of some of the different concepts
that we've talked about up to this point in the course.
We've covered linear algebra, solutions of systems
of nonlinear equations.
There's been a lot of fundamentals involved in that.
Things are going to get more applied
as the sorts of problems we're trying
to solve get more complicated.
Complicated problems require very sophisticated methods
to solve them.
It's not going to be our expectation that you're
necessarily able to implement all
of these sophisticated methods.
So for example, you implemented Newton's method
on your last homework assignment.
That'll be the last time we ask you to do that by hand.
You'll utilize tools built into Matlab
to solve other systems of nonlinear equations going
forward.
You understand how it works and what the limitations are.
You understand how to build different tools that
will be inputs to these functions
like the function you're trying to find the root for,
or an explicit expression for the Jacobian.
You even understand how to set certain parameters associated
with these methods, and what those particular parameter
values mean.
Like setting the function norm tolerance,
or setting the step norm tolerance
associated with these methods.
Setting the number of iterations you're
going to let the method go.
You've implemented it, so you understand
what these things mean, what the consequences of choosing them
are, how they might increase the time associated
with computations that you're going to conduct.
So we've done a lot of fundamentals.
Maybe that was difficult for some of you
who haven't seen this before.
Now things are becoming more about the practice
of numerical methods.
How do we actually solve things reliably,
and that's what today's lecture will talk about.
We discussed solving systems of nonlinear equations
and needing good initial guesses for our methods.
Bad initial guesses can lead to unpredictable results.
Good initial guesses, we know the Newton-Raphson method
is going to double the number of accurate digits each iteration,
and we should converge very rapidly to a solution.
So it seems like having good initial guesses is
really ideal.
And this topic, homotopy and bifurcation,
will cover precisely those issues,
and the practice of getting these good initial guesses.
I think is one of the most useful things you'll
learn in this class, actually.
It's a really nice way of going about solving
complicated problems.
Before we begin, I want to address a question
I got via email last night.
That's always good if these questions
go on [? Piaza ?] so everyone can read them,
but I thought it was a very nice question about how
do you find-- the question was how
do you find multiple roots using the Newton-Raphson method.
So if I have a polynomial, where the function
doesn't cross through the origin,
but it just touches the origin, and goes back up.
That's a root with multiplicity in the polynomial.
And can the Newton-Raphson method find that?
Because the derivative of the function is 0 at the root.
So who says the Newton-Raphson method will
be able to find those roots?
What do you think?
Yes?
And who says it won't be able to find those roots?
Or who doesn't know?
You don't know.
That's OK.
That's a good question, right.
Do you know or don't you know?
So we know the Newton-Raphson method
relies on knowing the derivative of the function
at different points in order to take
our next Newton-Raphson step towards that solution.
It may be the case, the derivative is 0 at the root,
but nearby the root, the derivative may not be 0.
And then the Newton-Raphson steps are going to be fine.
They're are always going to keep stepping down towards the root.
If you go through an analysis of how the algorithm converges,
you'll find out that you only get linear convergence
in that case, instead of quadratic convergence.
So there will be some penalty to pay for the fact
that the derivative is 0 at the root
itself, but it's not going to be so bad.
The method can still converge to those sorts of solutions.
Not all is lost there, but the rate of convergence
may be not as good as you expected it to be.
And the same is true in the multi-dimensional case as well.
It's a little hard to imagine what a multiple root is
going to look like.
It's something like tangent curves
in the solution of plane.
So if I have f1 of x1 and x2 equal 0,
that's a curve in the x1, x2 plane.
And I have f2 in that same plane equal to 0.
The place where they touch might be one of these multiple roots.
That's a place where the determinant of the Jacobian
may be equal to 0.
You'll have the same sort of circumstance.
Nearby, the Jacobian can be non singular,
and you can take steps towards the root,
but there's a worst case scenario where you're
coming in on a direction that's parallel
that belongs to the null space of the Jacobian at the root,
and then everything's going to slow down.
You can still converge, but it may
be very slow as a consequence.
I thought that was a really nice question.
It's a detailed question about how Newton-Raphson works.
It's not necessarily going to fail
under these weird circumstances.
They don't come up too often.
They're going to come up in our discussion today, actually.
It's not going to fail, but it may converge very slowly
to the solution.
OK?
Good.
I'm going to remind you, the last thing we
did with Newton-Raphson methods, we
talked about one way of fixing them.
In lieu of having good initial guesses,
we may not want to take the full Newton-Raphson step.
We may want to do these quasi Newton-Raphson methods instead.
One of those was this damp Newton-Raphson
where we introduced some damping parameter.
We don't take the full step, we take some shorter step,
and we try to make the absolute value over function
become as small as possible.
That's an optimization problem.
It's as hard as finding the root of the system of equations
that we're looking at.
So we don't solve that optimization problem exactly.
That's crazy.
Instead we use this backtracking line search method
where we change our damping parameter
in some systematic way until we get to a point
where we're satisfied.
Where we say OK, in the process of changing this,
we've actually found a smaller value of our function,
and we accept that as our best guess
for where the next iterate is supposed to go.
And this sort of approach can give you
global convergence for the algorithm, which is great,
but it's not globally convergent to roots.
It's globally convergent to minima, asymptotes, and roots
as well.
But it's globally convergent.
So it'll terminate.
The algorithm will stop, and when it stops,
you can ask is this a root or not?
Do I think this is sufficiently close to a root or not?
If it's not sufficiently close to a root,
then you know you terminated somewhere like a minima,
or you rushed out towards an asymptote,
and maybe you need to start with a new initial guess
for your solution.
These quasi Newton-Raphson methods are very useful.
This is what Matlab uses in order to--
it uses a more sophisticated version
of this called the Dogleg method,
but it's the same principle.
It's trying to make sure you're not
taking steps that are too big.
Trying to make sure you're taking
steps in the right direction so that you're approaching
the roots of your function.
Or at least reducing the absolute value of your function
as close to zero as you can get it.
Are there any questions about this?
Anything we need to follow up on in the previous lecture?
No.
OK, good initial guesses.
This is the key, really, for solving
systems of nonlinear equations.
It's the key for doing optimization problems too.
We need to have some idea of what our solution looks like,
and we want guesses that are close to those solutions
so that maybe we are in this region of local convergence
for the Newton-Raphson method, and everything
converges very quickly.
So where do they come from?
This is an important question, and I'll
show you how you can get them.
There's another issue too, which is related
to this, which is the non-linear equations,
they can have multiple roots.
If we do optimization problems-- that's our next topic--
they can have multiple minima or maxima
that we're trying to seek out.
Depending on our initial guess, our numerical method
may find one of these roots or another root.
It may find one of these minima or another minima.
Maybe we're looking for the global minimum,
so we need to look at them all.
Are there methods that we can use to find them all?
So we'll talk about continuation in homotopy,
and we'll talk about bifurcation.
And it turns out there are methods
one can use to find all the roots for a system
of non-linear equations.
In fact, the homework assignment you have this week
will ask you to do ternary phase equilibrium problem.
Vapor liquid phase equilibrium with three components.
There are five parts to this problem.
Four of them, I think, are relatively straightforward.
They're about finding the pure component roots
or randomly guessing some initial conditions
for your nonlinear solver in order
to find all of the different azeotropes in the system.
There's a fifth part that asks you to actually implement
this homotopy method, and try to find all the roots in one
go with your algorithm.
Like the other homeworks you've had,
that part's a little more tricky.
Usually these homeworks have one thing
that's a little bit more tricky than everything else.
That part is a little bit more tricky,
but it's really instructive to try to do it.
We'll learn how to do it today.
You can find all of the coexisting compositions
of this ternary phase diagram, simultaneously,
with one go through this algorithm.
So let's talk about continuation first.
This is an interesting idea, and it
relates to trying to find roots of these equations,
having good initial guesses.
Here's a simple example.
We have a third degree polynomial,
and it has three roots.
We can look, graphically, because this is one dimension,
and see that one of these roots is near minus 1 and 1/2,
another root is near a 1/2, and another root is near 1.
And so I can go to my non-linear equation solver,
and give it initial guesses, minus 1/2, 1/2, and 1,
and I'm probably going to find these roots pretty reliably.
But in many dimensions, we can't do these sorts
of graphical approaches.
You can't look at this function.
It may be difficult to actually get
enough data to even look at it if it's a two-dimensional set
of functions.
It's hard to see where these roots are,
and so continuation is a method of transforming the problem
into an easier one to solve.
So let's take this problem, this third order polynomial,
and let's transform it.
Let's replace the 2 with a 2 times the lambda,
and let's try to find the roots of this polynomials
as lambda grows from 0 to 1.
This is sort of a weird idea, but we actually
know when lambda is equal to 0, this goes away,
and there's only one root to this polynomial-- one
real root-- which is minus 1.
We actually know the answer when lambda is equal to 0.
If I make lambda a little bit bigger than 0,
there will be a root which is close to minus 1.
In fact, minus 1 is a good initial guess for that root,
and my non-linear equation solver
ought to converge very quickly to that value.
If I choose lambda to be small, then it
should be in that region of local convergence
of my Newton-Raphson method, and I very quickly
should find the next root.
So I can step lambda along.
I use the lambda equals 0 root as an initial guess
for the next root, and I keep increasing lambda
from 0 all the way up to 1, and when lambda is equal to 1,
I'm going to find the root to the initial problem I
was looking at.
Does that makes sense?
So I'm continuing from one set of solutions to the next.
I vary this parameter continuously,
and for previous values of that parameter,
they serve as initial guesses for the next solution.
Here's what a code that does that would look like.
I define my lambdas.
We don't know good ways to define the lambda.
I've got to figure that out for myself,
but here lambda is a vector that goes from 0 to 1,
in 1/100 increments.
Maybe that's too fine, but that's OK.
And my initial guess for the solution when lambda equals 0
will be minus 1, and I loop over all the values of lambda,
and I use f0, which is the one-dimensional non-linear
equation solver in Matlab, to find
the roots of this polynomial function--
x cubed minus 2 lambda x plus 1.
Here's the solution and that becomes my initial guess
for my next solution.
So this loop is going to go round and round until I
get lambda equals 1, and the solution that
results there should be the solution to the initial problem
I was interested in.
Here's what that looks like graphically.
I started with lambda equals 0.
The root there serves as an initial guess for lambda
equals 1/100, and so on, and these roots
converge eventually, to something that's
close to minus 1 and 1/2.
So the actual value is a little bit bigger or a little bit
smaller than minus 1.6.
Here's one of the roots.
So rather than having to have a precise initial guess,
I transformed my problem from an easy to solve one
to a more difficult to solve one,
and I used the path along that transformation
to give me rapid convergence to better and better
initial guesses to the problem that I was interested in.
Make sense?
There's another way of doing this too.
You don't have to go from 0 to 1.
Let's go from lambda equals something big back down to 1
instead.
So if lambda's big, this term dominates
the polynomial equation.
This is the biggest term in the polynomial equation,
and it can either balance x cubed, or it can balance 1.
It can't balance against both of them, it turns out.
That
So if we have this term large, and it balances against 1, then
when lambda's large, an approximation for the root
will be 1 over 2 lambda.
If x equals 1 over 2 lambda, and these are the two largest
terms in the equation, then these two will cancel,
and f will be very close to 0.
So that's one root when lambda is very large.
If I balance these other two terms against each other,
I'll see there are two other roots, which
are plus or minus the square root of 2 lambda.
So let's start with a large lambda instead, and these
as initial guesses for the roots--
and let's trace them back to lambda equals 1.
When lambda equals 1, we recover the problem
we were interested in before.
So here, I start with three different initial guesses,
and I trace them back, and I'm actually
able to find all three roots of the initial equation.
So there's something close to minus 1 and 1/2,
something close to 1/2, and something close to 1.
So in one go, I found all three roots using continuation.
Is this clear so far?
Sam.
STUDENT: Is that one go or is it three separate [INAUDIBLE]..
PROFESSOR: That's an excellent question.
So I had to initiate this process with three
separate guesses, but it turned out in this case,
all three guesses converged to three distinct roots
in the end.
it's one go in the sense that there's
one trajectory that I had to follow,
one passive lambda as I had to go through,
but I had to solve three times to follow
each of these three paths.
There's actually no reason to stop at lambda equals 1.
I stopped there because that's the solution to the problem.
I'm interested in, but there's actually something richer
to see here.
If I keep decreasing lambda, eventually these two branches
combine with each other, and then they jump down,
and they follow the original branch that I
had that went between minus 1 at lambda equals 0,
and this root about minus 1.6 at lambda equals 1.
So the paths that these trajectories follow
can be very, very complicated, actually.
A little bit unwieldy or unpredictable.
But oftentimes, they can result in predictions of all three
roots or multiple roots associated
with the system of equations.
Yeah.
STUDENT: Going back to the first [INAUDIBLE] solve
from lambda equals [INAUDIBLE]
PROFESSOR: Zero to one?
STUDENT: Zero to one.
So then would you find only one solution there?
PROFESSOR: Yes.
So going from 0 to 1, I found one solution.
Going from infinity--
10 is big enough.
10 to 1, I was able to find three.
I didn't know beforehand that I'm
going to be able to do that.
It happened to be the case that we could
with this by construction.
This is a particular problem where it works out that way.
Other questions?
So this is neat.
Like we've done before with iterative methods,
we turn hard to solve problems into a sequence of easier
to solve problems.
Good initial guesses converge fast with Newton-Raphson.
We have a sequence of better and better initial guesses, which
all converge very quickly to eventually, the solution
that we're looking for evaluated at lambda equals 1.
So this is one way of getting good initial guesses.
You have a sequence of good initial guesses.
Make sense?
So we're just here.
We had this polynomial problem and we injected lambda
in a particular place.
We injected it there because I knew it was going to work.
I knew the problem could be transformed
in this way in order to work out and distinguish this.
It's not always obvious how to transform a problem
to make it look this way, but oftentimes we
have physical problems in which there's a parameter
that we can do continuation with that is obvious.
So think about fluid mechanics problems.
Those are nonlinear partial differential equations.
So we're trying to solve say, the Navier-Stokes equations.
Rather than jump in and try to solve that right away,
we might do continuation.
We might find a solution at very small Reynolds numbers,
where the equations are almost linear,
exact solutions are known in that case,
and we might progressively step the Reynolds number up, and use
our solution set small Reynolds numbers as initial guesses
for increasingly large Reynolds numbers until we get to the one
that we're interested in.
You can imagine doing that in mass transfer problems
where you vary the Peclet number.
In multicomponent phase equilibria problems,
this might happen by varying temperature
or pressure until we get to the temperature or pressure
that we're interested in.
If we have reaction equilibria problems,
we might vary the different reaction rates
until we get to the right balance of reaction rates
to determine the equilibrium compositions.
So they're often natural parameters in the problem
that we're free to adjust.
We can transform from easy to solve problems
to hard to solve problems in a continuous fashion.
It's not always going to work out the way that I showed.
It's not always going to be so nice and neat.
It may be quite a bit more complicated,
and sometimes there are problems for which
there isn't a natural parameter available to vary.
In those cases, there's another strategy one can employ.
It's called homotopy.
This is the transformation from one problem
into another problem in a continuous fashion.
So oftentimes, we're seeking a set of roots.
x star, which are a function of a parameter lambda,
and they're the roots of an equation that looks like this.
So they're the roots of this function, h,
which is a linear superposition of a function
f and a function g.
When lambda equals 0, h is equal to g, and so the roots of h
are the roots of g.
And when lambda equals 1, h is equal
to f, and so the roots of h are the roots of f.
So we might transform-- we make very lambda continuously
between 0 and 1 in a construction like this.
So we might transform the roots from the roots of g
into the roots of f.
Maybe the roots of g are easy to find, but the roots of f
are hard for us to find.
There's a smooth transformation though, from g to f.
When we have those sorts of smooth transformations,
there's a chance things will be well-behaved,
or at least where we encounter bad behavior--
and I'll show you what I mean by bad behavior in a minute--
there's a reliable way to interpret
what's going on there.
So we vary lambda in small increments from 0 to 1,
and the solution x star for some lambda i
is used as the initial guess to find
the roots at some lambda, which is a little bit bigger
than lambda i--
our next iterate in lambda.
Does this idea make sense?
Where do f and g come from?
So this is a big question-- f is the problem
you want to solve, typically, and g
is some auxiliary problem.
It can be arbitrarily chosen.
There's nothing to tell you why one g should
be better than another one, but oftentimes, we
choose them in a physical way.
So there will be a physical connection
between the f problem and the g problem,
and I'll show you that right now.
So here's an example.
We looked at this example before.
We want to find the roots of the van der Waals
equation of state.
It's written in terms of the reduced pressure, the reduced
molar volume, and the reduced temperature,
and let's find those roots for a given pressure and temperature.
So this is a function f of the molar volume equals 0,
and we're looking for these three roots here.
There's a root down here, which is the molar volume is low.
Not a lot of volume per number of moles, so maybe this
is a liquid like root.
Out here, the molar volume is high.
We have a lot of volume, the number of moles and materials.
This is a vapor like root, and then there's
an unstable root in between them.
So three roots to find.
We'd like to be able to find them in one go,
or we'd like to have some reliable mechanisms for finding
them.
Homotopy is one way to do this.
So let's construct a new function h of the molar
volume, which is lambda times f.
We want to find the roots of f plus 1 minus lambda times g.
And as g, let's use something physical.
So let's let g be a function that represents the ideal gas
equation of state.
So PV is nRT, but I've written everything
in terms of the reduced pressure, molar volume,
and temperature of the van der Waals equation of states.
So you pick up this factor of 8/3 there.
But this is the ideal gas equation of state.
Its roots will be the roots of the ideal gas
equation of states.
I'm going to transform from something that's easy to solve.
I know how to find the molar volume here, there's just one,
to something that's hard for us to solve.
It has three roots--
the roots of them van der Waals equation of state.
So when lambda equals 0, I'll get the roots associated
with the ideal gas, and then lambda equals 1,
I'll get the roots associated with the van der Waals
equation.
Yes?
STUDENT: If you're starting with lambda as 0
and you get the roots of g, you only get one root,
but you're trying to find two roots [INAUDIBLE]..
So how does that work?
PROFESSOR: I'll show you.
Here's what the process might look like.
I've got to start with a reduced temperature and pressure.
Here's my initial guess for the molar volume.
It's the root of our ideal gas equation.
Here's the definition of f, and the definition of g,
and the definition of h.
l times f plus 1 minus l times g,
and I'm going to loop over all of my lambdas,
solve the equation subject to some initial guess,
and update my initial guess every time I update lambda.
That's fine.
That's the code you.
Can implement the code if you want or play around with it.
So what happens if I try to carry this out?
So I'm going to vary lambda between 0 and 1.
I start with my initial ideal gas guess for the root,
and as lambda grows to 1, the molar volume
shrinks until I get to lambda equals 1,
and I've found one root of my equation.
One root, good.
There's not necessarily any reason to stop at lambda
equals 1.
That's the root I wanted, but there's
no reason I have to stop there.
So I could continue.
If I keep pushing lambda up higher and higher,
eventually I'll get to a point where all of a sudden
the molar volume jumps down to a different place.
This is the result of the algorithm.
I just told lambda to go bigger.
But there will be a discontinuity or a cusp
here that I jump off of, and I find myself
on a different solution branch, which is sort of crazy.
So now I can do something clever.
I found myself on a different solution branch.
Why not try decreasing lambda along that solution branch
instead?
Maybe I'll find a different solution,
and it turns out that's what happens.
All right, so i jump down to this other solution branch,
and now I try decreasing lambda from 2 back down to 1,
and sure enough, I found the vapor molar volume,
and now I found the liquid molar volume.
There's no reason to stop and lambda equals 1.
So if I keep decreasing lambda, eventually I'll
hit another cusp, and I'll jump back up
to the original solution branch.
So I do my homotopy.
I vary lambda from 0 to 1, and I find one root,
but if I keep going, I might identify a cusp,
and jump off that cross.
I might reverse my direction with the homotopy,
and I could find another root.
I might hit another cusp as well.
Now, when I hit a cusp, maybe it's
best not to jump off the cusp, but instead
to reverse direction at the cusp,
and try to trace it out the other direction.
And if I do that, if I get right to this last point,
and I try to reverse direction again,
I can trace out one more solution path,
and I can find the third root to the equation.
So in one go, if I'm paying attention
to how my roots are changing-- if they change by too much,
I might switch the direction with which I change lambda,
and I could find all of these roots at once.
So three roots in one go.
Yes.
STUDENT: Is this like trial and error?
How do you know when to change?
PROFESSOR: This is a wonderful question.
So in this case, you can do it by trial and error.
You can detect these jumps.
That's even better.
You look for jumps in the value of the solution that
are of sufficient magnitude, and when you detect them,
you can reverse the direction.
There's an even more systematic way
to do this, which is called arclength continuation.
Give me one second, and I'll give you
slide on arclength continuation, which
is a systematic way of doing exactly this process of tracing
out this path in the solution plane.
Roots versus homotopy parameter.
These cusps here, are sometimes called turning points,
and they have a special property,
which is the determinant of the Jacobian of the homotopy--
the homotopy function-- is going to be equal to 0.
The Jacobian is going to be singular right at these cusps.
So you can detect them by tracking what
the Jacobian is doing right.
The Jacobian is going to be non-singular at most
of these other points, but when you
hit one of these turning points, the Jacobian will be singular.
There isn't a well-defined direction for the homotopy
to proceed along.
It hits this cusp where the slope
is infinite in this plane.
There isn't a well-defined direction to step along here.
The question was--
OK, it seems a little ad-hoc.
You've got this 2D plane, and it's
easy to guess when you jump from one solution to another,
and trace the path back.
That's fine.
There's a systematic way of doing
this, which is called arclength continuation, which
says I know that there's some curve in this solution plane.
There is some curve here.
I've shown it to you graphically.
Let's try to parametrize that curve in terms
of a measure of the arclength.
Call it s, the distance along the curve.
So the roots are a function of lambda,
which can be a function of the distance along the curve,
and lambda is a function of the distance along the curve,
and there is a constraint which says
s has to satisfy an arclength formula.
This is the formula for the arclength in this plane
of roots versus lambda.
The derivative of this function with respect to s
squared plus the derivative of lambda with respect to s
squared equals 1.
That defines s as a measure of length along this curve.
You learned arclength formulas in the context of integration
at one point.
This is the differential version of that same formula.
So I've got a differential equation.
Well, it's actually a differential algebraic
equation.
This is an algebraic constraint that
needs to be satisfied for all values of s,
and it depends on how the roots change with s.
We'll learn how to solve differential
algebraic equations in the two sections from now.
So we don't quite know how to solve something like this yet,
but if we know how to solve those sorts of equations,
then it's clear that we can trace out
this entire curve in the solution plane,
and then go back and try to find the places where the curve
intersects lambda equals 1.
So if we follow this path by solving this equation from s
equals 0 out to some very large s,
then we should be able to find all of the roots that
are connected to the path.
Yes.
STUDENT: If there's a turning point that [INAUDIBLE] then
does that imply multiplicity?
PROFESSOR: Yes.
So it's a real problem.
That turning point singularity is actually
built into this arclength formula as well.
This is the derivative of x with respect
to s squared, which takes the direction you're
moving in out of the problem.
This is a length, not a direction.
So when we get to these turning points,
the solution methods for these equations
are going to have to know that I was headed
in a downward direction first.
I can get to this turning point, and then the question
is which direction do I move in from the turning point?
Do I go up or do I go down?
Which way does the solution path change?
The slope isn't defined there, so it's quite challenging
to figure out exactly how to proceed.
But it is what you say.
It's like the solution has multiplicity at that point.
Curves are tangent in the solution plane.
It's a great question.
OK, so that's the notion of arclength continuation.
We may talk about this later when
we talk about differential algebraic equations,
and we'll talk about some equivalent sets of equations
that we know how to solve, but they
have to satisfy this constraint that s
is a measure of arclength along the curve.
But there you go.
You can find all the solutions, at least all of the solutions
connected to this path defined by the homotopy equation.
You can find all of them in one go using
some clever methodologies.
Oftentimes, we're not just worried about these
turning points or cusps, we encounter bifurcations
in the solution.
So we may have a path that has one root along it,
and then that path may split, and there
may be more than one root.
Here's a simple example.
Find the real roots of x cubed minus rx.
If r is negative, there will always be one real root,
and it'll be x equals 0.
If r is equal to 0, there will still be one real root, it's 0.
That root will have multiplicity 3 in that case.
And if r becomes positive, I'll suddenly
have three real roots instead of one.
So as I vary the parameter r from a negative number
to a positive number, I'll go from one real root to three.
In the solution plane, that will look like this.
Here are the roots as a function of r.
r is less than 0.
All the roots are equal to 0.
When I hit r equals 0, I have multiplicity
associated with the root, and the solution bifurcates,
and there will be three possible solutions now.
If I was trying to do a homotopy problem where I very r
continuously from some negative number to some positive number,
I'll hit this point, and I'll want
to follow all three paths to find solutions out here
at positive r.
But I can detect those points.
Just like I was said before, at these weird turning points
or bifurcations where we have multiplicity in the roots,
the Jacobian in is going to be singular.
The Jacobian of the homotopy is going to be singular.
Its determinant will be equal to 0.
So I can detect these points by checking the determinant
of the Jacobian.
And when I find one, I want to do
something smart to try to follow these paths.
Your homework this week has a problem like that.
You're trying to find all the azeotropes of this ternary
vapor liquid equilibrium system, and you'll
have a homotopy parameter that you change,
and it'll bifurcate.
OK, so as you change the parameter,
you'll go from one path to two, and then those two paths
will split again.
You'll be splitting from a single component solution
into a two component azeotrope, into a three component
azeotrope.
And you'll want to detect where these bifurcations occur,
and try to follow solution branches off
of those bifurcations, and there's
some good hints in the problem for how to do that.
So occasionally, we'll encounter problems
that switch from having one solution to having
many solutions.
You've seen how this happens in the discontinuous sense
with these turning points.
Here, for a given value of lambda,
all of a sudden it's clear that not only is there one solution
here, but there's another solution down there.
My path has to choose.
Do I want this solution or the one below?
I jump off this turning point.
That's a sort of bifurcation.
These are continuous bifurcations.
They're often easier to take care of.
So the bifurcations in a homotopy
let us find multiple solutions.
Every time we detect a bifurcation,
if we just track those solution branches,
we split our algorithm up to track each of the branches
separately, they'll terminate at different roots
to the homotopy equation.
And when we get to the value of lambda equals 1
in our homotopy, we're there.
We found multiple roots to the original equation.
These are often of great physical interest.
Do you know what sort of things these bifurcations represent?
One occurs in the van der Waals equation.
So if I start at very high temperatures and a given
pressure, the van der Waals equation
will have how many real roots?
One-- a gas.
And then as I turn the temperature down,
there is a propensity for liquefaction of that gas phase.
The formation of a liquid in coexistence with the vapor.
So I'll go from having one root to having three roots
all of a sudden, and there will be this bifurcation
point where that transition happens
in a continuous fashion.
STUDENT: Can I ask you a question?
PROFESSOR: Yes.
STUDENT: [INAUDIBLE]
PROFESSOR: Yeah.
STUDENT: So as you early put it out that the turning point,
the determinant of the Jacobian [INAUDIBLE] including lambda
is singular.
But here you're writing j of x.
Is that [INAUDIBLE]
PROFESSOR: Excuse me.
This is a typo.
So there should be a little h here.
This is the determinant of the Jacobian
of the homotopy equation at that root.
Where there's a bifurcation, the Jacobian will be singular.
I apologize.
That's a little bit--
STUDENT: And that lambda.
PROFESSOR: And that lambda.
STUDENT: Just to make it clear, you have n state variables.
Like x is how many unknowns you have,
and then you're adding another unknown to lambda,
another parameter.
This Jacobian has n plus 1 dimension.
This is also [INAUDIBLE] Is this correct?
PROFESSOR: Well, no.
Actually, the Jacobian with respect to x only
is also singular at these turning and bifurcation points.
STUDENT: [INAUDIBLE] as well.
PROFESSOR: Yes.
It's also true what you said, that if I have the Jacobian,
and I take the derivatives with respect
to lambda, that will also be singular,
but it will be singular because the subspace, which
is the derivatives with respect to x,
gives you a singular component.
The condition is the same as here.
So it's the Jacobian of h, taking the derivatives
with respect to x--
all the different x's of h.
The determinant of that matrix at a particular lambda
will be equal to 0, and the value of lambda
will be the lambda at which you have this turning point,
or which you have a bifurcation.
Here's what it looks like in a two-dimensional sense.
So I'll have my homotopy equation
where I'm trying to find the roots of h.
Say it's a two-dimensional-- x is dimensions.
So it has x1 and x2, and these curves
represent h1 equals 0 and h2 equals 0,
and they intersect at one point for a given value of lambda.
There's the root.
Now I increase the value of lambda.
The homotopy equations change.
So these curves change shape, and all of a sudden I
go from crossing curves to tangent curves,
and we know tangent curves in the plane
will have a Jacobian which is singular.
And this is the bifurcation point.
These two curves become tangent, and as I
go to a bigger value of lambda, they'll change shape again,
and I'll now have multiple roots.
I'll bifurcate from one solution to two,
or one solution to three, but the bifurcation
will occur through this transition
where I have tangent curves, and we know there,
the determinant of the Jacobian is 0.
There's a singularity under those conditions.
Does that makes sense?
OK, so a couple of things about the practice, and I'll
give you a simple example you can try out, which I think
is an interesting one, if pathological.
So in practice, it's hard to hit this bifurcation point exactly.
And I'm telling you that you should
be detecting these branches.
You should follow these branches,
but clearly, I've got to know where this point is in order
to follow the branches.
That's a problem.
If the determinant of the Jacobian
is 0 at the bifurcation point, then it's
going to change from positive to negative as I
cross the bifurcation point with lambda.
Actually, I don't necessarily have
to find the point where the determinant of the Jacobian
is equal to 0.
What I really should do is try to track
the sign of the determinant of the Jacobian,
and see when it goes from positive
signed to negative signed.
And when it does that, it will have
crossed through one of these bifurcation points
or one of these cusps.
That's the point at which I should stop,
and say, OK, there's possibly other solution branches
to follow.
So as I step from my smaller lambda to my bigger lambda,
the sign changes.
I want to step in some different directions
than the direction I was heading before,
and those different directions will trace out
the other solution branches.
It turns out those different directions
belong to the null space of the Jacobian.
So if you find vectors that are very
close to the null space of the Jacobian,
those directions are the directions
you should go to find these other solution branches.
So it's very interesting.
Yeah.
STUDENT: Can you have a Jacobian that goes from positive down
to 0 and [INAUDIBLE]?
PROFESSOR: That may be possible.
You won't encounter that in the problem that you're doing.
That would be a very difficult problem
to do branch detection on.
Let's Suppose that was the case and you
want to find exactly where that bifurcation point is.
These bifurcation points are of physical significance.
In the van der Waals equation of state,
that bifurcation point is the critical point.
It's the critical temperature at which that phase separates.
Maybe I want to know that exactly.
Well, all I need to do is solve an augmented system
of equations.
So let's solve a system of equations,
which is our original hematology equation equal to 0,
and determinant of the Jacobian of this homotopy equation
equal to 0.
This is if h had dimension n.
This is now n plus 1 equations, and let's
solve it for n plus 1 unknowns.
Let's solve it for x and for lambda. x has dimension n,
and we add one more unknown--
the value of the homotopy parameter at which
this bifurcation occurred.
So we can find that point exactly.
So there may be cases where this happens
where it doesn't continuously transition
from positive to negative.
It may hit 0 and go back up, but that's
a really challenging case to handle, but you can handle it.
You can find precisely the point at which
this bifurcation occurs by solving
these augmented equations.
It's just another non-linear equation to solve.
Yes.
STUDENT: I'm not sure about the previous slide.
So in the previous slide, why is the [INAUDIBLE]
PROFESSOR: Let me go back to here.
So let's let the blue line here, be
the place where the first element of h is equal to 0.
So this is all the values of x at which
the first element of our homotopy top equation
is equal to 0.
And the red curve here, let's let
that be all the values of x at which the second element of h
is equal to 0.
And so h itself, is only equal in 0 where both of these curves
are equal to 0, and that's their intersection.
So they intersect here, and this star is the solution.
It is the root of this non-linear equation up here.
Now, I increase lambda.
I make lambda a little bit bigger.
That changes h.
h is a function of lambda.
So h changes, and these curves change shape in the plane.
And at that bifurcation point, these two curves
will become tangent to each other.
There will be one root, but in this solution plane,
the two curves will be tangent.
And we saw before that when we get
these sort of tangency conditions
or the equivalent in higher dimensions,
the Jacobian of our function we're
trying to find the roots of will be singular.
It's like a multiple root.
And when we increase lambda further,
curves change shape again.
I'm not changing the red one because that's complicated.
I just changed the blue one for you
as a graphical illustration.
So the blue one changes shape, and now I've
gone from having one solution to having three.
So in order to go from one solution to three,
in a continuous fashion, I have to go
through this tangency state, and that defines the bifurcation
point.
Does that make sense?
Here's a simple example you should try out.
So I'll present it to you, and then you
should try to solve it in Matlab.
You actually know the solution, geometrically,
so it's not too hard.
So here's a function.
The first element to that function,
it's a function of two unknowns-- x1 and x2.
The first element to that function
is an equation for a circle of radius r with center at minus 3
minus 1, and the second element is
an equation for a circle of radius r with center at 2, 2.
So there's two circles.
We want to find the radius where the circles just touch,
and that point is a bifurcation point.
If r is perfect, these two just touch.
If r is a little too big, the two circles overlap,
and there are two solutions to this problem.
And if r is a little too small, there
are no solutions to this problem.
So the perfect r for touching is a bifurcation point.
It has all the properties of a bifurcation point,
just like we discussed.
It's a multidimensional equation and it bifurcates.
You know the solution to this problem, geometrically.
You could work it out in a couple of minutes,
but you could also solve it using
these augmented equations.
So we know that point at which the circles just touch,
the value of r at which those circles just touch,
will be the place where f of x is equal to 0
at which these equations define a circle,
and at which the Jacobian of these equations--
the derivatives of f with respect to x, it's determinant
is also equal to 0.
So we want to solve this system of three equations for three
unknowns.
The values of x1 and x2, where the circles just touch,
and the value of r, the radius of these circles
at which they just touch.
That would find the bifurcation point.
That will find the touching of the circles.
Can you see that, geometrically?
All these bifurcation problems work exactly that way.
There's some homotopy parameter.
In this case, that's r.
But there's some parameter we're varying as we
pass through that bifurcation.
We'd like to know the values of our unknowns
and the value of this parameter at that bifurcation point
or at that critical point.
So you find this point by solving
the augmented equations.
We have a non-linear equation-- three nonlinear equations
for three unknowns, and you can do
that with a Newton-Raphson iteration.
It's a non-linear equation.
We know Newton-Raphson is the way to solve these things.
To do Newton-Raphson, that means we need to compute
the Jacobian--
this is confusing.
We've got to compute the Jacobian
of these augmented equations, the derivatives
of these augmented equations with respect
to our augmented unknowns.
So that will be the derivatives of f with respect
to x, the derivatives of f with respect
to r, the derivatives of the determinant of the Jacobian
with respect to x.
That's this gradient here.
The derivatives of the determinant of the Jacobian
with respect to r.
This is our augmented Jacobian.
The augmented Jacobian times the Newton-Raphson step.
That's equal to minus the function
we're trying to solve evaluated at the previous iteration.
Solve the system of equations, you'll
get your Newton-Raphson step, and you can find the solution
to the augmented equations.
It'll converge quadratically.
You should commit yourself, you can figure this out.
Don't try to compute all of these derivatives.
That's complicated.
You're going to make a mistake.
I'm going to make a mistake if I try to do that.
What do we do instead to compute the Jacobian?
How do we compute these derivatives?
What's that?
Finite difference, right.
Don't play games with it.
Try to use finite difference to compute the elements of this
instead.
You know the exact solution too, so you
can use that as an initial guess,
or something close to the exact solution as an initial guess
to see if you can find this bifurcation point.
So homotopy bifurcation, there are
ways to get good initial guesses by continuously
changing our problem.
We get a succession of problems that
will converge faster because we have good initial guesses,
and that will give us a reliable convergence, hopefully, to one
of the roots of the equation we're after.
If we're able to pick up solution branches, we can track
these cusps, or find these bifurcation points,
we might be able to find multiple solutions.
There are problems for which depending
on your choice of g in this homotopy function,
maybe some of the solution branches
are disconnected from each other.
Boy, that's crazy and complicated.
We're not going to worry about those sorts of cases,
but they can pop up and happen.
But the point is there are a robust strategies
you can employ to find multiple roots
of these equations or multiple minima in optimization problems
using precisely this technique.
You'll have a chance to practice this on this week's homework
assignment, and I think you should
take a look at this problem.
It's simple.
It looks complicated , but if you can understand this from,
then you'll understand everything you need to know
about bifurcation and homotopy.
Any questions?
Yes.
STUDENT: At a bifurcated point, we
go from having one root to three roots, per se.
And you said that when that occurs,
we need to follow each of these solutions, right?
PROFESSOR: Yes.
STUDENT: And you said that we need to step
in a different direction.
So what does that exactly mean?
In a 2D sense, we have this stepping forward
or stepping backwards.
PROFESSOR: This is a wonderful question.
So let me go backwards here.
Here's the bifurcation point, and when
we go through that bifurcation point,
there are now going to be three solutions.
And I'm going to need some initial conditions
for my solver to find each of these three solutions.
I can't give them all the same initial condition.
I'm going to try to track all three of these branches
with different initial conditions.
If I give them all the same initial condition,
and my algorithm meets the minimum threshold for what
an algorithm should be, they're all
going to converge to the same root.
So that's useless.
So I need to give them different initial conditions.
So what are those supposed to be?
In your problem, I'll give you explicit directions
for how to find those.
It turns out those initial guesses
should be perturbed in directions
that belong to the null space of the Jacobian.
So I want to perturb from that one solution, which was here,
in a direction that's parallel to the tangency
of those curves.
Those two curves were tangent with each other.
In the other solutions, they're going
to sit off in the direction that's tangent to those curves.
So I've got to find the vector that's in the null space,
or very, very close to the null space,
and that will tell me how to get better initial guesses
for picking up these other solution branches.
It's not trivial.
It's actually quite challenging.
You'll have a chance to practice.
It's part five of five parts on this problem.
The first four, they're relatively
straightforward and easy.
Part five is a little bit tricky.
You get explicit directions, and the TAs
know that that's the most difficult part of that problem.
So we'll help you through it.
It's something interesting to think about trying to do.
You may find you're only able to pick up some of these solution
branches, but not all of them.
That's OK.
There are ways to pick them all up if your algorithm is
designed well enough.
You may just pick up some of them,
and that may be the result you turn in.
You can look at the solution to figure out
how the other ones are picked up as well.
But you're looking at the direction
along these tangent curves.
That's a vector in the null space of the Jacobian.
And if you perturb in those directions,
then you can trace out these other roots
so you can find the bifurcation.
It's a wonderful question.
You'll get a chance to practice on it.
OK, last question.
STUDENT: [INAUDIBLE]
PROFESSOR: Yes.
STUDENT: Is this supposed to be jrx?
PROFESSOR: That's the Jacobian of f.
That's the derivatives of f with respect to x.
STUDENT: Not including r?
PROFESSOR: Yes.
Not including derivatives with respect to r.
Just derivatives with respect to x.
So this is the Jacobian of f.
The derivatives of f with respect to x.
The determinant of that will be 0 at the bifurcation point.
It'll turn out that--
this will be the last thing and then
I have to get off the stage here.
It'll turn out that this matrix is also singular
at the bifurcation point.
It will be singular because the subspace j is singular.
So it'll turn out to be singular as well.
Don't worry about it.
The condition for the bifurcation point
is that the determine of the Jacobian is equal to 0,
and that point is a root of the equations you're looking at.
OK, thank you.