2016-03-14 39 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混合,並在代碼中使用錯誤的代碼,這導致了令人難以置信的小數值! –

相關問題