How to solve equations self-consistently

Question

I want to solve the following equation self-consistently. So,

H.u = e.u

{{1, d}, {d, 1}}.{u1, u2} = e.{u1, u2}

1. I guess an initial value for d
2. Diagonalize H and get eigenvectors and eigenfunction u and e->Eigenvector[H] Eigenvalue[H]
3. Define d -> d = u1 + u2
4. Plug d in the matrix again
5. Diagonalize H... once and again until I get certain convergence, i.e., until the value I get for d does not change within a certain convergence range dmin.

Of course, this a really simplified kind of problem of my actual problem I have to solve with a huge matrix.

This is really easy to do in Fortran90, but I wonder how do you do this kind of things in Mathematica.

Answer 1

Yes, you can use FixedPoint or possibly NSolve or FindRoot with one of its methods, if you are willing to specify the values of the non-numeric parameters here.

But be aware that you need to specify the equations to be solved appropriately. First you need Equal (==, a double equals sign) not Set (=, single equals sign) to specify your equations. Second you need to set your equations up so that Mathematica recognises them as actual equations. Consider:

{{1, d}, {d, 1}}.{u1, u2} == e.{u1, u2}

{u1 + d u2, d u1 + u2} == e.{u1, u2}

This is not a list of equations, and Dot doesn't work with a mix of vectors and scalars. Now try the following.

Thread[{{1, d}, {d, 1}}.{u1, u2} == e {u1, u2}]

{u1 + d u2 == e u1, d u1 + u2 == e u2}

Solve[Thread[{{1, d}, {d, 1}}.{u1, u2} == e {u1, u2}], {u1, u2}]

{{u1 -> 0, u2 -> 0}}

 Solve[Thread[{{1, d}, {d, 1}}.{u1, u2} == e {u1, u2}], {e, d}]

{{e -> 1, d -> 0}}

Your real problem is of course bigger than this but it seems to me that if it is linear, you don't necessarily need to iterate to the solution unless you are doing a numerical methods course and you have to complete an assignment that proves it can be done that way.

If you were trying to solve a $m.x=b$ type structure for $x$, then you would use LinearSolve.

Edit to cover FixedPoint

You are right that different disciplines use different terminology to get the same thing. Here is an example using FixedPointList, which shows the intermediate iterations for your information. If you don't need these, just use FixedPoint with the exact same syntax.

If I understood the list in your question properly, you want the vector u to be the eigenvector, but this is a matrix (list of eigenvectors) not a vector. Perhaps you want the first or last ones.

FixedPointList[Plus @@ Last@Eigenvectors[{{1, #}, {#, 1}}] &, -0.2, SameTest->delta]

To explain this code:

• Last@Eigenvectors[...] gives the last eigenvector, which I assume is what you want for u. (Maybe First instead, it depends on what you're looking for.)
• Plus @@ takes that list representing u and sums it. You could also use Total. You did say that the next iteration of d was u1 + u2.
• Notice the use of Slot (#) and pure functions to define how d is updated.
• The second argument is the starting value for d, and the third is the convergence test. Obviously delta actually needs to be a small positive number.