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
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
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"):
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