#include <stdio.h> #include <math.h> #define numfre 100 // 函数：生成以10为底的等比数列 void logspace(double *y, double start, double end, int num) { double delta = (log10(end) - log10(start)) / (num - 1); int i; for (i = 0; i<num; i++) { y[i] = pow(10, log10(start) + i*delta); } } int main() { double fre[numfre], kn[5], Z0m[5], rho[5] = { 1, 1, 1, 1, 1 }, Zm, hm, Lm, rhoa[numfre]; const double wPI = 3.14159265357898E0; const double MU_0 = wPI * .0000004E0; int nf, n, m, nc = 5; double z[5] = { 0, 500, 1000, 1500, 2000 }; // 生成频率数组 logspace(fre, 0.0001, 10000, numfre); //循环频率 for (nf = 0; nf<numfre; nf++) { //计算每一层的本征阻抗 for (n = 0; n<nc; n++) { kn[n] = sqrt(-n*2 * wPI*fre[nf] * MU_0 / rho[n]); Z0m[n] = -n*2 * wPI*fre[nf] * MU_0 / kn[n]; } Zm = Z0m[nc - 1]; //最后一层本征阻抗 for (m = nc - 2; m >= 0; m--) { hm = z[m + 1] - z[m]; Lm = (Z0m[m] - Zm) / (Z0m[m] + Zm); Zm = Z0m[m] * (1 - Lm*exp(-2 * kn[m] * hm)) / (1 + Lm*exp(-2 * kn[m] * hm)); } rhoa[nf] = (1 / (2 * wPI*fre[nf] * MU_0))*pow(fabs(Zm), 2); } //打印结果 for (nf = 0; nf<numfre; nf++) { printf("频率%fHz时，电阻率为：%f\n", fre[nf], rhoa[nf]); } return 0; }上述代码运行结果错误，导致算出的电阻率极大，找出原因并修改

在第16行循环频率时，应该使用循环变量nf而不是n，因为变量n在第19行的循环中已经被使用了。修改后的代码如下：

#include <stdio.h> 
#include <math.h>  

#define numfre 100  

void logspace(double *y, double start, double end, int num) {
    double delta = (log10(end) - log10(start)) / (num - 1);
    int i;
    for (i = 0; i<num; i++) {
        y[i] = pow(10, log10(start) + i*delta);
    }
}

int main() {
    double fre[numfre], kn[5], Z0m[5], rho[5] = { 1, 1, 1, 1, 1 }, Zm, hm, Lm, rhoa[numfre];
    const double wPI = 3.14159265357898E0;
    const double MU_0 = wPI * .0000004E0;
    int nf, n, m, nc = 5;
    double z[5] = { 0, 500, 1000, 1500, 2000 };

    logspace(fre, 0.0001, 10000, numfre);

    for (nf = 0; nf<numfre; nf++)
    {
        for (n = 0; n<nc; n++)
        {
            kn[n] = sqrt(-n*2 * wPI*fre[nf] * MU_0 / rho[n]);
            Z0m[n] = -n*2 * wPI*fre[nf] * MU_0 / kn[n];
        }

        Zm = Z0m[nc - 1];

        for (m = nc - 2; m >= 0; m--)
        {
            hm = z[m + 1] - z[m];
            Lm = (Z0m[m] - Zm) / (Z0m[m] + Zm);
            Zm = Z0m[m] * (1 - Lm*exp(-2 * kn[m] * hm)) / (1 + Lm*exp(-2 * kn[m] * hm));
        }

        rhoa[nf] = (1 / (2 * wPI*fre[nf] * MU_0))*pow(fabs(Zm), 2);
    }

    for (nf = 0; nf<numfre; nf++)
    {
        printf("频率%fHz时，电阻率为：%f\n", fre[nf], rhoa[nf]);
    }

    return 0;
}