Поиск корней большого числа функций с одной переменной

Я работаю с Python / numpy / scipy, чтобы написать небольшой трассировщик лучей. Поверхности моделируются как двумерные функции, дающие высоту над нормальной плоскостью. Я уменьшил задачу нахождения точки пересечения луча и поверхности с нахождением корня функции с одной переменной. Функции непрерывны и непрерывно дифференцируемы.

Есть ли способ сделать это более эффективно, чем просто перебирать все функции, используя scipy root finders (и, возможно, используя несколько процессов)?

Изменить: функции представляют собой разницу между линейной функцией, представляющей луч, и функцией поверхности, ограниченной плоскостью пересечения.

One Solution collect form web for “Поиск корней большого числа функций с одной переменной”

В следующем примере показано вычисление корней для 1 миллиона копий функции x ** (a + 1) – b (все с разными a и b) параллельно с использованием метода деления пополам. Здесь занимает около 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.where (). Это может продолжаться до тех пор, пока все корни не будут найдены с требуемой точностью или поочередно до тех пор, пока не будет найдено достаточное количество корней, что стоит удалить их из проблемы и продолжить меньшую проблему, исключающую эти корни.

Функции, которые я выбрал для решения, являются произвольными, но это помогает, если функции хорошо себя ведут; в этом случае все функции в семействе монотонны и имеют ровно один положительный корень. Кроме того, для метода деления пополам нам нужны догадки для переменной, которые дают разные знаки функции, и здесь также довольно легко придумать (начальные значения x0 и x1).

В приведенном выше коде используется, пожалуй, самый простой корневой искатель (bisection), но тот же метод можно легко применить к Ньютону-Рафсону, Риддеру и т. Д. Чем меньше условностей существует в методе поиска корней, тем лучше ему подходит. Тем не менее, вам нужно будет переопределить любой алгоритм, который вам нужен, нет возможности напрямую использовать существующую функцию корневого поиска библиотеки.

Вышеприведенный фрагмент кода написан с ясностью, а не скоростью. Избегая повторения некоторых вычислений, в частности, оценивая функцию только один раз на итерацию, а не в 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) 
Python - лучший язык программирования в мире.