2012-10-26 28 views
10

मैं एक छोटे रे ट्रैसर लिखने के लिए पाइथन/numpy/scipy के साथ काम कर रहा हूं। सतहों को सामान्य विमान के ऊपर ऊंचाई प्रदान करने वाले द्वि-आयामी कार्यों के रूप में मॉडलिंग किया जाता है। मैंने एक चर के साथ एक समारोह की जड़ को खोजने के लिए किरण और सतह के बीच चौराहे के बिंदु को खोजने की समस्या को कम कर दिया। कार्य निरंतर और लगातार भिन्न हैं।एक चर के साथ बड़ी संख्या में कार्यों की जड़ों को ढूंढना

क्या scipy रूट खोजकर्ताओं (और शायद कई प्रक्रियाओं का उपयोग कर) का उपयोग करके सभी कार्यों पर लूपिंग करने से कहीं अधिक कुशलता से ऐसा करने का कोई तरीका है?

संपादित करें: फ़ंक्शन एक रेखीय फ़ंक्शन के बीच अंतर है जो किरण और सतह फ़ंक्शन का प्रतिनिधित्व करता है, जो चौराहे के एक विमान के लिए बाध्य है।

+2

फ़ंक्शन क्या है? क्या यह संभव है कि इसका एक विश्लेषणात्मक समाधान हो? – mgilson

+1

सतह कार्यों को मनमाने ढंग से चुना जा सकता है, मैं इसे लचीला होना चाहता हूं। एक विशिष्ट कार्य (यानी चेबिशहेव बहुपदों की एक सुपरपोजिशन) के लिए एक विश्लेषणात्मक समाधान मौजूद है, फिर भी इसमें कई पैरामीटर शामिल हो सकते हैं। विशिष्ट सतहों के लिए रैखिक समीकरणों की एक प्रणाली को हल करके चौराहे ढूंढना संभव होना चाहिए। – mikebravo

+0

रे/विमान, किरण/क्षेत्र, किरण/त्रिभुज चौराहे खोजने के मानक तरीके हैं। क्या आप अपनी सतह को त्रिभुज जाल के रूप में मॉडल कर सकते हैं? एक विश्लेषणात्मक समाधान या आपके सतह समारोह के लिए एक ज्यामितीय अनुमान के बिना, मुझे नहीं पता कि कार्यों के माध्यम से केवल क्रैंकिंग की तुलना में एक और अधिक प्रभावी तरीका है। – engineerC

उत्तर

3

निम्नलिखित उदाहरण फ़ंक्शन x ** (ए + 1) - बी (सभी ए और बी के साथ सभी) की बीमारियों की विधि का उपयोग करते हुए समानांतर में जड़ों की गणना करता है। यहाँ लगभग 12 सेकंड लेता है।

import numpy 

def F(x, a, b): 
    return numpy.power(x, a+1.0) - b 

N = 1000000 

a = numpy.random.rand(N) 
b = numpy.random.rand(N) 

x0 = numpy.zeros(N) 
x1 = numpy.ones(N) * 1000.0 

max_step = 100 
for step in range(max_step): 
    x_mid = (x0 + x1)/2.0 
    F0 = F(x0, a, b) 
    F1 = F(x1, a, b) 
    F_mid = F(x_mid, a, b) 
    x0 = numpy.where(numpy.sign(F_mid) == numpy.sign(F0), x_mid, x0) 
    x1 = numpy.where(numpy.sign(F_mid) == numpy.sign(F1), x_mid, x1) 
    error_max = numpy.amax(numpy.abs(x1 - x0)) 
    print "step=%d error max=%f" % (step, error_max) 
    if error_max < 1e-6: break 

मूल विचार बस एक रूट खोजक के सभी सामान्य कदम एक समारोह है कि मानकों के चर का एक वेक्टर और बराबर वेक्टर (ओं) पर मूल्यांकन किया जा सकता का उपयोग कर, चर का एक वेक्टर पर समानांतर में चलाने के लिए है जो व्यक्तिगत घटक कार्यों को परिभाषित करता है। कंडीशनल को मास्क और numpy के संयोजन से बदल दिया जाता है। कहीं()। यह तब तक जारी रख सकता है जब तक कि सभी जड़ों को आवश्यक परिशुद्धता, या वैकल्पिक रूप से पर्याप्त जड़ें नहीं मिल जाती हैं, जब तक कि उन्हें समस्या से दूर करने के लिए लायक नहीं है और उन जड़ों को छोड़कर एक छोटी समस्या के साथ जारी रहें।

जिन कार्यों को मैंने हल करने के लिए चुना है वे मनमानी हैं, लेकिन यदि कार्य अच्छी तरह से व्यवहार किए जाते हैं तो इससे मदद मिलती है; इस मामले में परिवार में सभी कार्य monotonic हैं और वास्तव में एक सकारात्मक रूट है। इसके अतिरिक्त, बिसेक्शन विधि के लिए हमें चर के लिए अनुमानों की आवश्यकता होती है जो फ़ंक्शन के विभिन्न संकेत देते हैं, और यहां भी साथ आने के लिए काफी आसान होता है (x0 और x1 के प्रारंभिक मान)।

उपरोक्त कोड शायद सबसे सरल रूट खोजक (उधार) का उपयोग करता है, लेकिन उसी तकनीक को न्यूटन-रैफसन, रिडर के इत्यादि पर आसानी से लागू किया जा सकता है। कम सशर्त रूट रूट विधि में हैं, बेहतर अनुकूल है इसके लिए। हालांकि, आपको इच्छित एल्गोरिदम को फिर से लागू करना होगा, मौजूदा लाइब्रेरी रूट खोजक फ़ंक्शन का उपयोग करने का कोई तरीका नहीं है।

उपरोक्त कोड स्निपेट दिमाग में स्पष्टता के साथ लिखा गया है, गति नहीं। केवल एक बार के बजाय 3 बार की यात्रा प्रति समारोह का मूल्यांकन कुछ गणना की पुनरावृत्ति से बचना, विशेष रूप से, इस, 9 सेकंड तक गति इस प्रकार है:

... 
F0 = F(x0, a, b) 
F1 = F(x1, a, b) 

max_step = 100 
for step in range(max_step): 
    x_mid = (x0 + x1)/2.0 
    F_mid = F(x_mid, a, b) 
    mask0 = numpy.sign(F_mid) == numpy.sign(F0) 
    mask1 = numpy.sign(F_mid) == numpy.sign(F1) 
    x0 = numpy.where(mask0, x_mid, x0) 
    x1 = numpy.where(mask1, x_mid, x1) 
    F0 = numpy.where(mask0, F_mid, F0) 
    F1 = numpy.where(mask1, F_mid, F1) 
... 

तुलना के लिए, scipy.bisect() का उपयोग कर एक को खोजने के लिए एक समय में रूट ~ 94 सेकंड लेता है:

for i in range(N): 
    x_root = scipy.optimize.bisect(lambda x: F(x, a[i], b[i]), x0[i], x1[i], xtol=1e-6) 
+0

मैं अभी एक scipy रूट खोजक का उपयोग कर रहा हूँ ... और मैं थोड़ा डंबफॉल्ड हूँ। यह कभी मेरे लिए कभी नहीं हुआ कि एल्गोरिदम का एक taylored कार्यान्वयन पुस्तकालय की तुलना में उस नौकरी तेजी से कर सकता है। यह वहां परिमाण का एक क्रम है। और मैंने पहले से ही अपने सभी वेक्टर बीजगणित को numpy प्रसारण का उपयोग कर किया है ... अगली बार जब मैं अच्छी तरह से प्रलेखित एल्गोरिदम के पुस्तकालय कार्यान्वयन का उपयोग करता हूं, तो मुझे इस बारे में सोचना होगा! – mikebravo