在数值积分中求根(积分极限)

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:元组不支持的操作数类型：-和float

使用bb：

NameError:未定义名称'bb‘

我刚刚开始学习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
>>> 

关于Python的一个真正的好处是你可以在shell中尝试一些东西，看看它们是如何工作的，或者不工作的。