Найти корень (предел интегрирования) в численном интегрировании

Question

Я пытаюсь переписать код Mathematica для построения равнозаселены кольца:

Nr = 5; (*radial modes*) 


DF0[JJ_] := Exp[-JJ]; (*distribution function of long action*) 
Jmax = 20; (* max action for numerical cuts*) 
CF = NIntegrate[DF0[II], {II, 0, Jmax}]; 
DF[JJ_] := DF0[JJ]/CF; 

bJ = Array[0, Nr + 1]; 
bJ[[1]] = 0.; 
bJ[[Nr + 1]] = Jmax; 

For[ir = 2, ir < Nr + 1, ir++, 
    bb = bJ[[ir - 1]]; 
    bJ[[ir]] = 
    Jroot /. 
    FindRoot[ 
    Integrate[DF[JJ], {JJ, bb, Jroot}] == 1/Nr, {Jroot, bb + 1/(Nr DF[bb])}]; 
    ]; 

bJ = Re[bJ]; 


(* Finding actions aJ of the thin rings: *) 

aJ = Array[0, Nr]; 
For[ir = 1, ir < Nr + 1, ir++, 
    aJ[[ir]] = NIntegrate[J Nr DF[J], {J, bJ[[ir]], bJ[[ir + 1]]}]; 
    ]; 

aJ = Re[aJ]; 
rr = Sqrt[2 aJ]; (*radii of the rings*) 


NHTRingsPlot = 
ParametricPlot[Evaluate[Table[{rr[[i]] Cos[u], rr[[i]] Sin[u]}, {i, 1, Nr}]], {u, 0, 2 Pi}, PlotStyle -> {Blue},(*PlotLabel\[Rule]"Nested Rings",*) AxesLabel -> {"z", "\!\(\*SubscriptBox[\(p\), \(z\)]\)"},LabelStyle -> Directive[{FontSize -> 16}]] 

в Python:

import numpy as np 
import scipy 
import math 
from scipy.integrate import quad 

Nr = 5.  
Jmax = 20. 

CF = math.ceil(scipy.integrate.quad(lambda II : math.exp(-II), 0, Jmax)[0]) 

bJ = np.array([0.,0.,0.,0.,0.,0.]) 
bJ[0] = 0. 
bJ[Nr] = Jmax 

def func(): 
    func = quad(lambda JJ : (math.exp(-JJ)/CF), bb, Jroot) - 1/Nr 
    for i in bJ: 
     i=1 
     if i < Nr: 
      i += 1 
      bb = bJ[i] 
      bJ = np.roots(func, Jroots, bb + CF/(Nr * math.exp(-bb))) 
      bJ.append() 
     return bJ 

Проблема с FUNC :

TypeError: неподдерживаемый тип операнда (ов) для -: 'кортеж' и 'поплавком'

в d с ббом:

NameError: название «бб» не определен

Я только начал изучать Python, и буду признателен, если кто-нибудь из вас может мне помочь с этим.

Answer 1

Существует ряд проблем с тем, что вы написали. Поскольку я не совсем понимаю, что вы пытаетесь выполнить, это не будет ответ, но я укажу на проблемы, которые вижу. Вероятно, вы должны провести некоторое время с The Python Tutorial в документации.

CF = math.ceil(scipy.integrate.quad(lambda II : math.exp(-II), 0, Jmax)[0]) 

Правая сторона задания будет вычисляться 1.0 так CF всегда будет 1.0

def func(): 
    func = ..... 

Вы начинаете определение функции для функции с именем func, то вы сразу же назначить что-то другое имя func. Взгляните на Naming and binding и Defining Functions.

TypeError: unsupported operand type(s) for -: 'tuple' and 'float' 

scipy.integrate.quad возвращает кортеж, вы не можете выполнить операцию по математике (- 1/Nr) на кортеж.

NameError: name 'bb' is not defined 

Вы не назначили ничего bb но так оно не определено. Если вы хотите создать часть частичной функции для использования позже с различными значениями для bb, взгляните на functools.partial. Другой альтернативой было бы определить функцию позже в коде после того, как bb будет программно изменен.

>>> print bb 

Traceback (most recent call last): 
    File "<pyshell#34>", line 1, in <module> 
    print bb 
NameError: name 'bb' is not defined 
>>> # assign 1 to bb 
>>> bb = 1 
>>> print bb 
1 
>>> 

Это очень трудно расшифровать, что вы пытаетесь сделать с вашим for цикла (for i in bJ:), но вы не используете его правильно, for Statements, The for Statement - вот краткий пример:

>>> a = np.array([1.0, 2.0, 3.0, 4.0]) 
>>> a 
array([ 1., 2., 3., 4.]) 
>>> for n in a: 
    print n 


1.0 
2.0 
3.0 
4.0 
>>> 

for i in bJ: присваивает последовательные элементы/элементы от bJ до i на каждой итерации, а затем в наборе циклов for вы назначаете разные значения i, что неверно.Вы делаете это

>>> for c in 'abcde': 
    print 'c is ', c 
    c = 2 
    print 'now c is ', c 

c is a 
now c is 2 
c is b 
now c is 2 
c is c 
now c is 2 
c is d 
now c is 2 
c is e 
now c is 2 
>>> 

One действительно хорошая вещь о Python, что вы можете попробовать вещи в раковине и посмотреть, как они работают или не работают.