Mark's profile - activity

for an answer to this and related other questions see http://thingwy.blogspot.de/

Further, your suggestion that MM "knows" how to perform this integral analytically and does so is an egregious error. The Fermi-Dirac integral of negative half order (used in your example) must also be done numerically (quick google search). Even the answer given by MM, a number and not an expression, suggests that it is not employing a symbolic solution, but is indeed switching to a numerical method

A) unrelated to sage: as I already said in my previous post, I'd be careful making all-out statements. (Even more so based on 'quick' google seraches). So once again, to falsify your claims:

[1] The Analytical Evaluation of the Half-Order Fermi-Dirac Integrals Jerry A. Selvaggi, and Jerry P. Selvaggi The Open Mathematics Journal, 2012, 5, 1-7

[2] Solutions to the Fermi-Dirac Integrals in Semiconductor Physics Using Polylogarithms M. Ulrich, W. Seng, P. Barnes Journal of Computational Electronics 1: 431–434, 2002

and several other papers. You should also realize, that there is a difference between scientific papers, like the preceding, which have undergone peer review, and some unpublished and unrefereed 'notes', which pop up on the internet - like the one you point to @ nanohub.org

From the preceeding, I guess it is quite clear, where or to whom to apply vocabulary like 'egregious'. Moreover, it is Ref. [2] which tells you that things work down even to y>-1. Just pls. read the papers before posting.

B) related to sage: MM is not '..programmed to recognize Fermi-Dirac integrals and respond accordingly..'. Basically it uses eqn. (5) of Ref. [2], which is a paper free to be read with results free to be used by anyone - also outside of academia. Obviously MM has simply cared to include this result into its code - while sage has not. As simple as that.

with my tail between my legs

C) related to this forum I like this forum because it tends to be free of such statements. Pls. keep it that way.

I'd be the last to put the brake on your sage enthusiasm ... and I can't answer your question on SAGE cloud usage ... so, admittedly my answer is kind of unrelated to this forum ... but if you are really up to writing scientifically competitive numerical Monte Carlo production code including - as your question suggests - parallelism, i.e. most likely code that makes use of MPI, then, to my humble understanding you may want to reconsider, if a CAS like sage is really what you want to use for that - before you get stuck. Surely you can try doing your best at using essentially the C(P)ython subset of sage including libs like mpi4py. But even then - and I have lots of examples for that - it is very very likely, that your MC code will perform poorly, to say the least, as compared to any straightforward C++ implementation.

There are two comments in order here:

A) as you state yourself, the integrands of your expression, e.g.

U1*phi1[0]
y |--> (y - 1)*((y - 1)*(y + 1)*a12*y + (y - 1)*(y^2 + 1)*a13*y + (y - 1)*a11*y)*y

contains lots of symbolic varibles: a12,a13,a11 which are no valid floats, even after substituting a valid float for the integration variable y. That's what the error tells you ;-) . SciPy's quad has to know the numerical values of such parameters - how else could it do a numerical integral? How to supply float parameters to quad is described at http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad E.g.:

f(y,a12,a13,a11)=U1*phi1[0]
quad(f,0.0,1.0,args=(1,1,1,))
(0.1261904761904762, 1.4009957215507927e-15)

B) having said that, and not knowing what your actual problem is, I however presume that the preceeding is not what you actually want, and that rather you want to set up some kind of linear equation (eq1 ?) with symbolic variables and numerically evaluated 'constants'. If that is so, I would disect the numerical integration of each of those constant off from the symbolic part of the code.

@whomeveritmayconcern: I'd dare to disagree with the notion that a matrix is not sparse or dense by itself. I'd rather go with the majority of mathematical definitions on that, which can also be found publicly, e.g. http://en.wikipedia.org/wiki/Sparse_matrix, or http://www.phy.ornl.gov/csep/la/node14.html, or many others. Also, the 'difference' between sparse and dense matrices ist not how you store them ... I believe matrices exist even without computer RAM ;-) . Rather, once a matrix has been identified to be sparse, then you may want to consider storing it into various existing sparse matrix formats and may want to use various sparse matrix algorithms ... yet, even if you would stay away from doing so, the matrix would still remain sparse by itself.

For a matrix $M$ of rank $N\times N$, the number of non-zero elements should scale $\sim N$ or, at most, $\sim N \log(N)$ for $M$ to be sparse, otherwise most sparse matrix algorithms are of rather little help. (You can surely refine this statement for $N\times M$ matrices with $N\neq M$ ...)

In turn, I would not call 'your' matrix sparse.

First, when asking questions, it might speed up a potential answer if you would care to give enough detail about your problem ... so, I only guess that this is what you want $$F(y,\mu) = \int_0^\infty \frac{x^{y}}{e^{x+\mu}+1} dx$$

Second, I'd be very careful about www-ing statements like

"...which can only be done numerically..."

since, to my recollection $F(y,\mu)$ is just another way of writing the polylog for all $y>-1$, i.e down to those values of $y$ where the singularity of the integrand remains harmless. Therefore, at least I do not understand your claim, since in 1D the DOS in a semiconductor is $\propto x^{-1/2}$ (which you also did not specify that in your question). So it can be done analytically ...

Third, if still you really want to do this numerically in sage ... just look here http://www.sagemath.org/doc/reference/calculus/sage/gsl/integration.html and e.g. do this

def F(y,mu):
    ul = 100+mu # increase the 100 if you really need
    return numerical_integral(x^y/(e^(x+mu)+1),0,ul)

F(-1/2,0)
(1.072154920510084, 8.286087609477639e-07)

Finally, and only to my very humble understanding: sage is a chain of existing tools, i.e. it is kind of as strong as it's weakest link. Since in plain sage things like Maxima are the link for many of those calculations which 90%-99% of non-mathematician might call 'symbolic', plain sage frequently is of little help if compared to MM or Maple once you come to such 'symbolic' problems.

Applied to your case this means:

sage: var('y,mu')
(y, mu)
sage: assume(y+1>0)
sage: assume(y,'integer')
sage: integral(x^y/(e^(x+mu)+1),x,0,oo)
integrate(x^y/(e^(mu + x) + 1), x, 0, +Infinity)

so, nothing ... or if you try special values

sage: integral(x^(-1/2)/(e^(x+mu)+1),x,0,oo)
integrate(1/((e^(mu + x) + 1)*sqrt(x)), x, 0, +Infinity)

sage: integral(x^1/(e^(x+mu)+1),x,0,oo)
polylog(2, -e^mu) + limit(1/2*x^2 - x*log(e^(mu + x) + 1) - polylog(2, -e^(mu + x)), x, +Infinity, minus)

again, not too helpful. Maybe there are 'tricks' to get further but I don't know.

If however you use MM then, out of the box:

sage: mathematica('Integrate[x^a/ (Exp[x + b] + 1), {x, 0, Infinity}, Assumptions -> {Re[a] > -1}]')
-(a*Gamma[a]*PolyLog[1 + a, -E^(-b)])

so MM 'knows' that your integral can be done analytically and provides the answer, and

sage: mathematica('-(a*Gamma[a]*PolyLog[1 + a, -E^(-b)])/. {a -> -1/2, b -> 0} // N')
1.072154929940191

which is the numerical value we got previously.

In conclusion: To answer if that implies you should better be using MM, there are more powerful experts around here.

If you just want to have a numerical solution to a small set of (nonlinear) equations, why not simply 'fsolve' them?

from scipy.optimize import fsolve

def equations(p):
    a, b, c = p
    return (c + b*abs(a)**2 - 10, c + b*abs(a)**4 - 6, c + b*abs(a)**5 - 5)

a,b,c = fsolve(equations, (1, 1, 1))

print a,b,c
print equations((a,b,c))
0.640388203202 16.5345651552 3.2192235936
(-2.5757174171303632e-13, -2.3625545964023331e-13,-1.9895196601282805e-13)

The 'abs()' is there to provide at least some additional constraining on a,b,c

I don't know if this is elegant - but it is straightforward.

It is well known that for $n\in\mathbb{N}$ and $n>0$ (an maybe even for more than these restrictions): $$I_n = \int_0^\infty\frac{x^n}{e^x-1}dx = \zeta(n+1)n!$$ which, analytically can be shown easily by expanding $1/(1-e^{-x})$ into a geometric series, which leads to trivial integrals, and by using $\zeta(n+1)=\sum_{l=1}^\infty l^{-(n+1)}$. So, eg.: $$I_1 = \pi^2/6$$ $$I_2 = 2\zeta(3)$$ a.s.o...

Now, if I try even the simplest case with sage, I get this 'nifty' little results

sage: integrate(x/(exp(x)-1),x,0,oo)

-1/6*pi^2 + limit(-1/2*x^2 + x*log(-e^x + 1) + polylog(2, e^x), x,
+Infinity, minus)

Is there any trick to simplify this down to the final result, or is this about as far as I can get with sage alone?

PS.: it is probably needless to say that (once again ... :( ...)

In[1]:= Integrate[x/(Exp[x] - 1), {x, 0, Infinity}]
Out[1]:= Pi^2/6

If I have two plain formal functions, like

sage: var('q k')
sage: function('u')
sage: function('v')

it seems I can always get more 'involved' formal functions, like the following h(k,q), by symbolic operations

sage: h(k,q)=u(k)*v(q)
sage: h(17,sin(k))
u(17)*v(sin(k))       # OK.

My naive hope was, that this would straightforwardly apply also to matrix construction from formal functions, like so

sage: B(k,q) = matrix([[u(k),v(q)],[v(q),u(k)]])
# but
sage: B(17,sin(q))
[u(k) v(q)]
[v(q) u(k)]       # :(

So, B(k,q) does not substitute the variables. After some attempts, it seems that this is a valid way to construct a callable symbolic matrix from formal functions

sage: B=matrix(CallableSymbolicExpressionRing((k,q)),[[u(k),v(q)],[v(q),u(k)]])
# at least
sage: B(x^2,cos(q))
[   u(x^2) v(cos(q))]
[v(cos(q))    u(x^2)]        # OK.

My 1st question is, if there is not a 'less cryptic' way to construct B from u and v?

My 2nd question concerns doing matrix arithmetic with B. If I try to do something similar to that for the plain formal function arithmetic for h(k,q) with my formal matrix function B it 'does not work'. Eg.:

sage: H(k,q)=B(k,q)*B(q,q)
sage: H(sin(17),pi)
[   u(k)*u(q) + v(q)^2 u(k)*v(q) + u(q)*v(q)]
[u(k)*v(q) + u(q)*v(q)    u(k)*u(q) + v(q)^2]  # :(

I.e., again, upon matrix multiplication, the object B(k,q)*B(q,q) seems to not substitute the arguments into the formal functions anymore. However, if I do

sage: H=matrix(CallableSymbolicExpressionRing((k,q)),B(k,q)*B(q,q))

Then

sage: H(sin(17),pi)
[    u(pi)*u(sin(17)) + v(pi)^2 u(pi)*v(pi) + u(sin(17))*v(pi)]
[u(pi)*v(pi) + u(sin(17))*v(pi)     u(pi)*u(sin(17)) + v(pi)^2]    # OK.

But it seems awkward if one would have to do matrix arithmetic like that. So, how does one do proper matrix operations with matrices made from formal functions?

Ok., this makes (some) sense. Still, would you mind educating me on if this is just syntactic candy, or rather something absolutely standard that I should find in my high school linear algebra textbooks.

given a vector in sage

sage: a,b=var('a,b')
sage: sv=vector(SR,[1,a,b^2])

then why is this a meaningful operation

sage: sv/sv
1

To my (utterly humble) understanding, I always thought, that if vector multiplication '*' is defined via the scalar product - as it seems to be in sage

sage: sv*sv
a^2 + b^4 + 1

then, division cannot be defined meaningfully?

[Just as a sideremark along that line: I have no problem with numpy's element wise array-arithmetic

In [1]: v=array([1.,2.,3.]) # vector
In [2]: e=array([1.,1.,1.]) # unity is a vector
In [3]: e*v == v # multiplication by unity
Out[4]: array([ True,  True,  True], dtype=bool)
In [5]: vi=v**(-1) # inverse is a vector
In [6]: e/v == vi # unity/vector == inverse
Out[7]: array([ True,  True,  True], dtype=bool)
In [8]: e == v*vi # vector * inverse == unity
Out[9]: array([ True,  True,  True], dtype=bool)

In my layman's world, this type of division makes perfect sense. (Moreover I'm curious why sage seems to not adopt this pythonic way of array arithmetic)]

(Maybe) that helps. Let me ask one for clarification. If I get you right, kind of your point is, that callable symbolic expression should allow for meaningful usage with symbolic operations like multiplication, addition a.s.o. But while a scalar product, as f*f may make sense (even though I wasn't even looking for that) the division f(x)/f(x) is also a symbolic operation, which does not make sense for me in case of vectors - although sage simply return f(x)/f(x) = 1 with my definition (1)

Can you provide me with an example of what you mean by 'unsafe' or 'not meaningful' if the callable symbolic expression would actually return what I have asked for, namely a list?

I was (wrongfully) expecting the return type of all of the following function calls to be the same, namely a list. Why does the callable symbolic expression (1) return a vector, and how to rewrite (1) such that it does return a list?

# (1)
sage: f(x) = [x,x]
# (2)
sage: g = lambda x: [x,x]
# (3)
sage: def h(x):
....:     return [x,x]

sage: type(f(x))
<class 'sage.modules.vector_symbolic_dense.Vector_symbolic_dense'>
sage: type(g(x))
<type 'list'>
sage: type(h(x))
<type 'list'>

I would not agree with you that these solutions are 'wrong'. They are just presented in a very user unfriendly way. The 1st one, i.e.:

[x == 1, -2 > 0]

definitely is a solution - it just specifies a range of measure zero, since

sage: -2>0
False

so, it is not 'wrong' as you state - it is just 'totally useless' ;) The other solution is also perfectly ok since

sage: bool(-(t - 2)^2 + 2*(t - 2)*t - t^2 + 6)
True

so, the solution you want is

x == (t - 2)/t

which I think is ok ..... however I would definetely agree with you that

In[1]:= f[x_, t_] = Exp[x t] (x - 1)^2;
In[2]:= Solve[{D[f[x, t], x] == 0 , D[f[x, t], {x, 2}] < 0, {x, t}]
Out[2]= {{t -> ConditionalExpression[-(2/(-1 + x)), (t | x) \[Element] Reals]}}

from a well known competitor is far more user friendly

I'd like to look at the following: the set of formal functions with 2 variables f(x,y) and the real numbers a,b,c..., including an addition and a non-commutative multiplication, such that things like

expand((a+f(x,y))*(b+c*f(u,v))) = a*b + a*c*f(u,v) + b*f(x,y) + c*f(x,y)*f(u,v)

are possible (and vice versa), and with the multiplication of the functions

f(x,y)*f(u,v) != f(u,v)*f(x,y)

being non-commutative, however with the multiplication of the functions by the real scalars

a*f(x,y) == f(x,y)*a

still commutative. Can I construct something like that with sage?

Yes, I was expecting f(x)==f(x) to automatically simplify to True. (Mathematica will actually do so.) The bool()-typing workaround you suggest is kind of ok for me, yet, as I read your answer, it seems really strange that x-x==0 does not result in True but rather in 0==0 ... and moreover, if you type 0==0 at the sage prompt you will obviously get True ...

Why does the following not result in 'True'?

var('x')
f=function('f',nargs=1)
f(x)==f(x)

or, speaking differently, can I define formal functions such that the preceding will result in 'True'?

1. 'simple' is a matter of definition (or taste)
2. even in cases where it is rather obvious what the 'simple' answer would be, sage will frequently not live up to that ... see this one

anyway, I guess what you want is

u.expand()

ehem, guys, for a newbs 1st-in-a-life-time three-liner typed into sage this is a little over my head. If this nobrainer integral results in a bug already, is your point that I should use integrate() only in conjunction with simplify_full()? Or is this a 'rare' case? Or am I doing something wrong?

why is this integral performed:

forget()
var('n')
assume(n>0)
integrate(1/sqrt(1+x^2*n),x,1,2)

i.e. I get:

-arcsinh(sqrt(n))/sqrt(n) + arcsinh(2*sqrt(n))/sqrt(n)

while this is not

forget()
var('n')
assume(n>0)
integrate(1/sqrt(1+x^2/n),x,1,2)