Fusion809's profile - activity

I ran:

derivative(acosh(x), x)

and the result was:

1/(sqrt(x + 1)*sqrt(x - 1))

I personally would prefer 1/sqrt(x^2-1) as it's an expanded form and I'm wondering how I can get SageMath to show this. I've tried adding .expand() afterwards and that made no difference to the result.

This specific integral isn't symbolically integratable at all; I originally meant to write the bounds as -infinity to infinity, in that case the solution is the Airy $2\pi Ai(z)$ function. Complex integrals that are analytically integratable (e.g. $\int_{0}^{2\pi} e^{ix} dx$) can be analytically (or symbolically) solved in SageMath, but not numerically solved.

Well, it's a roundabout way of solving it, but sure it's satisfactory. Although, if anyone has a more direct way of numerical integration with complex numbers I'll be happy to accept such an answer.

If you look at my question the conditions are: $y(0)=y(\infty)=0$ (granted, I added that like 30 mins before you published this answer, so maybe you missed it). I didn't really expect desolve to work on this, otherwise my DuckDuckGo search probably would have detected a mention of its use for this purpose (SLEs are a different kettle of fish to usual ODEs in terms of how they are solved, so I would have expected a special mention of them, if desolve could solve them). Nonetheless, I gave it a go, here's the code I used and I placed the error in a comment: https://gist.github.com/fusion809/b3e....

z=1; numerical_integral(exp(i*(x*z+(x^3)/3)), 0, 1, max_points=100) returns the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-a00a879c8f90> in <module>()
----> 1 z=Integer(1); numerical_integral(exp(i*(x*z+(x**Integer(3))/Integer(3))), Integer(0), Integer(1), max_points=Integer(100))

/usr/lib/python2.7/site-packages/sage/calculus/integration.pyx in sage.calculus.integration.numerical_integral (build/cythonized/sage/calculus/integration.c:2943)()
    279                 to_sub = dict(zip(vars[1:], params))
    280                 func = func.subs(to_sub)
--> 281             func = func._fast_float_(str(vars[0]))
    282         except (AttributeError):
    283             pass

/usr/lib/python2.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._fast_float_ (build/cythonized/sage/symbolic/expression.cpp:62032)()
  11928         """
  11929         from sage.symbolic.expression_conversions import fast_float
> 11930         return fast_float(self, *vars)
  11931 
  11932     def _fast_callable_(self, etb):

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in fast_float(ex, *vars)
   1704         1.4142135623730951
   1705     """
-> 1706     return FastFloatConverter(ex, *vars)()
   1707 
   1708 #################

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    224             return self.tuple(ex)
    225         else:
--> 226             return self.composition(ex, operator)
    227 
    228     def get_fake_div(self, ex):

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in composition(self, ex, operator)
   1679         """
   1680         f = operator
-> 1681         g = [self(_) for _ in ex.operands()]
   1682         try:
   1683             return f(*g)

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    216                 div = self.get_fake_div(ex)
    217                 return self.arithmetic(div, div.operator())
--> 218             return self.arithmetic(ex, operator)
    219         elif operator in relation_operators:
    220             return self.relation(ex, operator)

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in arithmetic(self, ex, operator)
   1651             from sage.functions.all import sqrt
   1652             return sqrt(self(operands[0]))
-> 1653         fops = map(self, operands)
   1654         if operator == add_vararg:
   1655             operator = _operator.add

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    215             if getattr(self, 'use_fake_div', False) and (operator is _operator.mul or operator is mul_vararg):
    216                 div = self.get_fake_div(ex)
--> 217                 return self.arithmetic(div, div.operator())
    218             return self.arithmetic(ex, operator)
    219         elif operator in relation_operators:

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in arithmetic(self, ex, operator)
   1651             from sage.functions.all import sqrt
   1652             return sqrt(self(operands[0]))
-> 1653         fops = map(self, operands)
   1654         if operator == add_vararg:
   1655             operator = _operator.add

/usr/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    206         except TypeError as err:
    207             if 'self must be a numeric expression' not in err.args:
--> 208                 raise err
    209 
    210         operator = ex.operator()

TypeError: unable to coerce to a real number

Likewise:

z=1; integral(exp(i*(x*z+(x^3/3))), x, 0, +Infinity)

returns:

integrate(e^(1/3*I*x^3 + I*x), x, 0, +Infinity)

. If I add a .n() after it, to numerically evaluate it, I get much the same error as I showed earlier. So, how can I do numerical integration of complex integrands in SageMath? I know that simple complex integrands (e.g. exp(i*x)) can be analytically integrated by SageMath, but numerical integration does not seem to be an option.

Incorrect, because it does not take the boundary conditions into account, nor does it take into account that lambda is an eigenvalue, not just any old parameter. The solution only satisfies the boundary conditions when lambda is a zero of the Airy Ai(x) function.

Thanks, shall do.

I've conducted DuckDuckGo searches for "Sturm-Liouville sagemath" and received no usable results, and as a result I am here to ask, how can one solve Sturm-Liouville problems with SageMath? I know how to do so numerically with Julia, GNU Octave and Scilab, using matrices, but is there any way to solve at least some of them analytically in SageMath (e.g. $\frac{d^2 y}{dx^2}-xy = \lambda y$ has solutions $y_n=C\cdot \textit{Ai}(x+\lambda_n)$, where $\lambda_n $ are the zeros of the Airy $\textit{Ai}(x)$ function, if $y(0)=y(\infty)=0$)? Failing this, is there a way to numerically solve them in SageMath (without resorting to its interfaces with the languages I just mentioned)?

If it helps, here's the code I use in GNU Octave to solve the aforementioned SLE:

# Standard values are for transformation x=(y-L)/(y+L), y in 0, inf
# *lin are for x=2/L*y-1 (linear transformation)
clear all
format long e;
N=5000;
perc=0.266
Nfrag=round(perc*N);
n=0:N;
# 75 at N=1k; 162 for N=3k; 225 for N=5k;
L=226;
t=pi+n'*pi/N;
x=cos(t);
y=L*(cot(t/2)).^2;
ysub=y(2:N);
ylintrans = L/2*(x+1);
# Chebyshev polys of the first kind
T         = cos(t*n);
# Now for arrays that do not include the endpoints
xsub      = x(2:N);
Tsub      = T(2:N,:);
Usub      = diag(1./sqrt(1-xsub.^2))*sin(acos(xsub)*n);
dTsub     = Usub*diag(n);
dT        = [-(n.^2).*(-1).^(n); dTsub ; (n).^2];
dTL       = diag(2*L./((y+L).^2))*dT;
dTlin     = 2/L*dT;
D1        = dTL/T;
D1lin     = dTlin/T;
D2        = D1^2;
D2lin     = D1lin^2;
H         = D2-diag(y);
Hlin      = D2lin-diag(ylintrans);
H         = H(2:N,2:N);
Hlin      = Hlin(2:N,2:N);
[Y, Lam]  = eig(H);
[Ylin, Lamlin] = eig(Hlin);
Lam       = diag(Lam);
Lamlin    = diag(Lamlin);
[Lam, IX] = sort(Lam, 'ascend');
[Lamlin, IXlin] = sort(Lamlin, 'ascend');
Y         = Y(:,IX);
Ylin      = Ylin(:,IXlin);
aizerochk = abs(airy(0,Lam));
aizerochklin = abs(airy(0,Lamlin));
L
errrms    = sqrt(aizerochk(1:Nfrag)'*aizerochk(1:Nfrag)/(Nfrag+1))
errrmslin    = sqrt(aizerochklin(1:Nfrag)'*aizerochklin(1:Nfrag)/(Nfrag+1))
ain=quad("airynorm",0,100);
figure(1)
plot(ysub(1:3*N/5), Y(1:3*N/5,1));
figure(2)
plot(ysub(1:3*N/5), Y(1:3*N/5,2));
figure(3)
plot(ysub(1:3*N/5), Y(1:3*N/5,Nfrag));

. It uses a Chebyshev spectral method on the extrema grid, and then analyses the results.

I have been attempting to perform a 3D wireframe plot of the solution to the Lorenz equations, which is stored in the variables X, Y and Z. This is what I am presently using (unsuccessfully I might add):

 # Plot the result
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.figure(1)
plt.xlabel('X(t)')
plt.ylabel('Y(t)')
plt.zlabel('Z(t)')
axes3d.plot_wireframe(X,Y,Z)
plt.show()

How do I get this plot to work? I'm guessing you'll probably be able to guess that X, Y and Z are ndarrays, produced from this SageMath script:

x,y,z=var('x,y,z')

# Next we define the parameters
sigma=10
rho=28
beta=8/3

# The Lorenz equations
lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]

# Time and initial conditions
N=250000
tmax=250
h=tmax/N
t=srange(0,tmax+h,h)
ics=[0,1,1]
sol=desolve_odeint(lorenz,ics,t,[x,y,z],rtol=1e-13,atol=1e-14)
X=sol[:,0]
Y=sol[:,1]
Z=sol[:,2]

I would like to implement the square root of 2 power series:

$ \sqrt{2} = -\sum \limits_{n=0}^{\infty}{\frac{(-1)^n (2n-3)!!}{2^n \times n!}} $

(which I obtained from the Maclaurin series of $\sqrt{1+x}$ with $x=1$) in SageMath but I cannot seem to find the double factorial function in the SageMath docs. Is there one? I suppose if a double factorial function is not available I could use this method of finding the double factorial that I found on Wolfram MathWorld:

$\Gamma(n+\frac{1}{2}) = \frac{(2n-1)!!}{2^n}\sqrt{\pi}$

Ya I'm well aware convergence is an issue. This method with N (number of terms) = 10,000 is only accurate to 6 decimal places. I know a far more convergent method of estimating $\sqrt{2}$ is Newton's method. Just wanted to give this a try.