RosettaCodeData/Task/Resistor-mesh/Euler-Math-Toolbox/resistor-mesh-2.euler

45 lines
1.1 KiB
Plaintext

function makeRectangleX (n:index,m:index)
## Make the incidence matrix of a rectangle grid in compact form.
## see: makeRectangleIncidence
K=zeros(n*(m-1)+m*(n-1),3);
k=1;
for i=1 to n;
for j=1 to m-1;
K[k,1]=(i-1)*m+j; K[k,2]=(i-1)*m+j+1; K[k,3]=1;
k=k+1;
end;
end;
for i=1 to n-1;
for j=1 to m;
K[k,1]=(i-1)*m+j; K[k,2]=i*m+j; K[k,3]=1;
k=k+1;
end;
end;
H=cpxzeros([n*m,n*m]);
H=cpxset(H,K);
H=cpxset(H,K[:,[2,1,3]]);
return H;
endfunction
function solvePotentialX (A:cpx, i:index ,j:index)
## Solve the potential problem of resistance in a graph.
## This functions uses the conjugate gradient method.
## A is a compressed incidence matrix.
## Return the potential u for the nodes in A,
## such that u[i]=1, u[j]=-1, and the flow
## to each knot is equal to the flow from the knot,
## and the flow from i to j is (u[i]-u[j])*A[i,j].
## see: makeIncidence
n=size(A)[1];
b=ones(n,1); f=-cpxmult(A,b);
h=1:n; B=cpxset(A,h'|h'|f);
B=cpxset(B,i|h'|0);
B=cpxset(B,[i,i,1]);
B=cpxset(B,j|h'|0);
B=cpxset(B,[j,j,1]);
v=zeros(n,1); v[i]=1; v[j]=-1;
u=cpxfit(B,v);
f=(-f[i])*u[i]-cpxmult(A,u)[i];
return {u,2/f}
endfunction