iconEuler Home

Large and Sparse Systems

In this notebook, we demonstrate how to use the cpx functions of Euler
for compressed large matrices.
>n=5; A=id(n)*2+diag(n,n,-1,1)+diag(n,n,1,1); fracprint(A);
        2         1         0         0         0 
        1         2         1         0         0 
        0         1         2         1         0 
        0         0         1         2         1 
        0         0         0         1         2 
This matrix can be compressed into a data structure containing
only a list of the elements different from 0.
>R=cpx(A)
Compressed 5x5 matrix
    1     1                    2
    1     2                    1
    2     1                    1
    2     2                    2
    2     3                    1
    3     2                    1
    3     3                    2
    3     4                    1
    4     3                    1
    4     4                    2
    4     5                    1
    5     4                    1
    5     5                    2
This really pays in the case of huge matrices.
>n=1000; A=id(n)*2+diag(n,n,-1,1)+diag(n,n,1,1);
>R=cpx(A);
Multiplication of two such matrices is much faster than multiplication
of the orginal matrices.
>cpxmult(R,R);
Using iterative methods, we can solve such large systems. For
a demonstration, we call the cpxfit function, which uses the conjugate
gradient method for R'Rx=R'b. It works even for 1000 equations
of 1000 unknowns.

We use for b the sum of the rows of A. The correct solution would
be x[j]=1 for all j.
>b=cpxsum(R); ...
 x=cpxfit(R,b); ...
 totalmax(cpxmult(R,x)-b)
2.27965202271e-010

Resistances on a Grid

Here is an example, where we compute the resistance of a rectangular
grid of size 50x50 from corner to corner. The grid has 2500 grid
points, which will be the number of variables of the linear systems.

Each connection has a resistance of 1. The flow of current between
two points is the difference in voltage divided by the resistance (1
in this case). The sum of the flows must be 0. So the voltage in each
point must be the average of the voltages of the points connected to
this point.
>n=50;
The function rectangleIncidenceX returns the incidence matrix
of the grid. It is of size n^2 times n^2 and has 1 in R[i,j],
if point i is connected with point j.
>R=rectangleIncidenceX(n,n);
We compute the numbers of edges going from each point, and divide
the rows by that number.
>c=cpxsum(R); R=cpxmultrows(1/c,R);
Now we set the diagonal to -1. The system R.x=0 is now satisfied
by all voltages x, such that each point is the average of all
neighboring points.
>R=cpxsetdiag(R,0,-1);
Now we add the condition that the first point has voltage 1, and
the last has voltage 0.
>R=cpxsetrow(R,1,0); R=cpxset(R,[1,1,1]); ...
 R=cpxsetrow(R,n^2,0); R=cpxset(R,[n^2,n^2,1]); ...
 b=zeros(n^2,1); b[1]=1;
We solve for R.x=b, using the conjugate gradient method with thin
matrices for 

19 - Introduction to Large Systems

It takes a few seconds for this large matrix.
>x=cpxfit(R,b);
Now plot the lines of equal voltage.
>plot2d(redim(x,n,n),a=1,b=50,c=1,d=50,hue=1,niveau="auto"):

19 - Introduction to Large Systems

In this case the total resistance is 1 over the current flow out
of the corner point. Flowing along each connection is a current
equal to the differences of voltages, since the resistance is
1. Note that the connected points have numbers 2 and n+1.
>1/(2*x[1]-x[2]-x[n+1])
5.05836762787
Here is another matrix example with a random matrix. We first
define an empty 1000x1000 matrix.
>H=cpxzeros(1000,1000);
Then we 10000 create random indices, which we want to set to a
random number later.
>ind=intrandom(10000,2,1000);
Now set those indices in H to 1.
>H=cpxset(H,ind|normal(rows(ind),1));
Set the diagonal of H to a large number to avoid non-regular matrices.
>H=cpxsetdiag(H,0,20);
Define b equal to the sum of the rows in H.
>b=cpxsum(H);
Solve the system. The solution should be all ones, so the differnce
should be 0.
>totalmax(cpxfit(H,b)-1)
0

Euler Home