The Ising model is a well-known and well-studied model of magnetism. Ising solved the model in one dimension in 1925. In 1944, Onsager obtained the exact free energy of the two-dimensional (2D) model in zero field and, in 1952, Yang presented a computation of the spontaneous magnetization. But, the three-dimensional (3D) model has withstood challenges and remains, to this date, an outstanding unsolved problem.
There is a result I only heard about recently: it has been proven that computing partition functions for the Ising-model in dimensions > 2 is NP-complete. (The paper can be found at http://www.cs.brown.edu/people/sorin/pdfs/Ising-paper.pdf; a more readable one is here http://www.siam.org/pdf/news/654.pdf - both can be found on the Wikipedia on the Ising model). I'm far from an expert on this, but the main idea is that a certain NP-complete graph theory problem on finding maximal sets of edges can be mapped to ground states of Ising-3D. Roughly, this means that you can't find ground states in polynomial time, and as most physicists know, if the difficulty of your problem scales exponentially, solving something exactly for large systems quickly becomes impossible.
The 3d Ising model can be "solved" in a certain sense, it can be recast as the problem of a discrete lattice Fermionic string. This method is explained in detail in the last chapter of Polyakov's "Gauge Fields and Strings", and it is the natural generalization of Onsager's method to 3d.
This method does not produce analytical computable critical exponents in 3d as yet, but not because the 3d model is intractable. The proofs that you have intractability in calculating the free energy for an arbitrary sublattice of the 3d model is interesting, but it also works in 5d or 6d, where the critical exponents are mean field, and so exactly computable. This proof only shows that the general solution, in the sense of computing all correlation functions in the presence of arbitrary external fields, is going to be difficult. But it does not mean that the 2 point function is uncomputable in the long-distance limit.
The only precise meaning I can see to the statement that a statistical model is solvable is saying that the computation of the correlation functions can be reduced in complexity from doing a full Monte-Carlo simulation. In this regard, knowing that the configurations of the 3d model are described by Polyakov lattice Fermionic strings does help, because you can simulate noninteracting strings enclosing volumes instead of spins on each site. The issue is that the strings are Fermionic, so it might not be possible to actually simulate a typical configuration using Polyakov's transformation any more simply than the usual way, because of the Fermion sign problem.
This is all investigated by Polyakov from time to time, and there is still a reasonable hope for a new idea which will lead to progress, the computational intractability results nonwithstanding.
Two-dimensional theories simply have much more mathematical structure that makes many such models mathematically solvable - integrable.
In particular, in the long-distance limit, one obtains a scale-invariant theory that is typically conformally invariant as well. Two-dimensional conformal symmetry is - unlikely any higher-dimensional symmetry - infinite-dimensional. This fact plays a very important role in string theory which has 2-dimensional world sheets. In some sense, the 2D systems lead to "infinitely many conserved quantities" which often makes their physics solvable.
Perturbative string theory reflects much of the special mathematics that makes problems tractable in two dimensions. Needless to say, 1D systems may be as solvable as 2D systems or more so. Integrable - analytically solvable - systems have also included spin chains. All these things are parts of string theory in one way or another. On the other hand, problems in three or more dimensions are qualitatively harder and most of the questions about the Ising model in 3D and similar models are not analytically solvable.
Another question is whether one may understand a model qualitatively. Of course, an analytic understanding gives one a superior tool to answer this question as well. When it doesn't, it's still possible to gain some qualitative understanding - numerically or by various approximation schemes - and the fact that as of 2011, it hasn't been done, is just a historical accident that is more likely because it's a difficult problem.
Exaxt solvability has nothing to do with NP-completeness.
For equations on a lattice or a continuum, exact solvability happens to be equivalent with having enough symmetries to allow the solution to be determined by exploiting these. (To a large extent, this even holds for ordinary differential equations in more than a few variables.)
The reason that a few (classical or quantum) systems are integrable therefore comes from the fact that they have a much larger (infinite-dimensional) symmetry group, and hence infinitely many conservation laws, while a typical system has only a small, low-dimensional symmetry group. This is the (modern) explanation why Onsager's solution works, while there is no analogous solution in the next dimension.
If one looks at lists of integrable systems (e.g., the one at http://en.wikipedia.org/wiki/Integrable_system#Exactly_solvable_models which for the classical case seems fairly complete) one sees that they get very scarce in higher dimensions. There are just not enough possible large symmetry groups around....
I solved the Ising model:
I. A. Stepanov. Exact Solutions of the One-Dimensional, Two-Dimensional, and Three-Dimensional Ising Models. – Nano Science and Nano Technology: An Indian Journal. 2012. Vol. 6. No 3. 118 - 122. (The paper is on the Journal site with a free access)