Iterative Methods and Radiosity

Section 1: Building a Radiosity Matrix

In this demonstration, we do a simple radiosity calculation using various iterative methods to see how the results change over time. The basic idea is that we have a couple of spheres, one of which is glowing, and we study the distribution of light across the scene. We do this by breaking the spheres up into patches and solving for the radiation and reradiation of light between patches. The transfer of light energy between two patches of surface area is controlled by the radiosity equation:

(Subscript[B, i]) = (Subscript[E, i]) + (Subscript[R, i]) ∑(Subscript[F, ij])(Subscript[B, j])

where (Subscript[B, i]) is the radiosity of patch i, (Subscript[E, i]) is the light emitted by patch i, and (Subscript[F, ij]) is a ``view factor’’ describing the relative geometry of patch i and patch j.

We start with some code to build a sphere at a given radius and location as a polygonal graphicscomplex. We end up with a list in the form {Center, Radius, <List of 3d coordinates on faces in 3d coords>, <List of unit outward normals of faces>, <List of centers of Faces in 3d coords>, <List of Areas of Faces>,<GraphicsComplex displaying the whole thing>}.

Needs["VectorAnalysis`"];

PolygonArea[P_] := Module[{N,NewP,A},
NewP = Append[P[[1]],P[[1,1]]];
N = Normalize[CrossProduct[NewP[[3]] - NewP[[2]],NewP[[1]]- NewP[[2]]]];
A =(1/2)*Sum[ScalarTripleProduct[NewP[[i]] - NewP[[1]],NewP[[i+1]] - NewP[[1]],N],{i,2,Length[NewP]-1}];
A
];

BuildSphere[Center_,Radius_] :=Module[{CR,GC,FaceTable,NormalTable,SpaceCenters, FaceAreas,SpaceFT,SpaceVerts,SpaceGC},

CR = N[PolyhedronData["Icosidodecahedron","Circumradius"]];
GC = N[PolyhedronData["Icosidodecahedron" ,"Faces"]];
FaceTable = Table[GC[[1]][[GC[[2,1,i]]]],{i,1,Length[GC[[2,1]]]}];
NormalTable = Normalize /@ (Mean /@ FaceTable);
SpaceCenters = Map[(Radius/CR)*# + Center &, NormalTable,{1}];
SpaceFT = Map[(Radius/CR)*# + Center &,FaceTable,{2}];
SpaceVerts = Map[(Radius/CR)*# + Center &, GC[[1]],{1}];
SpaceGC = Normal[GraphicsComplex[SpaceVerts,GC[[2]]]];
FaceAreas = Map[PolygonArea,SpaceGC[[1]]];
{Center,Radius,SpaceFT,NormalTable,SpaceCenters,FaceAreas,SpaceGC}

];

Show[Graphics3D[{White,BuildSphere[{2,0,0},5][[7]]}],Lighting->"Neutral"]

GC = BuildSphere[{2,0,0},5][[7]];

We now test the face area calculations. Unfortunately, Mathematica doesn’t give us the areas for the faces, nor a function for the area of a polygon (really, Wolfram?), but we do get the surface area of the entire thing. Since we know that there are 32 faces, we can compute our surface as 32 * the mean of the face areas and subtract the stored SurfaceArea data to compare. We should get 0.

32*Mean[BuildSphere[{2,0,0},5][[6]]] - (5/PolyhedronData["Icosidodecahedron","Circumradius"])^2*PolyhedronData["Icosidodecahedron","SurfaceArea"]

We can add colors to each of the polygons by enclosing each polygon in a list in the form {color,Polygon[ ]}. Here’s a little function to take a list of colors and apply them to the faces of our polyhedron. We wrap that into a larger DisplaySolution function which take a list of RadiosityLevels and applies them to the faces of a collection of spheres. We keep in mind that the overall radiosity solution gives the TOTAL energy emitted by each face, so we must divide by face area in order to get the glow color of the face in DisplaySolution.

WarmColor[x_] := Module[{xScaled},
xScaled = Min[1.0,x];
RGBColor[xScaled * 0.9777676050965133`,xScaled*1.`,xScaled*0.04327458609903105`]
];

ApplyColors[Colors_,GC_] := Module[{ColorList,FaceColors,PadList},
ColorList = Map[Glow,Colors];
PadList = PadRight[ColorList,Length[GC[[1]]],Colors];
FaceColors = Table[GrayLevel[0.8],{Length[GC[[1]]]}];
Transpose[{PadList,FaceColors,GC[[1]]}]
];

DisplaySolution[RadiosityLevels_,Spheres_] := Module[{RPS,GC,i},
RPS = Partition[RadiosityLevels,32];
GC = Table[ApplyColors[Table[WarmColor[0.02+(1/Spheres[[i,6,j]])*RPS[[i,j]]],{j,1,32}],Spheres[[i,7]]],{i,1,Length[Spheres]}]
];

Now we can display the results. We need to turn lighting off so that we just see the color we applied to each face.

ApplyColors[Table[GrayLevel[i/32],{i,1,Length[GC[[1]]]}],GC];

Graphics3D[ApplyColors[Table[GrayLevel[i/32],{i,1,Length[GC[[1]]]}],GC],Lighting->"Neutral",Boxed->False]

Now we can start computing view factors. This is made more complicated by the fact that some polygons are going to be blocked by others in the scene. The first simplifying observation that we can make is that all the bodies in the scene are convex, so if the line joining the centers of the faces makes an obtuse angle with the (outward pointing) normal, the form factor should be zero. After that, we do a more time-consuming test which detects an occluding sphere in between the faces in question. We will operate on a “face” format which stores the center, normal, and area of the face.

FacingQ[{{A_,NormalA_,AreaA_,indA_},{B_,NormalB_,AreaB_,indB_}}] := Dot[NormalA,B-A] > 0 && Dot[NormalB,A-B] > 0;

ParallelNormal[A_,B_,{Center_,Radius_}] :={Dot[Normalize[B-A],Center -A], 
Norm[(Center - A) - Dot[Normalize[B-A],Center-A]*Normalize[B-A]] - Radius};

ClosestQ[{{A_,NormalA_,AreaA_,indA_},{B_,NormalB_,AreaB_,indB_}},ListOfSpheres_] := Module[{PN,SphereIntersectsAB,CenteredSpheres},
PN = Map[ParallelNormal[A,B,#]&,ListOfSpheres];
SphereIntersectsAB = Select[PN,#[[2]] < 0&];
Min[SphereIntersectsAB[[All,1]]] > Norm[A-B]
];

FacingQ[{{{1,0,0},{1,0,0},1,1},{{5,0,0},{1,0,0},1,2}}]

This code deserves a little explanation. We first note that by convexity, for each sphere that intersects the line segment AB, there is a perpendicular from the line to the center of the sphere inside the intersection. The closest intersecting sphere has this intersection point closest to A. We start by computing the parallel/normal function for each sphere-- we take the vector from the center of the sphere to A, resolve it into parallel and normal components to the vector B - A, and record the length of the parallel component and the difference between the radius of the sphere and the normal component (if this is negative, the sphere intersects the line).

The minimum parallel distance among the intersecting spheres is the closest sphere. If that sphere is closer than B, we return FALSE. Otherwise, we return TRUE.

For example:

B = {1.5,0,0};
A = {1,0,0};
ListOfSpheres = {{{2,4,0},5},{{4,7,0},3},{{10,3,0},7}};
ClosestQ[{A,{1,1,1},2},{B,{1,1,1},4},ListOfSpheres]

B = {5,0,0};
A = {1,0,0};
ListOfSpheres = {{{2,4,0},5},{{4,7,0},3},{{10,3,0},7}};
ClosestQ[{A,{1,1,1},2},{B,{1,1,1},4},ListOfSpheres]

At this point, we basically want to take all pairs of faces for which FacingQ and ClosestQ are both TRUE and compute

ViewFactor[{{A_,NormalA_,AreaA_,indA_},{B_,NormalB_,AreaB_,indB_}}] := {indA,indB} -> -0.9*AreaA*AreaB*DotProduct[NormalA,B - A]*DotProduct[NormalB,A - B]/(Pi*(Norm[A-B]^4));

We start by computing a block of the view factor matrix. We will take two distinguished spheres in the BuildSphere format and a list of the other spheres in the scene and compute a block of the ViewFactor matrix corresponding to the pairs of faces in these spheres.

VMBlock[{SphereA_,SphereB_},OtherSpheres_] := Module[{SList,AFaces,BFaces,VFMatrix,Tups},

If [Norm[SphereA[[1]] - SphereB[[1]]] < 0.01,SparseArray[IdentityMatrix[32]], 

AFaces = Transpose[{SphereA[[5]],SphereA[[4]],SphereA[[6]],Range[32]}];
BFaces = Transpose[{SphereB[[5]],SphereB[[4]],SphereB[[6]],Range[32]}];
SList = OtherSpheres[[All,1;;2]];

Tups = Select[Tuples[{AFaces,BFaces}],FacingQ[#] && ClosestQ[#,SList] &];
SparseArray[Map[ViewFactor,Tups],{32,32}]

]

];

VMBlock[{BuildSphere[{0,0,0},5],BuildSphere[{10,0,0},3]},{BuildSphere[{5,0,0},1],BuildSphere[{20,5,7},3]}]

We now need to assemble the blocks into the large matrix. We can keep a couple of things in mind: first, the radiosity matrix is symmetric, so we really only need to compute the upper triangle. Second, the blocks along the diagonal should all be identity blocks as they record the coefficient of the (Subscript[B, i]) in the radiosity equation. Last, we can construct a block matrix using ArrayFlatten as in the example below, which preserves sparsity as well.

A = SparseArray[Partition[Range[4*4],4]];
ArrayFlatten[{{SparseArray[IdentityMatrix[4]],A},{A,SparseArray[IdentityMatrix[4]]}}]

ViewFactorMatrix[SphereList_] :=Module[{SphereTups,i,j},
SphereTups = Table[{SphereList[[i]],SphereList[[j]]},{i,1,Length[SphereList]},{j,1,Length[SphereList]}];
ArrayFlatten[Map[VMBlock[#,Complement[SphereList,#]] &,SphereTups,{2}]]
];

We now need some geometry to work with. In order to get a nice-looking arrangement of spheres, we need to define an energy function which encourages them to cluster, but not to get too close. We try the function

DistScore[x_] := 1/(0.001*x)  +x^3;
Plot[DistScore[x],{x,0,8}]

Now we minimize the function using Mathematica to create some nice geometry. This is an old trick for constructing interesting looking objects; just decide on some general parameters you like, invent a crude metric for scoring various constructions, and then minimize the score numerically. In this case, we want a fairly dense, but not actually touching, collection of unit spheres:

CenterScore[CentersList_] := Module[{i,j,Centers},
Centers = Partition[CentersList,3];
N[Sum[DistScore[Norm[Centers[[i]] - Centers[[j]]]],{i,1,Length[Centers]},{j,i+1,Length[Centers-1]}]]];

InitialCenters = Flatten[Table[RandomReal[{0,8},3],{16}]];
CenterScore[InitialCenters]

MinFinder[InitCenters_] := Module[{nv,mStr,i,Soln},
nv= Length[InitCenters];
mStr = "FindMinimum[CenterScore[{"<>StringJoin[Riffle[Table[ToString[StringForm["t``",i]],{i,1,nv}],","]] <> "}],{" <> StringJoin[Riffle[Table[ToString[StringForm["{t``,``}",i,InitCenters[[i]]]],{i,1,nv}],","]] <> "},AccuracyGoal->2]";
Soln = ToExpression[mStr];
{Soln[[1]],ToExpression["{" <> StringJoin[Riffle[Table[ToString[StringForm["t``",i]],{i,1,nv}],","]] <> "}"] /. Soln[[2]]}

];

Centers = Partition[MinFinder[InitialCenters][[2]],3]

SphereList = Map[BuildSphere[#,1] &,Centers];
Graphics3D[SphereList[[All,7]]]

The MinFinder code uses another good trick; some Mathematica functions require you to spell out a list of variables in advance. But you might not know exactly how many variables are involved! (In this case, I wanted the same MinFinder code to work regardless of the number of spheres.) Luckily, you can use Mathematica code to write Mathematica code; you just construct the input you want as a string and then use ToExpression to evaluate the strings. 

So that we can repeat the experiment, we save the PrettyCenters:

PrettyCenters ={{2.2886596070268563`,2.2859998280192886`,5.486666749520893`},{5.044385830436865`,2.512269847058062`,6.018705678310594`},{3.2083771129750436`,3.209289532159709`,1.3162766956480554`},{1.3417157692440724`,2.9389371956748085`,3.1682909672020747`},{3.1206579350669994`,4.425697006063052`,6.9685633897201775`},{6.138421996057725`,6.374031046956715`,3.5618204679617693`},{1.657560697500216`,5.63095803853559`,2.2840062010783035`},{5.907346224648009`,3.36637183427523`,2.068280899489525`},{3.5190735815373686`,7.376142315786063`,3.249307985935117`},{3.171446676261871`,6.896692218665062`,5.69751297554516`},{5.597553196113725`,5.626619665916129`,6.1910476348409835`},{1.0690220944927311`,5.190237336888103`,4.991314319318154`},{4.253285562115383`,1.5060765767428508`,3.4609002276945504`},{6.76686983998841`,3.7348359629235914`,4.45336879175596`},{4.526627013863203`,5.673526650614165`,1.3241291958391936`},{3.844682904018713`,4.4528463560431355`,4.026129204451894`}};

PrettySpheres = Table[BuildSphere[PrettyCenters[[i]],1.0],{i,1,16}];

Finally, we can generate a nice radiosity matrix:

A = ViewFactorMatrix[PrettySpheres]

MatrixPlot[A]

N[35114 / (512*512)]

As we can see, it’s about 13% full, so fairly sparse, but not amazingly so.

Section 2: The Exact Solution to the Radiosity Model

This 512 x 512 matrix is big enough that we will be able to watch the various iterative methods converge, but not so large that we can’t get the solution in a classical way to compare how fast the iterative methods are converging. So we start by solving using LinearSolve.

We need to start with a RHS. In this case, I want the faces of sphere 7 to glow.

GlowSphere[n_,PrettySpheres_] := Module[{nSpheres,b},
b = Table[0.0,{32 * Length[PrettySpheres]}];
b[[32*(n-1)+1;;32*n]] = Table[5.0,{32}];
b
];

b = GlowSphere[16,PrettySpheres];

Now we can solve the system:

Rad = LinearSolve[A,b];
Graphics3D[DisplaySolution[Rad,PrettySpheres],Lighting->{{"Ambient",GrayLevel[0.4]}}]

Just for fun, here’s a Manipulate where we can play with which sphere is illuminated.

Manipulate[b = GlowSphere[n,PrettySpheres];
Rad = LinearSolve[A,b];
Graphics3D[DisplaySolution[Rad,PrettySpheres],Lighting->{{"Ambient",GrayLevel[0.4]}}],{n,1,16,1}]

Section 3. Iterative Methods.

Our job now is to apply the iterative methods from class: Jacobi, Gauss-Seidel, and SOR(w) to the radiosity matrix that we’ve derived above. Once we see these converging to the correct solution, we’ll use them to construct a much larger radiosity model (which we probably couldn’t do in reasonable time using Mathematica’s standard linear solver). We start with a general form for iterative methods and work from there:

IterativeSolve[{R_,c_},x0_,n_] := NestList[R.# + c &,x0,n];

Jacobi’s Method

Recall that the Jacobi splitting is given by:

A = M - K = D - ((OverTilde[L]) + (OverTilde[U])), so (Subscript[R, J]) = (M^(-1))K = (D^(-1))((OverTilde[L]) + (OverTilde[U]))

and the iteration for solving Ax = b is defined by

(Subscript[X, m+1]) = (Subscript[R, J]) (Subscript[X, m]) + c

where c = (M^(-1))b. In the case of the radiosity matrix, it’s relatively easy to construct these iterates, since the diagonal part of a radiosity matrix is just the identity matrix!

JacobiMethod[A_,b_] := Module[{D,Lhat,Uhat,R,c},

D = DiagonalMatrix[Diagonal[A]];
Lhat = -LowerTriangularize[A,-1];
Uhat = -UpperTriangularize[A,1];

R = Inverse[D].(Lhat + Uhat); 
c = Inverse[D]. b;

{R,c}

];

{Rj,Cj} = JacobiMethod[A,GlowSphere[7,PrettySpheres]];

JacobiIterates = IterativeSolve[{Rj,Cj},b,50];

We start by looking at a collection of Jacobi iterates applied to our original picture to see how the solution changes as we continue to iterate Jacobi’s method:

Manipulate[
Graphics3D[DisplaySolution[JacobiIterates[[n]],PrettySpheres],Lighting->{{"Ambient",GrayLevel[0.4]}}],{n,1,Length[JacobiIterates],1}]

Wow! We can see that by iteration 12 this looks really convincing. To get a sense of what this means in practice, let’s collect a little bit of timing data.

b = GlowSphere[7,PrettySpheres];
Timing[LinearSolve[A,b]][[1]]

{Rj,Cj} = JacobiMethod[A,GlowSphere[7,PrettySpheres]];
Timing[JacobiIterates = IterativeSolve[{Rj,Cj},b,20]][[1]]

Of course, before we get too excited, we should make a plot of the error in the solution:

JacobiIterates = IterativeSolve[{Rj,Cj},b,50];
TrueSolution = LinearSolve[A,b];

JacobiErrorData = -Log[10,Table[Norm[JacobiIterates[[i]] - TrueSolution],{i,1,Length[JacobiIterates]}]];

ListLinePlot[JacobiErrorData,Filling->Axis,PlotMarkers->Automatic,AxesLabel->{Style["Jacobi Iteration",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct digits in Jacobi solution",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]

The Gauss-Seidel Method.

Recall that the Gauss-Seidel splitting is given by:

A = M - K = (D - (OverTilde[L])) - (OverTilde[U]), so (Subscript[R, J]) = (M^(-1))K = (((D-(OverTilde[L])))^(-1))(OverTilde[U])

and the iteration for solving Ax = b is defined by

(Subscript[X, m+1]) = (Subscript[R, J]) (Subscript[X, m]) + c

where c = (M^(-1))b.

GaussSeidelMethod[A_,b_] := Module[{D,Lhat,Uhat,R,c},

D = DiagonalMatrix[Diagonal[A]];
Lhat = -LowerTriangularize[A,-1];
Uhat = -UpperTriangularize[A,1];

R = Inverse[D - Lhat].(Uhat); 
c = Inverse[D - Lhat]. b;

{R,c}

];

{Rgs,Cgs} = GaussSeidelMethod[A,GlowSphere[7,PrettySpheres]];

GaussSeidelIterates = IterativeSolve[{Rgs,Cgs},b,50];

We start by looking at a collection of Gauss-Seidel iterates applied to our original picture to see how the solution changes as we continue to iterate the Gauss-Seidel method:

Manipulate[
Graphics3D[DisplaySolution[GaussSeidelIterates[[n]],PrettySpheres],Lighting->{{"Ambient",GrayLevel[0.4]}}],{n,1,Length[GaussSeidelIterates],1}]

This one looks even better. By iteration 6, this looks really convincing. Again, let’s collect some timing data.

b = GlowSphere[7,PrettySpheres];
Rawtiming = Timing[LinearSolve[A,b]][[1]]

{Rgs,Cgs} = GaussSeidelMethod[A,GlowSphere[7,PrettySpheres]];
GStiming  = Timing[GaussSeidelIterates = IterativeSolve[{Rgs,Cgs},b,12]][[1]]

{Rj,Cj} = JacobiMethod[A,GlowSphere[7,PrettySpheres]];
Jtiming = Timing[JacobiIterates = IterativeSolve[{Rj,Cj},b,20]][[1]]

We can see that the GS method offers some savings over the Jacobi method, taking only 75% of the time. (This fraction is likely to be implementation-dependent, so be a little careful about taking this kind of thing literally.)

100*(GStiming/Jtiming)

Of course, before we get too excited, we should make a plot of the error in the solution:

GaussSeidelIterates = IterativeSolve[{Rgs,Cgs},b,20];
TrueSolution = LinearSolve[A,b];

GaussSeidelErrorData = -Log[10,Table[Norm[GaussSeidelIterates[[i]] - TrueSolution],{i,1,Length[GaussSeidelIterates]}]];

ListLinePlot[{JacobiErrorData,GaussSeidelErrorData},Filling->Axis,PlotMarkers->Automatic,AxesLabel->{Style["Iteration",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct digits in Jacobi and GS solutions",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica")]

Successive Overrelaxtion Method

Recall that the SOR(ω) splitting is given by:

A = M - K = (1/ω)(D - ω(OverTilde[L])) - (1/ω)((1 - ω) D + ω(OverTilde[U])), so (Subscript[R, J]) = (M^(-1))K = (((D- ω(OverTilde[L])))^(-1))((1-ω)D + ω(OverTilde[U]))

and the iteration for solving Ax = b is defined by

(Subscript[X, m+1]) = (Subscript[R, J]) (Subscript[X, m]) + c

where c = (M^(-1))b. This is straightforward to code:

SORMethod[A_,b_,w_] := Module[{D,Lhat,Uhat,R,c},

D = DiagonalMatrix[Diagonal[A]];
Lhat = -LowerTriangularize[A,-1];
Uhat = -UpperTriangularize[A,1];

R = Inverse[D - w *Lhat].((1 - w)*D + w *Uhat); 
c = w Inverse[D - w *Lhat]. b;

{R,c}

];

We face a new wrinkle here: how should we pick ω? We recall that we proved that in general SOR(ω) should converge for 0 < ω < 2, but it is equal to Gauss-Seidel for ω = 1, and only overrelaxed when ω > 1. So we will try three versions of SOR: SOR(1.25), SOR(1.5) and SOR(1.75).

{Rsor85,Csor85} = SORMethod[A,GlowSphere[7,PrettySpheres],0.85];
{Rsor125,Csor125} = SORMethod[A,GlowSphere[7,PrettySpheres],1.25];
{Rsor15,Csor15} = SORMethod[A,GlowSphere[7,PrettySpheres],1.5];
{Rsor175,Csor175} = SORMethod[A,GlowSphere[7,PrettySpheres],1.75];
{Rsor105,Csor105} = SORMethod[A,GlowSphere[7,PrettySpheres],1.05];
{Rsor1,Csor1} = SORMethod[A,GlowSphere[7,PrettySpheres],1];
{Rsor102,Csor102} = SORMethod[A,GlowSphere[7,PrettySpheres],1.02];

SOR125Iterates = IterativeSolve[{Rsor125,Csor125},b,20];
SOR15Iterates = IterativeSolve[{Rsor15,Csor15},b,20];
SOR175Iterates = IterativeSolve[{Rsor175,Csor175},b,20];
SOR102Iterates = IterativeSolve[{Rsor102,Csor102},b,20];

We start by looking at a collection of SOR iterates applied to our original picture to see how the solution changes as we iterate the SOR method at various values of ω:

Manipulate[
Iterates = Algorithm /. {"Jacobi"->JacobiIterates,"GS" -> GaussSeidelIterates, "SOR(1.02)" -> SOR102Iterates,"SOR(1.25)" -> SOR125Iterates,"SOR(1.5)" -> SOR15Iterates, "SOR(1.75)" -> SOR175Iterates};
Graphics3D[DisplaySolution[Iterates[[n]],PrettySpheres],Lighting->{{"Ambient",GrayLevel[0.4]}}],{{n,1,"Iteration"},1,Length[GaussSeidelIterates],1},{Algorithm,{"Jacobi","GS","SOR(1.02)","SOR(1.25)","SOR(1.5)","SOR(1.75)"}},ContinuousAction->False]

This one certainly looks different for the various ω values, but not neccesarily better for the larger ones. So let’s make a plot of the error in the solution and try to see how it depends on ω. This needs to be a manipulate, because we are interested in seeing how the plot changes as we alter ω:

Manipulate[
{Rsor,Csor} = SORMethod[A,GlowSphere[7,PrettySpheres],w];
SORiterates = IterativeSolve[{Rsor,Csor},b,20];
SORErrorData = -Log[10,Table[Norm[SORiterates[[i]] - TrueSolution],{i,1,Length[SORiterates]}]];
ListLinePlot[{JacobiErrorData,GaussSeidelErrorData,SORErrorData},Filling->Axis,PlotMarkers->Automatic,AxesLabel->{Style["Iteration",Larger,FontFamily->"Helvetica",Bold],Style["Digits",Larger,FontFamily->"Helvetica",Bold]},PlotLabel->Style["Correct digits in Jacobi, GS, and SOR(ω) solutions",Larger,FontFamily->"Helvetica",Bold],AxesStyle->(FontFamily->"Helvetica"),ImageSize->600],{{w,1,"ω"},0,2},ContinuousAction->False]

This result seems really fascinating to me! As far as I can tell, going to SOR is quite flaky on this problem, and indeed only faster than GS for very small and very carefully selected ω values. I don’t think that this is a bug (it’s kind of hard to see a bug in the SOR code, since it’s so simple, so I’m really interested to see the explanation for this phenomenon).

3. The Big Example

Having seen how fast these methods run on this toy problem, it’s interesting to try them out on a larger-scale radiosity calculation. For this, of course, we’ll need a larger collection of spheres.

InitialCenters = Flatten[Table[RandomReal[{0,8},3],{64}]];
CenterScore[InitialCenters]
Centers = Partition[MinFinder[InitialCenters][[2]],3]
SphereList = Map[BuildSphere[#,0.65] &,Centers];
Graphics3D[SphereList[[All,7]]]

BigPrettyCenters = {{2.2740950895299554`,4.740075868263721`,4.495052839535714`},{5.093056599210722`,6.900477637778688`,2.1434252251072583`},{6.2351392606847265`,1.736315199643782`,3.641470091146464`},{7.1660354113960025`,3.4001851467260154`,4.48599979494305`},{3.9978193949384138`,4.1308651943776145`,3.7785077882682354`},{3.131690379998294`,3.050478161402028`,0.7497003305739139`},{4.765504378038178`,5.643697911629409`,2.9405844989794097`},{2.448847882582285`,2.352512245695755`,6.139891350534459`},{3.2610838331234833`,5.107678717298819`,0.696350277791504`},{3.7206396716307335`,5.916587272303049`,4.456676768992831`},{3.2789879677578306`,1.64157732167617`,1.7882584850070657`},{4.635616353459288`,2.375514146396928`,1.146938716965054`},{6.3959310443165664`,3.7475739528144065`,5.887209794392837`},{2.1379709205325486`,1.5750326263236067`,3.0427837013158907`},{4.719475673668796`,3.96110609726016`,0.5367623193952286`},{6.4643318581785865`,6.252735909356443`,3.5499244123585467`},{1.9471672001429452`,4.08153131953299`,6.35106887347859`},{5.403468942610716`,4.949714900557432`,4.631701001532983`},{6.040430536809812`,2.0681854475990513`,5.381285369742642`},{5.184690950888709`,7.203327063075141`,3.854312312984883`},{4.422720522496958`,6.185429473933341`,6.302612041758114`},{5.706626939228236`,6.5635189370180544`,5.217247384995836`},{4.96730399153369`,1.2060987134412167`,2.586952865109839`},{2.179463191189814`,6.8422923003673874`,4.496515955918441`},{5.677853414654981`,3.199243362886643`,4.080669783689459`},{1.1400526330071101`,5.829518045408407`,3.54451114570676`},{4.11041208835816`,2.4604399088344167`,3.0881668638571074`},{1.1212278780027585`,3.296229133003629`,5.267388426449819`},{2.7161298171412622`,6.294578385619316`,5.949222074869983`},{4.7420879621759005`,4.571881158229435`,6.992488941856955`},{3.678472801409677`,7.370729238198435`,3.248166564895779`},{2.915201912859241`,5.338861428805338`,2.786929163064681`},{4.878072379834104`,5.547442831452923`,1.00658145331402`},{7.089588235050072`,4.756722290361729`,3.1379861533794733`},{6.214629797206607`,4.005333461810647`,1.2706705696482146`},{6.99043836903991`,5.142266389240784`,4.778454018926507`},{4.415714755745422`,2.799259788065253`,5.100976867085535`},{5.193723981689503`,3.0403121950979237`,6.636617595323946`},{4.212854983476078`,1.7148798011729425`,6.143313790286217`},{5.93902491869871`,5.264347921450217`,6.210566394440857`},{0.9345089340739213`,4.397759585098591`,2.5117690205797407`},{0.659758374962512`,4.405262287232518`,4.197056422537024`},{3.4863028646599052`,3.2826695385238693`,6.933670651967107`},{3.9140299804921077`,7.306510000876147`,5.056622420856835`},{1.661620326404187`,1.802330395021509`,4.561530167717098`},{5.464719329045474`,4.054098433626446`,2.578362430240529`},{3.509511114063874`,6.541336055558709`,1.5828563757852534`},{1.7371733531629296`,2.605721195564736`,1.8002847051410857`},{3.5400350640705747`,0.8885289580738025`,3.422116005720985`},{5.951336680823175`,2.3700250968759335`,1.939445944907746`},{6.293265882050611`,5.596994012517881`,1.932197657787541`},{3.062599012048457`,5.071502178997133`,6.902485332218903`},{3.18859713265176`,1.1902157955214108`,4.944752693441965`},{3.821662789674708`,3.993597474997181`,1.8783198086663604`},{2.3922111179670034`,3.497959232161391`,2.9361107747398996`},{0.90836588129929`,2.887718829692632`,3.4240143142860253`},{3.811035563989462`,4.4413327573896675`,5.616196958764813`},{4.8909970034322985`,1.1043197464066854`,4.48500432593918`},{1.9595651745934095`,4.153493489156728`,1.1846120472715664`},{2.754994874898982`,2.9527622124241653`,4.604921127690306`},{7.003000214467353`,3.1795040082087813`,2.8133061903710743`},{1.3403502574381987`,5.466686161556347`,5.482892632962236`},{1.8125145196478838`,5.703135839603253`,1.7140145416671526`},{2.262724478313045`,6.872616582442896`,2.8470195339289135`}};

BigPrettySpheres = Map[BuildSphere[#,0.65] &,BigPrettyCenters];

"Timing[BigA = ViewFactorMatrix[BigPrettySpheres]][[1]]"

It turns out to take almost an hour for Mathematica to compute the radiosity matrix itself. 

Timing[BigA = ViewFactorMatrix[BigPrettySpheres]][[1]]
2763.448156000001`

so we’re commenting out the line. This is basically an artifact; we’re using the world’s least efficient hidden surface removal code, so we shouldn’t pay too much attention to matrix construction at this point. Still, there’s a life lesson to learn here, which is that often constructing the matrices involved in a linear algebra problem is significantly more expensive than solving the resulting problem. So it can make a lot of sense to work hard on optimizing that portion of your calculation (whatever it is) rather than implementing or buying a really exotic linear algebra solver.

Dimensions[BigA]

In order to avoid repeating this process, we exported the resulting matrix to the Harwell-Boeing sparse matrix format:

Export["/Users/cantarel/courses/math4500/BigRadiosityMatrix.rsa",BigA]
"/Users/cantarel/courses/math4500/BigRadiosityMatrix.rsa"

"Export["/Users/cantarel/courses/math4500/BigRadiosityMatrix.rsa",BigA]"

And now we’ll import the data back from disk in order to save time for this classroom demo. Note that you’ll probably have to worry about filenames if you try this at home:

BigADisk = Import["/Users/cantarel/courses/math4500/BigRadiosityMatrix.rsa"]

We are finally ready to try out our iterative solver! We’ll give GS a crack at it:

Bigb =Sum[RandomReal[{0.25,1.25}]*GlowSphere[RandomInteger[{1,64}],BigPrettySpheres],{5}];
{BigRgs,BigCgs} = GaussSeidelMethod[BigADisk,Bigb];

Timing[IterativeSolve[{BigRgs,BigCgs},Bigb,70]][[1]]

BigGaussSeidelIterates = IterativeSolve[{BigRgs,BigCgs},Bigb,70];

We start by looking at a collection of Gauss-Seidel iterates applied to our original picture to see how the solution changes as we continue to iterate the Gauss-Seidel method:

Manipulate[
Graphics3D[DisplaySolution[BigGaussSeidelIterates[[n]],BigPrettySpheres],Lighting->{{"Ambient",GrayLevel[0.2]}},Boxed->False],{{n,20},1,Length[BigGaussSeidelIterates],1}]