# How to solve equations self-consistently

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.

-
You can use FixedPoint. – b.gatessucks Feb 26 '13 at 9:55
Specifically, FixedPoint[Total[Eigenvalues[{{1, #}, {#, 1}}]] &, d] possibly with a SameTest setting to account for your dmin. – Daniel Lichtblau Feb 26 '13 at 14:34

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.
-
Thanks for the answer, but I indeed need to iterate since what I want to do, is to solve the so called Bogoliubov deGene equations self-consistently(this is related to Superconductivity), so I do need to iterate until d keeps more or less the same. My steps must be: Diagonalization->define d->Diagonalization->redefine d->Diagonalization->redefine d->...until the change in d is smaller than a certain delta. Thanks again. – Mencia Feb 26 '13 at 12:14
Well, I'm an economist, so I've no idea what the Bogoliubov deGene equations are. :) But I'll update my answer with a toy example. – Verbeia Feb 26 '13 at 12:17
Thanks for the answer, but I indeed need to iterate since what I want to do, is to solve the so called Bogoliubov deGene equations self-consistently(this is related to Superconductivity), so I do need to iterate until d keeps more or less the same. My steps must be: Diagonalization->define d->Diagonalization->redefine d->Diagonalization->redefine d->...until the change in d is smaller than a certain delta. Thanks again. – Mencia Feb 26 '13 at 12:19
Oh! And by the way Verbeia, I am not trying to be rude, but I do not know wether economist know what diagonalization means. It means just to get the eigenvectors(u) and eigenvalues(e) of a matrix(H). I am sure you know it ;). Thanks again. – Mencia Feb 26 '13 at 12:23
No worries, we call it something different. I'm not sure I fully understood what you are trying to do, but it should be enough to get you started. – Verbeia Feb 26 '13 at 12:33