关于C语言和python语言数学运算的问题
问题背景：本人初学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)/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 +=Ab * Y0js * 0.01
b += 1
a += 1
#cycle 2
a, b = 0 ,0
while a <= 9:
b = a
while b <= 8:
_2nd += _1st * 0.01
b += 1
a += 1
#cycle 3
a, b = 0 ,0
while a <= 9:
b = 0
while a >= b:
_3rd += _2nd * 0.01 * (1/Ib)
b += 1
a += 1
#cycle 4
a, b = 0 ,0
while a <= 9:
b = 0
while a >= b:
_4th += _3rd * 0.01
b += 1
a += 1
a = 0
while a <= 9:
Y0sj = _4th / _4th
Y0xxx = _4th
a += 1
dev = 0
a = 0
while a <= 9:
dev += abs(Y0js - Y0sj)
a += 1
a = 0
while a <= 9:
Y0js = Y0sj
a += 1
return 5000 * _4th **-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 * Y0xxx **2 * 0.01
i += 1
i = 0
while i <= 9:
c21 += 2850 * Ab * Y0xxx * Y2js * 0.01
i += 1
a21 = c21 / b11
i = 0
while i <= 9:
Y2 = Y2js - a21 * Y0xxx
i += 1
i = 0
while i <= 9:
Y2js = Y2
i += 1
i = 0
while i <= 9:
Y2sj = Y2 / Y2
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 +=Ab * Y2sj * 0.01
b += 1
a += 1
#cycle 2
a, b = 0 ,0
while a <= 9:
b = a
while b <= 8:
_2nd += _1st * 0.01
b += 1
a += 1
#cycle 3
a, b = 0 ,0
while a <= 9:
b = 0
while a >= b:
_3rd += _2nd * 0.01 * (1/Ib)
b += 1
a += 1
#cycle 4
a, b = 0 ,0
while a <= 9:
b = 0
while a >= b:
_4th += _3rd * 0.01
b += 1
a += 1
return (E / rho / _4th) ** 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 = Y3
# i += 1
b12, b21, b22, c31, c32, a31 = 0, 0, 0, 0, 0, 0
i = 0
while i <= 9:
b11 += 2850 * Ab * Y0xxx **2 * 0.01
i += 1
i = 0
while i <= 9:
b12 += 2850 * Ab * Y0xxx * Y2 * 0.01
i += 1
i = 0
while i <= 9:
b21 += 2850 * Ab * Y0xxx * Y2 * 0.01
i += 1
i = 0
while i <= 9:
#b22 += 2850 * Ab * Y2 * Y2 * 0.01
b22 += 2850 * Ab * Y2 **2 * 0.01
i += 1
i = 0
while i <= 9:
c31 += 2850 * Ab * Y0xxx * Y3js * 0.01
i += 1
i = 0
while i <= 9:
c32 += 2850 * Ab * Y2 * Y3js * 0.01
i += 1
a31 = c31 / b11
a32 = c32 / b22
i=0
while i <= 9:
Y3 = Y3js - a31 * Y0xxx - a32 * Y2
i += 1
cvg = abs(Y3js - Y3)
_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 +=Ab * Y3 * 0.01
b += 1
a += 1
#cycle 2
a, b = 0 ,0
while a <= 9:
b = a
while b <= 8:
_2nd += _1st * 0.01
b += 1
a += 1
#cycle 3
a, b = 0 ,0
while a <= 9:
b = 0
while a >= b:
_3rd += _2nd * 0.01 * (1/Ib)
b += 1
a += 1
#cycle 4
a, b = 0 ,0
while a <= 9:
b = 0
while a >= b:
_4th += _3rd * 0.01
b += 1
a += 1
return (E / rho / _4th) ** 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={0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10};
float A={0.00017,0.000146,0.000126,0.000109,0.000096,0.000086,0.000077,0.000073,0.00007,0.000068,0.000068};
float I={0.000000000279,0.000000000212,0.000000000157,0.000000000108,0.000000000084,0.000000000061,0.000000000045,0.000000000037,0.000000000032,0.000000000030,0.000000000030};
float Ab;
float Ib;
int i=0;
while(i<=9)
{
Ab=((A+A)/2);
i=i+1;
}
i=0;
while(i<=9)
{
Ib=((I+I)/2);
i=i+1;
}
//第一个函数开始
float Y0js={0,0.1,0.2,0.4,0.5,0.6,0.7,0.8,0.9,1};
float Y0sj={0,0,0,0,0,0,0,0,0,0};
float Y0xxx={0,0,0,0,0,0,0,0,0,0};
float one={0,0,0,0,0,0,0,0,0,0};
float two={0,0,0,0,0,0,0,0,0,0};
float three={0,0,0,0,0,0,0,0,0,0};
float four={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=0;
two=0;
three=0;
four=0;
q=q+1;
} //这步是给数组清零，千万不能忘！
a=0;
b=9;
while(9>=a)
{
b=a;
while(9>=b+1)
{
one=one+(Ab*Y0js*0.01);
b=b+1;
}
a=a+1;
} //第一重循环
c=0;
d=0;
while(9>=c)
{
d=c;
while(9>=d+1)
{
two=two+(one*0.01);
d=d+1;
}
c=c+1;
} //第二重循环
e=0;
f=0;
while(9>=e)
{
f=0;
while(e>=f)
{
three=three+two*0.01*(1/Ib);
f=f+1;
}
e=e+1;
} //第三重循环
g=0;
h=0;
while(9>=g)
{
h=0;
while(g>=h)
{
four=four+three*0.01;
h=h+1;
}
g=g+1;
} //第四重循环
k=0;
while(9>=k)
{
Y0sj=four/four;
Y0xxx=four;
k=k+1;
} //求出实际y0，以便和假设yo对比迭代
wucha=0;
j=0;
while(j<=9)
{
wucha=wucha+fabs(Y0js-Y0sj);
j=j+1;
} //假设的y0与求出y0之间的误差
p=0;
while(9>=p)
{
Y0js=Y0sj;
p=p+1;
} //令实际值等于假设值，再次迭代运算
}
float omega=0;
omega=5000*sqrt((1/four));
printf("一阶固有静频为：%.5fHZ\n",omega/(2*3.1415926));
//第二个函数开始
float Y2js={0,-0.1,-0.3,-0.5,-0.4,-0.2,-0.1,0.5,0.8,1};
float Y2={0,0,0,0,0,0,0,0,0,0};
float xiuzhen=4;
float Y2sj={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*Y0xxx*Y0xxx*0.01;
i++;
}
i=0;
while(i<=9)
{
C21=C21+2850*Ab*Y0xxx*Y2js*0.01;
i++;
}
a21=C21/b11;
i=0;
while(i<=9)
{
Y2=Y2js-a21*Y0xxx;
i++;
}
p=0;
while(9>=p)
{
Y2js=Y2;
p=p+1;
} //令实际值等于假设值，再次迭代运算
}
i=0;
while (i<=9)
{
Y2sj=Y2/Y2;
i++;
} //归一化
q=0;
while(9>=q)
{
one=0;
two=0;
three=0;
four=0;
q=q+1;
} //这步是给数组清零，千万不能忘！
a=0;
b=0;
while(9>=a)
{
b=a;
while(9>=b+2)
{
one=one+(Ab*Y2sj*0.01);
b=b+1;
}
a=a+1;
} //第一重循环
c=0;
d=0;
while(9>=c)
{
d=c;
while(9>=d+1)
{
two=two+(one*0.01);
d=d+1;
}
c=c+1;
} //第二重循环
e=0;
f=0;
while(9>=e)
{
f=0;
while(e>=f)
{
three=three+two*0.01*(1/Ib);
f=f+1;
}
e=e+1;
} //第三重循环
g=0;
h=0;
while(9>=g)
{
h=0;
while(g>=h)
{
four=four+three*0.01;
h=h+1;
}
g=g+1;
} //第四重循环
omega=sqrt((E/rou)*(1/four));
omega=omega*xz;
printf("二阶固有静频为：%.5fHZ\n",omega/(2*3.1415926));
//第三个函数开始
float Y3js={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1};
float Y3={0,0,0,0,0,0,0,0,0,0};
float Y3sj={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=Y3;
i++;
}
b12=0;
b21=0;
b22=0;
C31=0;
C32=0;
a31=0;
a32=0;
i=0;
while(i<=9)
{
b11=b11+2850*Ab*Y0xxx*Y0xxx*0.01;
i++;
}
i=0;
while(i<=9)
{
b12=b12+2850*Ab*Y0xxx*Y2*0.01;
i++;
}
i=0;
while(i<=9)
{
b21=b21+2850*Ab*Y0xxx*Y2*0.01;
i++;
}
i=0;
while(i<=9)
{
b22=b22+2850*Ab*Y2*Y2*0.01;
i++;
}
i=0;
while(i<=9)
{
C31=C31+2850*Ab*Y0xxx*Y3js*0.01;
i++;
}
i=0;
while(i<=9)
{
C32=C32+2850*Ab*Y2*Y3js*0.01;
i++;
}
a31=C31/b11;
a32=C32/b22;
i=0;
while(i<=9)
{
Y3=Y3js-a31*Y0xxx-a32*Y2;
i++;
}
shoulian=fabs(Y3js-Y3);
}
q=0;
while(9>=q)
{
one=0;
two=0;
three=0;
four=0;
q=q+1;
} //这步是给数组清零，千万不能忘！
a=0;
b=0;
while(9>=a)
{
b=a;
while(9>=b+2)
{
one=one+(Ab*Y3*0.01);
b=b+1;
}
a=a+1;
} //第一重循环
c=0;
d=0;
while(9>=c)
{
d=c;
while(9>=d+1)
{
two=two+(one*0.01);
d=d+1;
}
c=c+1;
} //第二重循环
e=0;
f=0;
while(9>=e)
{
f=0;
while(e>=f)
{
three=three+two*0.01*(1/Ib);
f=f+1;
}
e=e+1;
} //第三重循环
g=0;
h=0;
while(9>=g)
{
h=0;
while(g>=h)
{
four=four+three*0.01;
h=h+1;
}
g=g+1;
} //第四重循环
omega=sqrt((E/rou)*(1/four));
omega=omega*xz;
omega=sqrt((E/rou)*(1/four));
omega=omega*xiuzhen;
printf("三阶固有静频为：%.5fHZ\n",(2*omega)/(2*3.1415926));
system("pause");
}
问题描述：两个程序第一个函数和第二个函数产生数据相同，第三个函数从变量b12开始计算有误。结算结果两者相差1e9。 续帖：
与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 于 2020-4-15 23:33 编辑
所以你python的答案是什么？
我把你所有的float改成了double，因为我认为可能是精度的问题，下面是改了后c运行的结果
一阶固有静频为：406.98650HZ
二阶固有静频为：1795.49790HZ
三阶固有静频为：4594.51734HZ
都改为double：第九次循环b12=-4e-019
没改之前float:第九次循环b12=3e-010
耻思lhj 发表于 2020-4-15 23:09
所以你python的答案是什么？
我把你所有的float改成了double，因为我认为可能是精度的问题，下面是改了后c ...
谢谢带佬，原来是float和double的问题！
我这个python程序可能还有点问题，我再去看看。
页:
[1]