在数值积分中求根（积分极限）答案

【问题标题】：Find root (limit of integration) in numerical integration在数值积分中求根（积分极限）
【发布时间】：2015-03-30 22:58:34
【问题描述】：

我正在尝试重写 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: 不支持的操作数类型 -: 'tuple' 和 'float'

和与bb：

NameError: name 'bb' 未定义

我刚开始学习 Python，如果有人能帮助我，我将不胜感激。

【问题讨论】：

对于我们这些可能受到 Mathematica 挑战的人，请在您的 Mathematica 代码中添加一些 cmets 或描述它在做什么，或将我们指向您尝试实现的算法。此外，说 I'm stuck 也不是很有帮助 - 结果有什么问题。

标签： python numpy scipy wolfram-mathematica numerical-integration

【解决方案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 = .....

您开始为名为@987654332@ 的函数定义函数，然后立即为名称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 分配任何东西，所以它没有被定义。如果您想创建一个 partial 函数以供以后使用 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 中进行尝试，看看它们如何工作或不工作。

【讨论】：