2016-03-14 81 views
0

我模拟一个研究项目的点源,并使用我在这里的参数,我的代码应该输出值的顺序1.something +或 - 0.00something * j,我得到的值非常小。像3.8 * 10^-127一样。我跑到我的代码的输出溢出错误

import numpy as np 
from math import pi 

wvl = 1e-6 
R = 50e3 
k = (2 * pi)/wvl 
N = 10.0 
r0 = 0.1 
D2 = 0.5 
DROI = 4.0 * D2 
D1 = (wvl * R)/DROI 
x = np.linspace(-N/2, (N/2)-1, N) 
y = np.linspace((N/2)-1, -N/2, N) 
x1, y1 = np.meshgrid(x, y) 

rho = np.sqrt(x1**2 + y1**2) 

r1 = np.arctan2(y1, x1) 

pointOne = np.exp(-1j*k/(2*R) * (r1**2)) 
pointTwo = pointOne/((D1**2) * np.sinc(x1/D1) * np.sinc(y1/D1) * np.exp(-(r1/(4.0*D1))**2)) 

print pointTwo 

任何帮助,这将不胜感激!

回答

0

您是略有不正确的时候你说

像3.8 * 10^-127

我看到有exp(-986.96044011)这是......我不知道! :) 例如,exp(-745)等于4.9406564584124654e-324这是比你说的小。在numpy中最精确的类型是float64,上面的值是它可以采取的最小数量(https://en.wikipedia.org/wiki/Double-precision_floating-point_format

Soooo,我会说这是不可能用numpy来计算的。

我认为有没有其他的方法来计算这个除了使用decimal模块和计算与固定点:

>> decimal.Decimal(-986.96044011).exp() 
Decimal('2.336291362873238577842336912E-429') 

UPDATE

from decimal import Decimal 

factor1 = (D1**2) * np.sinc(x1/D1) * np.sinc(y1/D1) 
factor2 = -(r1/(4.0*D1))**2 

for i in range(10): 
    for j in range(10): 
     denom = Decimal(factor1[i][j]) * Decimal(factor2[i][j]).exp() 
     real = Decimal(pointOne[i][j].real)/denom 
     imag = Decimal(pointOne[i][j].imag)/denom 
     print('{} {:+}j'.format(real, imag)) 
    print 

我不会复制/粘贴整个输出,但是对于具有错误(i = 3)的块,其由numpy输出,如:

[    -inf    +infj    inf    +infj 
       -inf    +infj    -inf    +infj 
    -1.33430914e+277 +1.39038774e+276j 2.71287770e+126 -5.24351460e+126j 
    3.17716647e+062 -5.65262280e+062j 1.34513017e+045 -1.84387526e+045j 
    -3.44573750e+040 +7.75538991e+039j -3.43926477e+038 +2.50357768e+038j] 

该代码输出如下:

2.581492951664502006187030185E+411 -5.762227636689223315330742967E+411j 
2.122343066940105718447024076E+400 +1.467312608787517020408471212E+400j 
2.397507101683903303475691738E+381 -2.280696798129953463268938559E+380j 
-6.666386049218522221620946661E+346 +2.888168970661726622146331973E+347j 
-1.334309136463760612094105714E+277 +1.390387738631591203923444430E+276j 
2.712877702928865987066875589E+126 -5.243514596467982559177136476E+126j 
3.177166469663657952353553473E+62 -5.652622800145748169241881476E+62j 
1.345130165198935427036822984E+45 -1.843875257401926846810358724E+45j 
-3.445737499727613189278598753E+40 +7.755389906593324611267619237E+39j 
-3.439264771527389293832980697E+38 +2.503577682325662775720694139E+38j 

正如你所看到的,正确的尾部相匹配。 不知道它是否可以接受您的解决方案,因为您无法再使用该数据与numpy一起使用。 decimal模块缺乏对复杂数字的支持,因此您必须手动实现它或者为此目的找到一些第三方库

+0

谢谢哀思!你对exp函数的评论让我意识到了一些事情,并且在改变代码后我能够得到它来输出我想要的数据! –

+0

我恳请您阅读本帮助文章:http://stackoverflow.com/help/someone-answers – Grief

+1

您的回复非常有帮助,谢谢!我意识到我正在使用rho和r1混合,并在代码中使用错误的代码,这导致了令人难以置信的小数值! –