POJ 3532 基尔霍夫电流定律 + 高斯消元

题意

传送门 POJ 3532

题解

基尔霍夫电流定律($KCL$)：在集总参数电路中，任一瞬间，流入或流出电路中任一节点的电流代数和恒等于零。若规定流出为正，则流入为负，$KCL$ 表示为
$$\sum i_{出}-\sum i_{入}=0$$

根据欧姆定律，计算电导值
$$G=1/R,I=U/R=GU$$

对于 $2$ 个节点间存在多个电阻 $R_{i,j}$，根据电阻并联公式计算阻值，设此 $2$ 节点 $i,j$ 间等效电阻为 $r_{i,j}$

$$1/r_{i,j}=\sum 1/R_{i,j}$$
设各节点电势为 $u_i$ ，则可以联立非节点 $1,N$ 的节点电流方程
$$\{\begin{array}{l}{\sum }_{节点i与节点1通过等效电阻直接连接}(u_i-u_1)/r_{1,j}=0\\ \dots \\ {\sum }_{节点i与节点N-1通过等效电阻直接连接}(u_i-u_{N-1})/r_{N-1,i}=0\end{array}$$

考虑到解的唯一性，取节点 $0$ 为零电势点，取节点 $N$ 为单位电势，则相当于在节点 $0,N$ 间加入了单位电压的电源
$$\{\begin{array}{l}{u}_0=0\\ u_N=1\end{array}$$

联立节点方程，高斯消元求解即可。考虑到节点 $0,N$ 间电势差为单位电势，则节点 $0,N$ 间等效电阻为
$$R=1/\sum i_{0,v}$$

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
#define min(a, b) (((a) < (b)) ? (a) : (b))
#define max(a, b) (((a) > (b)) ? (a) : (b))
#define abs(x) ((x) < 0 ? -(x) : (x))
#define INF 0x3f3f3f3f3f3f3f3f
#define delta 0.85
using namespace std;

#define maxn 105
const double EPS = 1e-8;
int N, M;
double R[maxn][maxn];

typedef vector<double> vec;
typedef vector<vec> mat;
// Ax = b
vec gauss_jordan(const mat &A, const vec &b)
{
    int n = A.size();
    mat B(n, vec(n + 1));
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = A[i][j];
        }
    }
    for (int i = 0; i < n; i++)
    {
        B[i][n] = b[i];
    }
    for (int i = 0; i < n; i++)
    {
        int pivot = i;
        for (int j = i; j < n; j++)
        {
            if (abs(B[j][i]) > abs(B[pivot][i]))
            {
                pivot = j;
            }
        }
        swap(B[i], B[pivot]);
        if (abs(B[i][i]) < EPS)
            return vec();
        for (int j = i + 1; j <= n; j++)
        {
            B[i][j] /= B[i][i];
        }
        for (int j = 0; j < n; j++)
        {
            if (i != j)
            {
                for (int k = i + 1; k <= n; k++)
                    B[j][k] -= B[j][i] * B[i][k];
            }
        }
    }
    vec x(n);
    for (int i = 0; i < n; i++)
        x[i] = B[i][n];
    return x;
}

void solve()
{
    mat A(N, vec(N, 0));
    vec b(N, 0);
    A[0][0] = A[N - 1][N - 1] = 1.0;
    b[N - 1] = 1.0;
    for (int i = 1; i < N - 1; i++)
    {
        for (int j = 0; j < N; j++)
        {
            if (abs(R[i][j]) < EPS)
                continue;
            A[i][i] -= 1.0 / R[i][j], A[i][j] += 1.0 / R[i][j];
        }
    }
    vec U = gauss_jordan(A, b);
    double I = 0.0;
    for (int i = 0; i < N; i++)
    {
        if (abs(R[0][i]) < EPS)
            continue;
        I += U[i] / R[0][i];
    }
    printf("%.2f", 1.0 / I);
}

int main()
{
    while (~scanf("%d%d", &N, &M))
    {
        for (int i = 0; i < N; i++)
        {
            fill(R[i], R[i] + N, 0);
        }
        for (int i = 0; i < M; i++)
        {
            int x, y, r;
            scanf("%d%d%d", &x, &y, &r);
            --x, --y;
            R[x][y] += 1.0 / r;
            R[y][x] += 1.0 / r;
        }
        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < N; j++)
            {
                if (abs(R[i][j]) < EPS)
                    continue;
                R[i][j] = 1 / R[i][j];
            }
        }
        solve();
    }
    return 0;
}