[已解决]关于C语言和python语言数学运算的问题

xxfholic

问题背景：本人初学python，最近在完成学校里老师留的大作业，将祖传的C语言代码改成Python语言代码，在改的过程中发现两组代码产生的数据不同，想向fishc大佬们求助。
代码如下：

Python代码
import numpy as np
#from numba import jit

rho = 2850
E = 7.154*1e10
X = np.array([
    0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10
])
A = np.array([
    1.7*1e-4, 1.46*1e-4, 1.26*1e-4, 1.09*1e-4, 0.96*1e-4, 0.86*1e-4, 0.77*1e-4, 0.73*1e-4, 0.7*1e-4, 0.68*1e-4, 0.68*1e-4
])
I = np.array([
    2.79*1e-10, 2.12*1e-10, 1.57*1e-10, 1.08*1e-10, 0.84*1e-10, 0.61*1e-10, 0.45*1e-10, 0.37*1e-10, 0.32*1e-10, 0.30*1e-10, 0.30*1e-10
])
Y0js = np.array([
    0, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0
])
Y2js = np.array([
    0, -0.1, -0.3, -0.5, -0.4, -0.2, -0.1, 0.5, 0.8, 1.0
])
Y3js = np.array([
    0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0
])

Y0sj = np.zeros(10)
Y0xxx = np.zeros(10)
Y2 = np.zeros(10)
Y2sj = np.zeros(10)
Y3 = np.zeros(10)
Y3sj = np.zeros(10)

xz = 2.58

# midpoint function
def midpoint(x):
    return (x[:-1]+x[1:])/2.0

def _1st(Ab, Ib, Y0js, Y0sj, Y0xxx):
    dev = 10

    while dev >= 1e-6:
        
        _1st = np.zeros(10)
        _2nd = np.zeros(10)
        _3rd = np.zeros(10)
        _4th = np.zeros(10)
        
        #cycle 1
        a, b = 0 ,0 
        while a <= 9:
            b = a
            while b <= 8:
                _1st[a] +=  Ab[b+1] * Y0js[b+1] * 0.01
                b += 1
            a += 1
        
        #cycle 2
        a, b = 0 ,0 
        while a <= 9:
            b = a
            while b <= 8:
                _2nd[a] += _1st[b+1] * 0.01
                b += 1
            a += 1
       
        #cycle 3
        a, b = 0 ,0
        while a <= 9:
            b = 0
            while a >= b:
                _3rd[a] += _2nd[b] * 0.01 * (1/Ib[b])
                b += 1
            a += 1
        
        #cycle 4
        a, b = 0 ,0
        while a <= 9:
            b = 0
            while a >= b:
                _4th[a] += _3rd[b] * 0.01
                b += 1
            a += 1

        a = 0
        while a <= 9:
            Y0sj[a] = _4th[a] / _4th[9]
            Y0xxx[a] = _4th[a]
            a += 1
        
        dev = 0
        a = 0
        while a <= 9:
            dev += abs(Y0js[a] - Y0sj[a])
            a += 1
        
        a = 0
        while a <= 9:
            Y0js[a] = Y0sj[a]
            a += 1
    
    return 5000 * _4th[9] **-0.5 / 2 / np.pi

def _2nd(E, xz, rho, Ab, Ib, Y0xxx, Y2js, Y2, Y2sj):
    a21 = 10

    while abs(a21) >= 1e-8:
        b11 = 0
        c21 = 0
        a21 = 0

        i = 0
        while i <= 9:
            b11 += 2850 * Ab[i] * Y0xxx[i] **2 * 0.01
            i += 1
        
        i = 0
        while i <= 9:
            c21 += 2850 * Ab[i] * Y0xxx[i] * Y2js[i] * 0.01 
            i += 1

        a21 = c21 / b11

        i = 0
        while i <= 9:
            Y2[i] = Y2js[i] - a21 * Y0xxx[i]
            i += 1
        
        i = 0
        while i <= 9:
            Y2js[i] = Y2[i]
            i += 1
        
    i = 0
    while i <= 9:
        Y2sj[i] = Y2[i] / Y2[9]
        i += 1
    
    _1st = np.zeros(10)
    _2nd = np.zeros(10)
    _3rd = np.zeros(10)
    _4th = np.zeros(10)

    #cycle 1
    a, b = 0 ,0 
    while a <= 9:
        b = a
        while b <= 7:
            _1st[a] +=  Ab[b+2] * Y2sj[b+2] * 0.01
            b += 1
        a += 1
    
    #cycle 2
    a, b = 0 ,0 
    while a <= 9:
        b = a
        while b <= 8:
            _2nd[a] += _1st[b+1] * 0.01
            b += 1
        a += 1
    
    #cycle 3
    a, b = 0 ,0
    while a <= 9:
        b = 0
        while a >= b:
            _3rd[a] += _2nd[b] * 0.01 * (1/Ib[b])
            b += 1
        a += 1
    
    #cycle 4
    a, b = 0 ,0
    while a <= 9:
        b = 0
        while a >= b:
            _4th[a] += _3rd[b] * 0.01
            b += 1
        a += 1
    
    return (E / rho / _4th[9]) ** 0.5 / 2 / np.pi * xz

def _3rd(E, rho, Ab, Ib, Y0xxx, Y2js, Y2, Y2sj, Y3, Y3js):
    cvg = 10
    b11, b12, b21, b22, c31, c32, a31, a32 = 0, 0, 0, 0, 0, 0, 0, 10
    i = 10
    while cvg >= 1e-4: 
        
        #while i <= 9:
        #    Y3js[i] = Y3[i]
        #    i += 1

        b12, b21, b22, c31, c32, a31 = 0, 0, 0, 0, 0, 0

        i = 0
        while i <= 9:
            b11 += 2850 * Ab[i] * Y0xxx[i] **2 * 0.01
            i += 1
        
        i = 0
        while i <= 9:
                b12 += 2850 * Ab[i] * Y0xxx[i] * Y2[i] * 0.01
                i += 1 

        i = 0
        while i <= 9:
                b21 += 2850 * Ab[i] * Y0xxx[i] * Y2[i] * 0.01
                i += 1

        i = 0
        while i <= 9:
                #b22 += 2850 * Ab[i] * Y2[b] * Y2[i] * 0.01
            b22 += 2850 * Ab[i] * Y2[9] **2 * 0.01
            i += 1 

        i = 0
        while i <= 9:
                c31 += 2850 * Ab[i] * Y0xxx[i] * Y3js[i] * 0.01
                i += 1 

        i = 0
        while i <= 9:
                c32 += 2850 * Ab[i] * Y2[i] * Y3js[i] * 0.01
                i += 1
        
        a31 = c31 / b11
        a32 = c32 / b22
        
        i=0
        while i <= 9:
                Y3[i] = Y3js[i] - a31 * Y0xxx[i] - a32 * Y2[i]
                i += 1
        
        cvg = abs(Y3js[0] - Y3[0])

    _1st = np.zeros(10)
    _2nd = np.zeros(10)
    _3rd = np.zeros(10)
    _4th = np.zeros(10)

    #cycle 1
    a, b = 0 ,0 
    while a <= 9:
        b = a
        while b <= 7:
            _1st[a] +=  Ab[b+2] * Y3[b+2] * 0.01
            b += 1
        a += 1
    
    #cycle 2
    a, b = 0 ,0 
    while a <= 9:
        b = a
        while b <= 8:
            _2nd[a] += _1st[b+1] * 0.01
            b += 1
        a += 1
    
    #cycle 3
    a, b = 0 ,0
    while a <= 9:
        b = 0
        while a >= b:
            _3rd[a] += _2nd[b] * 0.01 * (1/Ib[b])
            b += 1
        a += 1
    
    #cycle 4
    a, b = 0 ,0
    while a <= 9:
        b = 0
        while a >= b:
            _4th[a] += _3rd[b] * 0.01
            b += 1
        a += 1
    
    return (E / rho / _4th[9]) ** 0.5 * cvg /np.pi

if __name__ == "__main__":
    
    Ab = midpoint(A)
    Ib = midpoint(I)

    omega_1 = _1st(Ab, Ib, Y0js, Y0sj, Y0xxx)
    omega_2 = _2nd(E, xz, rho, Ab, Ib, Y0xxx, Y2js, Y2, Y2sj)
    omega_3 = _3rd(E, rho, Ab, Ib, Y0xxx, Y2js, Y2, Y2sj, Y3, Y3js)
    print(omega_1)
    print(omega_2)
    print(omega_3)

C代码
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
int main(void)
{
        float rou=2850;
        float E=71540000000;
        float X[11]={0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10};
        float A[11]={0.00017,0.000146,0.000126,0.000109,0.000096,0.000086,0.000077,0.000073,0.00007,0.000068,0.000068};
        float I[11]={0.000000000279,0.000000000212,0.000000000157,0.000000000108,0.000000000084,0.000000000061,0.000000000045,0.000000000037,0.000000000032,0.000000000030,0.000000000030};
        float Ab[10];
        float Ib[10];
        
        int i=0;
        while(i<=9)
        {
                  Ab[i]=((A[i+1]+A[i])/2);
                  i=i+1;
        }
        
        i=0;
        while(i<=9)
        {
                  Ib[i]=((I[i+1]+I[i])/2);
                  i=i+1;
        }
        
        //第一个函数开始 
        float Y0js[10]={0,0.1,0.2,0.4,0.5,0.6,0.7,0.8,0.9,1};
        float Y0sj[10]={0,0,0,0,0,0,0,0,0,0};
        float Y0xxx[10]={0,0,0,0,0,0,0,0,0,0};
        float one[10]={0,0,0,0,0,0,0,0,0,0};
        float two[10]={0,0,0,0,0,0,0,0,0,0};
        float three[10]={0,0,0,0,0,0,0,0,0,0};
        float four[10]={0,0,0,0,0,0,0,0,0,0};
        float wucha=10;
        float xz=2.58; 
        int a=0;
        int b=0;
        int c=0;
        int d=0;
        int e=0;
        int f=0;
        int g=0;
        int h=0;
        int j=0;
        int k=0;
        int p=0;
        int q=0;
        while(wucha>=0.000001)                  //给定拟合精确度 
        {
                   q=0;
                   while(9>=q)
            {
                      one[q]=0;
                      two[q]=0;
                      three[q]=0;
                      four[q]=0;
                      q=q+1;
                   }                             //这步是给数组清零，千万不能忘！ 
                   
                a=0;
                   b=9;
            while(9>=a)
                {
                      b=a;
                      while(9>=b+1)
                        {
                          one[a]=one[a]+(Ab[b+1]*Y0js[b+1]*0.01);
                          b=b+1;
                      }    
                a=a+1;
                }                                         //第一重循环
                
                c=0;
            d=0;
                   while(9>=c)
            {
                      d=c;
                      while(9>=d+1)
                {
                              two[c]=two[c]+(one[d+1]*0.01);
                          d=d+1;
                }
                c=c+1;
            }                                 //第二重循环 

                e=0;
            f=0;
            while(9>=e)
            {
                        f=0;
                      while(e>=f)
                {
                          three[e]=three[e]+two[f]*0.01*(1/Ib[f]);
                          f=f+1;
                }
                e=e+1;
            }                                 //第三重循环  
   
                g=0;
            h=0;
            while(9>=g)
            {
                      h=0;
                      while(g>=h)
                {
                                 four[g]=four[g]+three[h]*0.01;
                                 h=h+1;
                }
                g=g+1;
            }                                 //第四重循环 
   
            k=0;
            while(9>=k)
            {
                      Y0sj[k]=four[k]/four[9];
                      Y0xxx[k]=four[k];
                          k=k+1;
              }                                //求出实际y0，以便和假设yo对比迭代 
   
            wucha=0;
            
            j=0;
            while(j<=9)
              {
                      wucha=wucha+fabs(Y0js[j]-Y0sj[j]);
                      j=j+1;
              }                                 //假设的y0与求出y0之间的误差  
   
            p=0;
            while(9>=p)
             {
                      Y0js[p]=Y0sj[p];
                p=p+1;
            }                                 //令实际值等于假设值，再次迭代运算 
   } 
        float omega=0;
        omega=5000*sqrt((1/four[9]));
        printf("一阶固有静频为：%.5fHZ\n",omega/(2*3.1415926));  
        
        //第二个函数开始 
        float Y2js[10]={0,-0.1,-0.3,-0.5,-0.4,-0.2,-0.1,0.5,0.8,1};
        float Y2[10]={0,0,0,0,0,0,0,0,0,0};
        float xiuzhen=4;
        float Y2sj[10]={0,0,0,0,0,0,0,0,0,0};
        float b11=0; 
        float C21=0;
        float a21=0;
        a21=10;
 
        while(fabs(a21)>=0.00000001)                  //给定拟合精确度 
        {
                   b11=0;
                   C21=0;
                   a21=0;
                
                i=0;
                while(i<=9)
                {
                        b11=b11+2850*Ab[i]*Y0xxx[i]*Y0xxx[i]*0.01;
                        i++; 
                } 
                
                i=0;
                while(i<=9)
                {
                        C21=C21+2850*Ab[i]*Y0xxx[i]*Y2js[i]*0.01;
                        i++; 
                }
                
                a21=C21/b11;

                i=0;
                while(i<=9)
                {
                        Y2[i]=Y2js[i]-a21*Y0xxx[i];
                        i++; 
                } 
 
                p=0;
            while(9>=p)
            {
                      Y2js[p]=Y2[p];
                p=p+1;
            }                                 //令实际值等于假设值，再次迭代运算 
   }
                
                i=0;
                while (i<=9)
                {
                        Y2sj[i]=Y2[i]/Y2[9];
                           i++;        
                   }                                     //归一化
                  
                q=0;
                   while(9>=q)
        {
                  one[q]=0;
                  two[q]=0;
                  three[q]=0;
                  four[q]=0;
                  q=q+1;
        }                             //这步是给数组清零，千万不能忘！ 
            
                a=0;
            b=0;
            while(9>=a)
            {
                      b=a;
                      while(9>=b+2)
                {
                          one[a]=one[a]+(Ab[b+2]*Y2sj[b+2]*0.01);
                          b=b+1;
                  }
                a=a+1;
              }                                 //第一重循环
 
            c=0;
            d=0;
            while(9>=c)
            {
                          d=c;
                      while(9>=d+1)
                {
                          two[c]=two[c]+(one[d+1]*0.01);
                          d=d+1;
                }
                c=c+1;
            }                                 //第二重循环 
  
            e=0;
            f=0;
            while(9>=e)
            {
                      f=0;
                      while(e>=f)
                {
                          three[e]=three[e]+two[f]*0.01*(1/Ib[f]);
                          f=f+1;
                }
                e=e+1;
            }                                 //第三重循环  
   
            g=0;
            h=0;
            while(9>=g)
            {
                      h=0;
                      while(g>=h)
                {
                                 four[g]=four[g]+three[h]*0.01;
                                 h=h+1;
                }
                g=g+1;
              }                                 //第四重循环
   
        omega=sqrt((E/rou)*(1/four[9]));
        omega=omega*xz;
        printf("二阶固有静频为：%.5fHZ\n",omega/(2*3.1415926));

        //第三个函数开始 
        float Y3js[10]={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1};
        float Y3[10]={0,0,0,0,0,0,0,0,0,0};
        float Y3sj[10]={0,0,0,0,0,0,0,0,0,0};

        b11=0;
        float b12=0; 
        float b21=0;
        float b22=0;
        float C31=0;
        float C32=0;
        float a31=0;
        float a32=0;
        a32=10;
        float shoulian=10;

        while(shoulian>=0.0001)                  //给定拟合精确度 
        {        
                while(i<=9)
                  {
                        Y3js[i]=Y3[i];
                        i++; 
                  } 
        
            b12=0;
            b21=0;
            b22=0;
            C31=0;
            C32=0;
            a31=0;
            a32=0;
                 
                i=0;
                while(i<=9)
                  {
                        b11=b11+2850*Ab[i]*Y0xxx[i]*Y0xxx[i]*0.01;
                        i++; 
                  } 
                
                i=0;
                while(i<=9)
                  {
                        b12=b12+2850*Ab[i]*Y0xxx[i]*Y2[i]*0.01;
                        i++; 
                  }

                i=0;
                while(i<=9)
                {
                        b21=b21+2850*Ab[i]*Y0xxx[i]*Y2[i]*0.01;
                        i++; 
                }  

                i=0;
                while(i<=9)
                {
                        b22=b22+2850*Ab[i]*Y2[b]*Y2[i]*0.01;
                        i++; 
                  }    

                i=0;
                while(i<=9)
                  {
                        C31=C31+2850*Ab[i]*Y0xxx[i]*Y3js[i]*0.01;
                        i++; 
                }

                i=0;
                while(i<=9)
                {
                        C32=C32+2850*Ab[i]*Y2[i]*Y3js[i]*0.01;
                        i++; 
                  }
                
                a31=C31/b11;  
                a32=C32/b22;   
                
                i=0;
                while(i<=9)
                {
                        Y3[i]=Y3js[i]-a31*Y0xxx[i]-a32*Y2[i];
                        i++; 
                } 
                  
                  shoulian=fabs(Y3js[0]-Y3[0]);
        }

    q=0;
           while(9>=q)
    {
        one[q]=0;
        two[q]=0;
        three[q]=0;
        four[q]=0;
        q=q+1;
    }                             //这步是给数组清零，千万不能忘！ 
    
        a=0;
    b=0;
    while(9>=a)
    {
              b=a;
              while(9>=b+2)
        {
                  one[a]=one[a]+(Ab[b+2]*Y3[b+2]*0.01);
                  b=b+1;
        }
        a=a+1;
    }                                 //第一重循环
 
    c=0;
    d=0;
    while(9>=c)
    {
              d=c;
              while(9>=d+1)
        {
                  two[c]=two[c]+(one[d+1]*0.01);
                  d=d+1;
        }
        c=c+1;
    }                                 //第二重循环 
  
    e=0;
    f=0;
    while(9>=e)
    {
              f=0;
              while(e>=f)
        {
                  three[e]=three[e]+two[f]*0.01*(1/Ib[f]);
                  f=f+1;
        }
        e=e+1;
        }                                 //第三重循环  
   
    g=0;
    h=0;
    while(9>=g)
    {
              h=0;
              while(g>=h)
        {
                 four[g]=four[g]+three[h]*0.01;
                 h=h+1;
        }
        g=g+1;
        }                                 //第四重循环
   
        omega=sqrt((E/rou)*(1/four[9]));
        omega=omega*xz; 
        omega=sqrt((E/rou)*(1/four[9]));
        omega=omega*xiuzhen;
        printf("三阶固有静频为：%.5fHZ\n",(2*omega)/(2*3.1415926));
        system("pause");
}

问题描述：两个程序第一个函数和第二个函数产生数据相同，第三个函数从变量b12开始计算有误。结算结果两者相差1e9。

最佳答案

月排行榜 / 总排行榜

耻思lhj

2020-4-15 23:09:10

所以你python的答案是什么？
我把你所有的float改成了double，因为我认为可能是精度的问题，下面是改了后c运行的结果
一阶固有静频为：406.98650HZ
二阶固有静频为：1795.49790HZ
三阶固有静频为：4594.51734HZ
都改为double：第九次循环b12=-4e-019
没改之前float:第九次循环b12=3e-010

跳转到最佳答案楼层

xxfholic

续帖：
与b12有关的数组：Ab，Y2，Y0xxx三组数值相同

python数组
Y2 = array([-0.00915524, 
-0.12828047,
 -0.3583434 ,
 -0.59971634,
 -0.55154928,
 -0.41236732, 
-0.37925495, 
 0.15134434,  
0.38194363, 
 0.51254293])

Y0xxx =
array([0.07180484,
 0.22180461,
 0.45758906, 
0.78207828, 
1.1886056 ,
1.66560334,
 2.19020504, 
2.73451689,
 3.27882874,
 3.82314059])

Ab= 
array([1.580e-04, 
1.360e-04, 
1.175e-04, 
1.025e-04,
 9.100e-05,
 8.150e-05,
7.500e-05, 
7.150e-05,
 6.900e-05,
 6.800e-05])

C数组
Ab = 
{0.00015800001,
 0.000136000002,
 0.0001175, 
0.000102500002, 
9.10000017e-005, 
8.14999948e-005,
 7.50000036e-005,
 7.1500006e-005, 
6.90000015e-005,
 6.80000012e-005}

Y0xxx = 
{0.0718048364,
 0.221804589, 
0.45758903,
 0.782078266, 
1.18860555, 
1.66560328, 
2.1902051, 
2.73451686, 
3.27882862, 
3.82314038}

Y2 = 
{-0.00915524177, 
-0.128280476, 
-0.358343422, 
-0.599716365, 
-0.551549315, 
-0.412367314, 
-0.379254967, 
0.151344314, 
0.381943643, 
0.512542963}

在执行循环时，前八次循环b12值相同，然而第九次循环结束后，两个代码b12的值就相差1e9，目前不知道为何，希望鱼C论坛的大佬们能够帮助我，比较着急，跪谢！

耻思lhj

所以你python的答案是什么？
我把你所有的float改成了double，因为我认为可能是精度的问题，下面是改了后c运行的结果
一阶固有静频为：406.98650HZ
二阶固有静频为：1795.49790HZ
三阶固有静频为：4594.51734HZ
都改为double：第九次循环b12=-4e-019
没改之前float:第九次循环b12=3e-010

xxfholic

耻思lhj 发表于 2020-4-15 23:09
所以你python的答案是什么？
我把你所有的float改成了double，因为我认为可能是精度的问题，下面是改了后c ...

谢谢带佬，原来是float和double的问题！
我这个python程序可能还有点问题，我再去看看。