One of my students pointed out this video showing how it’s done, from Ireland.

I was forced to broaden my horizons in my PhD pursuit.  That’s because, although I’ve done coding since I was eighteen, I had to acquire a deeper understanding for two things: linear algebra and numerical methods.  It’s no understatement to say that both of these are at the core of the advances wrought by computerisation, whether we’re talking about statistical analysis or (in my case) simulation.

After my initial boffo performance, I turned to my Iranian friends for more help.  So they let me use some of the books they found useful for study back in the “old country”.  One of those was a sizeable book entitled Applied Numerical Methods by Brice Carnahan, H.A. Luther and James O. Wilkes.  As was the case with their wedding video, the heart skipped a beat, because the middle author, Hubert A. Luther, was my Differential Equations teacher at Texas A&M, forty years ago this spring.

Applied Numerical Methods was, AFAIK, the first really comprehensive textbook which combined linear algebra, numerical methods, and coding (in their case, FORTRAN IV) in one text.  Although some of the methodologies have been improved since it was published in 1969, and languages have certainly changed, it’s still a very useful book, although a little dense in spots.  Many of the books on the subject that have come afterwards have learned from its mistakes, but still refer back to the original.

Dr. Luther taught me the last required math class in my pursuit of an engineering degree at Texas A&M.  It wasn’t an easy class, even after three semesters of calculus (which I did reasonably well at).  Although he was originally from Pennsylvania, he acclimated himself to the Lone Star State with western shirt, belt and string tie, the only professor I can remember who did so. The start to his course was especially rough; the textbook was terrible, he was a picky grader, the scores I got back were low.  I thought I was facing the abyss…until another one of those “aha” moments came along.

We (the engineering students) were standing outside our Modern Physics class, which came before Differential Equations.  I found out I wasn’t the only one having this problem.  But one of my colleagues, a Nuclear Engineering student who went on to become my class’ wealthiest member, had a simple suggestion.  Go visit his office, he said.  He’s lonely (he was nearing retirement) and likes the company.  Your grade will go up.

I wasn’t much for visiting my professors, but I was desperate enough to try anything.  I made a couple of office visits.  I’m not sure how helpful his advice was, but his grading became more lenient and I got through the course OK.

Today I’m on the other end of the visitation.  I spend a lot of time in the office with no student visits.  Part of the problem comes from scheduling, both theirs and mine.  But I’ve found out something else about student visits: the students that come to see you really care about what they’re supposed to be doing in your class.  Although there are still students who think it their duty to “tough it out” without asking questions, many others just want to get through in the quickest and least time-consuming way they can find.

I’m glad I took my classmate’s advice and made the office visits.  But there are two other lessons I have learned since that time.

The first is that I wish I had taken a numerical methods course taught by Dr. Luther, it would have prepared me for what I’ve been doing both before and during the time of my PhD pursuit.

The second is that, when I started my MS degree twenty years later, I took a course over basically the same material taught by a Russian.  I found out that there was a great deal I hadn’t learned from Dr. Luther, and that American math education leaves a lot to be desired of.  So sometimes making the way easier up front comes back to get you in the end.

Although it’s been out there for a long time, these days the clash between religion and science has been especially heated.  The simplest way to solve the problem would be to cancel elections and let the self-proclaimed know-it-alls run the show.  In this way they could ignore the religious “masses” and insure the continuity of their funding.  Funded scientists are happy scientists…

But what would happen if there was a synergy between the two?  Basically the same thing that is happening between religion and science now: an academic slugfest.  (Reminds one of Rabbi Jonathan Sacks’ joke about “the tradition”…)    That’s pretty much the bottom line of the life of the life of Georg Cantor (1845-1918) and his formulation of set theory.

He was born in St. Petersburg, Russia to parents who originally came from Denmark.  When he was eleven, they moved to Frankfurt, in the German Electorate of Hesse (the Hessians were the ones George Washington crossed the Delaware to defeat at Trenton).  As Carl Boyer notes in his A History of Mathematics:

His (Cantor’s) parents were Christians of Jewish background–his father had been converted to Protestantism, his mother had been born a Catholic.  The son Georg took a strong interest in the finespun arguments of medieval theologians concerning continuity and the infinite, and this militated against his pursuing a mundane career in engineering as suggested by his father.  In his studies at Zurich, Göttingen and Berlin the young man consequently concentrated on philosophy, physics and mathematics–a program that seems to have fostered his unprecedented mathematical imagination.

His central “claim to fame” is the elucidation of set theory (or, as the Germans are wont to call it, mengenlehre).  It’s not an understatement that set theory has come to dominate the teaching of mathematics and its conceptualisation, as I found out the hard way taking advanced linear algebra a couple of years ago, complete with the bizarre notation that has just about taken over math textbooks.  In the 1960’s it was the centrepiece of the “new math” that came into primary and secondary school curricula, and that was controversial, but a great deal more useful than its critics would admit.

The controversy didn’t end with the sets themselves.  Cantor realised that set theory forced him to consider something that mathematicians had danced around for almost two centuries: infinite quantities, or more precisely transfinite quantities.  Sets can have an infinite number of elements, but just what that means was something Cantor plunged very deeply into.

It’s easy to get lost in Cantor’s reasoning, as the concepts he proposed are very profound.  I’ll try to keep things as uncomplicated as I can, taking the risk that I may oversimplify the business.

Let us consider the set of integers.  We know instinctively that there are an infinite number of integers.  Now let us consider the number of even integers.  You’d think that there are half as many even integers as all integers, right?  But both quantities are in fact infinite, which means that dividing it by two doesn’t mean much.  In fact Cantor proved that, if we considered the set of all integers and the set of even integers,  we would have a one-to-one correspondence between each member of each set.  So the size of the two sets is equal, even though one set is a subset of the other.

Things get more complicated when we pass from integers and rational numbers to transcendental numbers like e and pi. Cantor proved that the number of transcendental numbers was larger than that of either/or or both/and the integers and rational numbers, even though all of them were infinite.  Cantor had shown, in effect, that not all infinities were equal to each other!

One device that Cantor, and just about anyone else who deals with transfinite numbers, uses is the limit.  But one major difference between Cantor and many of his contemporaries–and predecessors–is that Cantor showed that infinity was in fact an existing quantity, the problem with the transcendentals not withstanding.

That lit several fuses.  Before Cantor’s time the French mathematician Cauchy stated the following:

I protest against the use of an infinite magnitude as something completed, which is never permissible in mathematics.  Infinity is merely a way of speaking, the true meaning being a limit which certain ratios approach indefinitely close.

The most deadly grenade pin that Cantor pulled, however, was that of Leopold Kronecker (1823-1891), after whom the famous Kronecker Delta is named.  Kronecker, like Cantor of Jewish origins but a Christian, famously stated that “God made the integers, and all the rest is the work of man.” Kronecker made a career out of academically trashing Cantor, blocking appointments and delaying publications.  Cantor, not the scholarly pugilist the situation called for (he should have read Jerome with the mediaevals) had his first nervous breakdown in 1884.  After that time he published little, and died in a psychiatric hospital in 1918, although by then his work was receiving the recognition it deserved.   David Hilbert pithily stated that “From the paradise created for us by Cantor, no one will drive us out”.

So how did the mediaevals influence this revolution in mathematics? The problem of infinity wasn’t as far-fetched as you might think.  It had sat there since Newton and Leibniz set forth the calculus, which in turn hangs on infinitesimals.  Between two finite points there is an infinite number of infinitesimals.  Mediaevals have been jeered for wondering how many angels could dance on the head of a pin, but as long as they were infinitesimals, the answer is clear for finite pins.  It was only a matter of someone putting infinitesimals and infinities together, and that person (with help from others) was Cantor.

Anyone who has explored the philosophy of the scholastics with a mathematical background sooner or later will consider the relationship between their idea and the mathematics of infinity.  Coming off of a master’s degree, I found myself doing that in My Lord and My God.  Although I would not dare to rank myself any where near Cantor, I discovered that all infinities were not equal, and, although they could not have a finite ratio with finite quantities, they were not necessarily equal to each other.  That in turn helped me to see that subordination in God does not impair the deity of the subordinate persons, which solves many problems.  Unfortunately there are those who either can’t or–ahem–won’t see that relationship, and there is always the problem that prelates and seminary academics and are often mathematically challenged.

Today we live in a world where science and religion are forcibly bifurcated.  But it was not always so.  Cantor–and Kronecker and others for that matter–allowed the two to intermingle, and before that Euler was more religiously conservative than Voltaire.  And the Nineteenth Century in Europe was a golden age in mathematics, where advances came one after the other.

But there’s a price.  If you want to get into serious trouble, read the mediaevals, and that’s true for mathematicians and theologians alike.

Note: in addition to Boyer’s book, I used Jane Muir’s Of Men and Numbers: The Story of the Great Mathematicians (Dover Books on Mathematics) in writing this piece.

The purpose of this piece is to document the numerical integration of a function in a two-dimensional space using the conjugate gradient method.

Basic Problem Statement and Closed Form Solution

The boundary value problem is as follows:

where

A plot of the right hand side of this is shown below.

To solve this in closed form, we need to consider our solution method. Given the form of Poisson’s Equation above, it is tempting to use a double Fourier series (or alternatively a complex exponential solution with constant coefficients) with only a few (1-2) terms in the solution. Complicating this is the presence of the first order term for . One way to deal with this is to consider including the linear term in an assumed solution, as follows:

Transforming this into sines and cosines with a more general solution, we can write

Applying the operator to this yields

The idea now is to eliminate terms of sine and cosine combinations that do not appear in the second partial derivatives.

The simplest place to start is to state the following, based on the structure of same partial derivatives:

This harmonizes the sine and cosine functions. By inspection, we can also state that

since these terms do not appear in the problem statement. The values of and are not as obvious, but based on the previous statement our assumed solution (equated to the problem statement) reduces to

This can be reduced to a linear system as follows:

In matrix form, the equation becomes

Inverting the matrix and multiplying it by the right hand side yields

This yields our solution,

A plot is shown below.

Solution by Conjugate Gradient

The residual norm history is shown below for n=9, 19 and 39, where n is the number of nodes on each axis to produce a grid. The plot shows the decrease in the 2-norm of the residual as the iterations progress. The iterations are stopped when they either exceed 2000 or if the 2-norm falls below .

All of the discretisations converged within 40 steps. In any event the conjugate gradient method will converge in no more steps than the row/column number of the matrix, and in this case the performance of the method far exceeded that, doubtless in part to the fact that we employed a preconditioner to accelerate the solution.

As mentioned in the code comments, the method used is conjugate gradient as outlined in Gourdin and Boumahrat (2003). The method begins by computing the original residual:

We then solve the following equation

At this point we need to consider the conditioning matrix . It is defined as follows:

where is the lower diagonal portion of . As an aside, this is an incomplete Cholesky factorization. Obviously is the upper diagonal, being symmetric. Since we have an upper diagonal solver, it would make sense to compute and invert using .

To accomplish this, we invert this to yield

We can accomplish the first term of the right hand side with the information at hand. For the second term, considering

we transpose both sides to yield

and substituting

We can use the factoring method to invert , then take its transpose, then multiply by its transpose to yield , which is in essence our conditioning matrix and which can be used to solve for the vector .

Returning to the algorithm, to save calculations later we define the
vector

and solving

we equate

We now begin our iteration for . We first compute

for which purpose we developed a special subroutine. We then update our step as follows:

Using the conditioning matrix we developed, we compute the following

and (using a similar routine we used for )

and

From here we index and repeat the cycle until either we reach or our convergence criterion.  It should be noted that the efficiency of conjugate gradient was such that, particularly at the finer discretizations, the routine spent more time generating the preconditioning matrices than it did in actually going through the cycles!

Solution Code

The solution codes were written in FORTRAN 77, and is presented below.  The traditional FORTRAN column spacing got messed up in the cut and paste. Also presented is some of the code which was incorporated using the INCLUDE statement of Open WATCOM FORTRAN. The size of the grid was varied using the parameter statement at the beginning of the problem,and is generally set to 39 for the codes presented. It was varied to 9 and 19 as needed.

We should also note that the routine was run in single precision.

 c Solution of Two-Dimensional Grid Using
 c Conjugate Gradient Iteration
 c Implementation of conjugate gradient method
 c as per Gourdin and Boumahrat (2003)
 c Includes preconditioning by construction of matrix C
 c Matrix size changed via parameter statement
 include 'mpar.for'
 parameter(pi=3.141592654,istep=2000)
 c Define arrays for main array, diagonal block, off-diagonal block,
 c and solution and rhs vectors
 dimension a(nn,nn),d(n,n),c(n,n),u(nn),b(nn)
 c Define conditioning and related matrices
 dimension cond(nn,nn),t(nn,nn),tt(nn,nn)
 c Define residual vectors and norm for residual
 dimension r(nn),rold(nn),p(nn),s(nn),sold(nn),xnorm(istep)
 c Coordinate Arrays
 dimension x(n,n),y(n,n)
 c Error Vector for Final Iterate
 dimension uerr(nn)
 c Statement function to define rhs
 f(x1,y1)=-5*y1*sin(x1)*sin(2*y1)+4*sin(x1)*cos(2*y1)
 c Statement function to define closed form solution
 fexact(x1,y1)=y1*sin(x1)*sin(2*y1)
 c Write Header for Output
 write(*,*)'Conjugate Gradient Iteration to solve',
 &' two-dimensional grid'
 write(*,*)'Grid Size = ',n,' x ',n
 write(*,*)'Matrix Size ',nn,' x ',nn
 call tstamp
 c Initialise diagonal block array
 do 20 i=1,n,1
 do 20 j=1,n,1
 if(i.eq.j) then
 d(i,j)=4.
 elseif(abs(i-j).eq.1) then
 d(i,j)=-1.
 else
 d(i,j)=0.
 endif
 20 continue
 c Initialise off-diagonal block array
 do 30 i=1,n,1
 do 30 j=1,n,1
 if(i.eq.j) then
 c(i,j)=-1.
 else
 c(i,j)=0.
 endif
 30 continue
 c Compute grid spacing h
 xh=pi/float(n+1)
 c Initialise rhs
 do 60 j=1,n,1
 do 60 i=1,n,1
 x(i,j)=float(i)*xh
 y(i,j)=float(j)*xh
 ii=(j-1)*n+i
 b(ii)=-f(x(i,j),y(i,j))*xh**2
 write(*,61)i,j,ii,x(i,j),y(i,j),b(ii)
 61 format(3i5,3f10.3)
 60 continue
 c Insert block matrices into main matrix
 call blkadd(d,c,a)
 c Initialise result vector
 do 70 ii=1,nn,1
 70 u(ii)=1.0
 c Compute first residual vector, using residual vector as temporary storage
 call xmvmat(a,u,r)
 do 75 ii=1,nn,1
 r(ii)=b(ii)-r(ii)
 75 continue
 c Construct upper triangular matrix tt by stripping lower diagonal portion of
 c matrix a
 do 76 ii=1,nn,1
 do 76 jj=1,nn,1
 if(jj-ii)78,77,77
 77 tt(ii,jj)=a(ii,jj)
 goto 76
 78 tt(ii,jj)=0.0
 76 continue
 c Invert upper triangular matrix tt by factorisation method
 call uminv(tt)
 c Construct inverted matrix t by transposing inverted matrix tt
 call trnspz(tt,t)
 c Multiply t by tt to obtain preconditioning matrix c (which is actually
 c c**(-1))
 call xmvmul(tt,t,cond)
 c Compute initial vector p by multiplying cond (inverse of c) by residual
 call xmvmat(cond,r,p)
 c Set initial vector s to p
 do 79 ii=1,nn,1
 79 s(ii)=p(ii)
 do 80 kk=1,istep,1
 c set old values of vectors r and s
 do 190 ii=1,nn,1
 rold(ii)=r(ii)
 190 sold(ii)=s(ii)
 c Compute alpha for each step
 alpha1=alpha(a,p,r,s)
 c Update value of u
 do 181 ii=1,nn,1
 181 u(ii)=u(ii)+alpha1*p(ii)
 c Update residual for each step r = b-Au
 call xmvmat(a,u,r)
 do 82 ii=1,nn,1
 r(ii)=b(ii)-r(ii)
 82 continue
 c Update vector s
 call xmvmat(cond,r,s)
 c Compute value of beta for each step
 beta1=beta(r,rold,s,sold)
 c Update value of p vector
 do 195 ii=1,nn,1
 195 p(ii)=s(ii)+beta1*p(ii)
 c Invoke function to compute Euclidean norm of residual
 c Kicks the iteration out once tolerance is reached
 xnorm(kk)=vnorm(r,nn)
 if(xnorm(kk).lt.1.0e-04)goto 83
 write(*,*)kk,xnorm(kk)
 80 continue
 83 continue
 do 90 j=1,n,1
 do 90 i=1,n,1
 ii=(j-1)*n+i
 uerr(ii)=fexact(x(i,j),y(i,j))-u(ii)
 90 continue
 uerrf=vnorm(uerr,nn)
 write(*,*)'Number of iterations = ',kk
 write(*,*)'Euclidean norm for final error =',uerrf
 open(2,file='m5610p1d.csv')
 if(kk.gt.istep)kk=istep
 do 40 ii=1,kk,1
 write(2,50)ii,xnorm(ii)
 50 format(i5,1h,,e15.5)
 40 continue
 close(2)
 stop
 end

c Subroutine to insert nxn blocks into n**2xn**2 array
 c Assumes all diagonal blocks are the same
 c Assumes all off-diagonal blocks (upper and lower) are the same
 c Assigns zero values elsewhere in the array
 subroutine blkadd(d,c,a)
 include 'm5610par.for'
 dimension a(nn,nn),d(n,n),c(n,n)
 c Insert blocks into array
 do 20 ii=1,n,1
 do 20 jj=1,n,1
 c Determine row and column index in master array for corner of block
 ic=(ii-1)*n+1
 jc=(jj-1)*n+1
 do 30 i=1,n,1
 do 30 j=1,n,1
 iii=ic+i-1
 jjj=jc+j-1
 if(ii.eq.jj)then
 a(iii,jjj)=d(i,j)
 elseif(abs(ii-jj).eq.1)then
 a(iii,jjj)=c(i,j)
 else
 a(iii,jjj)=0.
 endif
 30 continue
 20 continue
 return
 end

c Subroutine to multiply a nn x nn matrix with an nn vector
 c a is the nn x nn square matrix
 c b is the nn vector
 c c is the result
 c nn is the matrix and vector size
 c Result is obviously an nn vector
 c Routine loosely based on Gennaro (1965)
 subroutine xmvmat(a,b,c)
 include 'm5610par.for'
 dimension a(nn,nn),b(nn),c(nn)
 do 20 i=1,nn,1
 c(i)=0.0
 do 20 l=1,nn,1
 20 c(i)=c(i)+a(i,l)*b(l)
 return
 end

c Subroutine to multiply a nn x nn matrix with another nn x nn matrix
 c a is the first nn x nn square matrix
 c b is the second nn x nn matrix
 c c is the result
 c nn is the matrix size
 c Result is obviously an nn x nn matrix
 c Routine based on Gennaro (1965)
 subroutine xmvmul(a,b,c)
 include 'm5610par.for'
 dimension a(nn,nn),b(nn,nn),c(nn,nn)
 do 20 i=1,nn,1
 write(*,*)'Multiplying Row ',i
 do 20 j=1,nn,1
 c(i,j)=0.0
 do 20 l=1,nn,1
 20 c(i,j)=c(i,j)+a(i,l)*b(l,j)
 return
 end

c Subroutine to invert a mm x mm upper triangular matrix using factor method
 c u is input and output matrix
 c t is result matrix, written back into the output matrix
 c mm is matrix size
 c Result overwrites original matrix
 c Based on Gennaro (1965)
 subroutine uminv(u)
 include 'm5610par.for'
 dimension u(nn,nn),t(nn,nn)
 c Zero all entries in matrix t except for diagonals, which are the reciprocals
 c of the diagonals of the orignal matrix
 do 20 i=1,nn,1
 do 20 j=1,nn,1
 if(i-j)11,10,11
 10 t(i,j)=1.0/u(i,j)
 goto 20
 11 t(i,j)=0.0
 20 continue
 c Compute strictly upper triangular entries of inverted matrix
 do 40 i=1,nn,1
 write(*,*)'Inverting Row ',i
 do 40 j=1,nn,1
 if(j-i)40,40,31
 31 l=j-1
 do 41 k=i,l,1
 t(i,j)=t(i,j)-t(i,k)*u(k,j)/u(j,j)
 41 continue
 40 continue
 c Write back result into original input matrix
 do 50 i=1,nn,1
 do 50 j=1,nn,1
 50 u(i,j)=t(i,j)
 return
 end

c Function to compute Euclidean norms of vector
 function vnorm(V,no)
 dimension V(no)
 znorm=0.
 do 100 ia=1,no,1
 100 znorm=znorm+V(ia)**2
 vnorm=sqrt(znorm)
 return
 end

c Function to determine scalar multiplier alpha for conjugate gradient method
 c Function returns a scalar which is multiplied by residual vector
 c Function uses subroutine xmvmat to multiply matrix A by residual
 c A is main matrix
 c r is residual vector
 c nn is size of matrix/vector
 c rgrad is scalar multiplier for gradient method
 function alpha(a,p,r,s)
 include 'm5610par.for'
 dimension a(nn,nn),r(nn),pt(nn),p(nn),s(nn)
 c Compute dot product of r and s for numerator
 rnumer=0.0
 do 10 ii=1,nn,1
 rnumer=rnumer+r(ii)*s(ii)
 10 continue
 c Compute matrix product of A * p
 call xmvmat(a,p,pt)
 c Premultiply matrix product by transpose of p vector for denominator
 rdenom=0.0
 do 20 ii=1,nn,1
 rdenom=rdenom+p(ii)*pt(ii)
 20 continue
 alpha = rnumer/rdenom
 return
 end

c Function to determine scalar multiplier beta for conjugate gradient method
 c Function returns a scalar which is multiplied by p vector
 c r,s are vectors is residual vector
 c nn is size of matrix/vector
 c rgrad is scalar multiplier for gradient method
 function beta(r,rold,s,sold)
 include 'm5610par.for'
 dimension r(nn),rold(nn),s(nn),sold(nn)
 c Compute dot product for numerator
 rnumer=0.0
 do 10 ii=1,nn,1
 rnumer=rnumer+r(ii)*s(ii)
 10 continue
 c Compute dot product for denominator
 rdenom=0.0
 do 20 ii=1,nn,1
 rdenom=rdenom+rold(ii)*sold(ii)
 20 continue
 beta = rnumer/rdenom
 return
 end

c Subroutine to transpose a matrix
 c Adapted from Carnahan, Luther and Wilkes (1969)
 subroutine trnspz(a,at)
 include 'm5610par.for'
 dimension a(nn,nn),at(nn,nn)
 do 14 ii=1,nn,1
 do 14 jj=1,nn,1
 14 at(jj,ii)=a(ii,jj)
 return
 end

include 'tstamp.for'

Included Routines

• tstamp.for calls a OPEN WATCOM function and time stamps the output.
• mpar.for sets the parameter statements and is shown below.

 parameter(n=39,nn=n*n)

References

• Carnahan, B., Luther, H.A., and Wilkes, J.O. (1969) Applied Numerical Methods. New York: Wiley.
• Gennaro, J.J. (1965) Computer Methods in Solid Mechanics. New York: Macmillan.
• Gourdin, A., and Boumahrat, M. (2003) Applied Numerical Methods. New Delhi: Prentice-Hall India.