The following solution was supplied by Brian Medley
SIAM Challenge Problem 7:
Since we only want the 1,1 element of the inverse, it is enough to find the
vector x (x_1, x_2, ... , x_20000) for which Ax = (1,0,0,...,0) and x_1 is
then the number we're seaching for. In other words, we are solving a system
of linear equations
(a_1,1)*(x_1) + ... + (a_1,20000)*(x_20000) = 1
(a_2,1)*(x_1) + ... + (a_2,20000)*(x_20000) = 0
...
(a_20000,1)*(x_1) + ... + (a_20000,20000)*(x_20000) = 0
20000 seems chosen to rule out exact methods of solving, so I used Gauss
Seidel iteration:
(x_1)_(i+1) = {1-(a_1,2)*(x_2)_i - ... - (a_1,20000)*(x_20000)_i}/(a_1,1)
(x_2)_(i+1) = {-(a_2,1)*(x_1)_(i+1) - (a_2,3)*(x_3)_i... -
(a_1,20000)*(x_20000)_i}/(a_2,2)
...
(x_20000)_(i+1) = {-(a_20000,1)*(x_1)_(i+1) - (a_20000,2)*(x_2)_i... -
(a_20000,19999)*(x_19999)_i}/(a_20000,20000)
Since we are dividing by the diagonal coefficients, which are large, the
iteration looked very likely to converge. Rather than stopping to prove
this rigorously, I just went ahead and tried it with the UBASIC program
below. (No special reason for using UBASIC other than it was what I had
been using from problem 1 and it allows data files to emulate random-access
displays - I don't use any language implementation that I can set up two
memory-resident 20000-element arrays in.)
I reproduce the program below.
I started with a random initial vector produced by a tiny program - also
below. Using 1 as the random seed this led to the following output from the
main program:
-3.101685058674775063991546630859375
1.710273258650500480449806836167656148068123254767
0.560106265987016364553922305092009103311631866233
0.741538307396110261569760522731287644420481334112
0.727161301082217180143278466932462304558772157168
0.725501745009314329275732596978493267328561251751
0.725230798731860692293883213712712880613034610284
0.725126486420130425012250877920733508252708585739
0.72509323457494246776813920109013824975105118134
0.725082979071599326956147319028195277956039262739
0.72507978513079630216341650050591990474619363878
0.725078793237438412095625958539543698236163698988
0.725078485192179158058420718868522210900182511916
0.725078389445418392453736043330416495740873577237
0.725078359686792640646582053236955268581406362956
0.725078350438552911781751129445750050671394858183
0.725078347564408346610120634196508619714284794022
0.725078346671177581713561177548403984913441100435
0.725078346393577146486975902529735007528585737095
0.725078346307303727624091991942759104297833856823
0.725078346280491425023719407631506224988754601814
0.725078346272158616142422742375693986900626554449
0.725078346269568919530814464724619474374685630618
0.725078346268764085274769487940763571564603245775
After several minutes I had
0.725078346268401167468687719251160968869180594472
Just to check, I did another little program to multiply my solution vector
by A, and obtained (1,0,0,0,....,0) with accuracy for the first 34+ decimal
places.
Here are the programs:
Main program for the Gauss-Seidel iteration:
(after seeing M's program for #2 I was shamed into adding some
explanation:- still not up to her standard tho)
10 point 10
20 '
30 ' Here are two files I prepared earlier, containing, respectively,
40 ' the primes for the matrix diagonal and a file of 20000 random
50 ' elements which are the initial x_i values for the iteration.
60 '
70 open "primes" as file1
80 open "elements" as file2
85 '
90 repeat
100 for I=1 to 20000
110 S=0
120 '
130 ' S will become the sum of (a_i,j)*x_j for i<>j
140 ' We will only bother picking out the non-zero (a_i,j)'s
150 ' i.e. the 1's
160 '
170 for J=0 to 14 '=the powers of 2 for offsetting the 1's
180 P=2^J ' Why P? I can't remember. It's the offset.
190 L=I-P:U=I+P ' L and U are two candidate positions for 1's...
195 '...but we'd better check they actually lie in the matrix...
200 if L>0 then S=S+file2(L)
210 if U<20001 then S=S+file2(U)
220 next J
230 '
240 ' Now to divide by the prime number on the diagonal to get the
250 ' new approximation to x_i. But don't forget the very first equation
260 ' has a 1 on the other side.
270 '
280 if I=1 then S=S-1
290 file2(I)=-1*S/file1(I)
300 if I=1 then print file2(I)
310 next I
320 until 0 'the user has to use "CTRL-BREAK" when bored
Program to generate a random initial vector:
5 point 10
10 randomize 1
20 open "elements" as file2(20000) word 12
30 for I=1 to 20000:file2(I)=rnd:next I
Program to check the answer:
5 point 10
10 open "primes" as file1
20 open "elements" as file2
30 for I=1 to 20000
35 S=file1(I)*file2(I)
40 for J=0 to 14
50 P=2^J
60 L=I-P:U=I+P
70 if L>0 then S=S+file2(L)
80 if U<20001 then S=S+file2(U)
90 next J
100 if S>0 then print I,S:A$=input$(1)
110 next I
>PS Gauss Seidel doesn't mean a thing to me...
I know you wrote this before seeing my write-up of problem 7, but in case that didn't help...
Gauss Seidel is an iterative approach to solving simultaneous equations (I've only seen it for linear equations though in principle it could work in other cases too). Suppose the unknowns are x_1, x_2, ... , x_n and the kth equation in a set of n is
a_1*x_1 + a_2*x_2 + ... + a_n*x_n = b
This gets rearranged to
x_k = (b - a_1*x_1 - a_2*x_2 - ... - a_(k-1)*x_(k-1) - a_(k+1)*x_(k+1) -
... -a_n*x_n)/a_k
We do this with all n equations, so that each x_i is the subject of the ith equation. Then we take an initial guess at x_1, ... , x_n (actually we don't need an x_1) and plug them into the first equation which overwrites the initial choice of x_1, then we plug values into the second equation and so on up to the nth equation. And then we go round again. Either the values converge on the solution or they don't. Convergence is most likely if we start with a 'strong leading diagonal' i.e. ensure that we are dividing by big a_k's rather than small ones.
Regards,,
Brian.