ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 29 Jan 2024 19:11:05 +0100How to take inverse of matrix with complex entrieshttps://ask.sagemath.org/question/75698/how-to-take-inverse-of-matrix-with-complex-entries/I am trying to find the similarity matrix between two 18 by 18 matrices `K_1` and `K_2` which are cospectral. The `is_similar()` function takes too long to run. We can compute the similarity matrix by finding the transition matrices `S` and Q which diagonalize both `K_1` and `K_2`. In other words, `S^(-1)*K_1*S = H = Q^(-1)*K_2*Q`, where `H` is a diagonal matrix. We then know that `Q*S^(-1)*K_1*S*Q^(-1) = K_2`. My problem now though is that we can't compute `Q^(-1)` or `S^(-1)` since both `Q` and `S` have nasty complex numbers and the inverse method in Sage doesn't run fast enough. Any thoughts on how to quickly compute these inverses?
fyi, I did try using the left and right eigenmatrices to find the similarity matrix, but the eigenvectors didn't line up.luke_mattgreenMon, 29 Jan 2024 19:11:05 +0100https://ask.sagemath.org/question/75698/Real and imaginary part of an expression with assumptions on variableshttps://ask.sagemath.org/question/73517/real-and-imaginary-part-of-an-expression-with-assumptions-on-variables/ Hello,
I have this expression:
Exp[-(x+I*y)*sinh(t/(1-t)) - nu * t/(1-t)] * 1/(1-t)
All variables are real except nu which is a complex number. I want the real part and the imaginary part of this expression. I tried with Wolfram Alpha and either it assumes all variables are complex, or, of I specify some Assumptions, it does not return the result. Is it possible with Sage ? I'm very new to Sage. I'm new to GIAC too and I didn't found how to force a variable to be real.
stlaWed, 20 Sep 2023 17:45:17 +0200https://ask.sagemath.org/question/73517/Getting exact answers when using eigenvectors_right()https://ask.sagemath.org/question/66147/getting-exact-answers-when-using-eigenvectors_right/ I am trying to compute complex eigenvectors for a matrix that has two distinct complex eigenvalues. Is there a way to get the exact coefficients?
For example:
A=matrix(QQ,2,2,[3, -13, 5,1])
A.eigenvectors_right()
Gives me:
[(2 + 8*I,
[
(1, 0.07692307692307692? - 0.6153846153846154?*I)
],
1),
(2 - 8*I,
[
(1, 0.07692307692307692? + 0.6153846153846154?*I)
],
1)]
Is there any way to make it give me integer or rational Re and Im for the eigenvectors?
Thank you!
TatiMon, 30 Jan 2023 19:09:59 +0100https://ask.sagemath.org/question/66147/How can I calculate limit(cos(x+I*y),y=infinity) ?https://ask.sagemath.org/question/65551/how-can-i-calculate-limitcosxiyyinfinity/ How can I calculate limit(cos(x+I*y),y=infinity) ?hichemeMon, 26 Dec 2022 12:17:32 +0100https://ask.sagemath.org/question/65551/complex limitshttps://ask.sagemath.org/question/65172/complex-limits/Given that $z$ is a complex number of the form $z = a +b i$, where $a, \, b \in \mathbb{R}$, and $\overline{z}$ is its conjugate, what is
$$\lim_{z \to 0} \frac{\left( \overline{z} \right)^2}{ z^2} =\text{?}$$
When I ask Sage to do this limit, it says it's 1. I am almost certain that's not right. Here's the Sage code:
`sage: z=var('z')
sage: assume(z,'complex')
sage: limit((conjugate(z))^2/z^2,z=0)
`
My analysis leads to DNE, so I wonder if I am misusing Sage.b@nnon.usSat, 03 Dec 2022 00:58:04 +0100https://ask.sagemath.org/question/65172/Solving a system of linear equations with complex numbers yields false solutionhttps://ask.sagemath.org/question/64344/solving-a-system-of-linear-equations-with-complex-numbers-yields-false-solution/I have equations for the analysis of an electric circuit. They contain at several places the term w*I (I is the imaginary unit). I get a solution for the currents I1, I2, I3, I4. Then I substitute the solution into the equations, but I see that two of the four equations are not fulfilled:
var('L1 L2 L3 L4 C1 C2 C3 C4 M23 I1 I2 I3 I4 U w p RL')
lsg = solve([I1*(w*I*L1 + 1/(w*I*C1)) + I2*(-1/(w*I*C1) ) == U, \
I2*(w*I*L2 + 1/(w*I*C1) + 1/(w*I*C2)) + I3*(w*I*M23 - 1/(w*I*C2)) + I1*(- 1/(w*I*C1)) == 0, \
I3*(1/(w*I*C3) + w*I*L3 + 1/(w*I*C2)) + I2*(w*I*M23 - 1/(w*I*C2)) - I4/(w*I*C3) == 0, \
I4/(w*I*C3) + I4*RL - I3/(w*I*C3) == 0], \
[I1, I2, I3, I4])
param = [w==1, U==1, C1==1, C2==1, C3==1, C4==1, L1==1, L2==1, L3==1, L4==1, RL==1, M23==0]
I1 = I1.subs(lsg).subs(param)
I2 = I2.subs(lsg).subs(param)
I3 = I3.subs(lsg).subs(param)
I4 = I4.subs(lsg).subs(param)
eqn = [I1*(w*I*L1 + 1/(w*I*C1)) + I2*(-1/(w*I*C1) ) == U, \
I2*(w*I*L2 + 1/(w*I*C1) + 1/(w*I*C2)) + I3*(w*I*M23 - 1/(w*I*C2)) + I1*(- 1/(w*I*C1)) == 0, \
I3*(1/(w*I*C3) + w*I*L3 + 1/(w*I*C2)) + I2*(w*I*M23 - 1/(w*I*C2)) - I4/(w*I*C3) == 0, \
I4/(w*I*C3) + I4*RL - I3/(w*I*C3) == 0]
print("I1=", I1)
print("I2=", I2)
print("I3=", I3)
print("I4=", I4)
[eq.subs(param) for eq in eqn]
Output (unexpected):
I1= 1
I2= -I
I3= -2
I4= I - 1
[1 == 1, (-I - 1) == 0, I == 0, 0 == 0]
I have copied and pasted the equations from within the solve() command into the eqn = ... statement. So they are guaranteed to be equal.
Because I have set M23=0 in the parameters, I could remove the terms with M23 in the equations. Then I get the correct solution. I don't understand why not when M23 is present.
Another way to get the correct solution is replacing w*I with p in the equations and setting p=I in the parameters:
var('L1 L2 L3 L4 C1 C2 C3 C4 M23 I1 I2 I3 I4 U w p RL')
lsg = solve([I1*(p*L1 + 1/(p*C1)) + I2*(-1/(p*C1) ) == U, \
I2*(p*L2 + 1/(p*C1) + 1/(p*C2)) + I3*(p*M23 - 1/(p*C2)) + I1*(- 1/(p*C1)) == 0, \
I3*(1/(p*C3) + p*L3 + 1/(p*C2)) + I2*(p*M23 - 1/(p*C2)) - I4/(p*C3) == 0, \
I4/(p*C3) + I4*RL - I3/(p*C3) == 0], \
[I1, I2, I3, I4])
param = [p==I, U==1, C1==1, C2==1, C3==1, C4==1, L1==1, L2==1, L3==1, L4==1, RL==1, M23==0]
I1 = I1.subs(lsg).subs(param)
I2 = I2.subs(lsg).subs(param)
I3 = I3.subs(lsg).subs(param)
I4 = I4.subs(lsg).subs(param)
eqn = [I1*(p*L1 + 1/(p*C1)) + I2*(-1/(p*C1) ) == U, \
I2*(p*L2 + 1/(p*C1) + 1/(p*C2)) + I3*(p*M23 - 1/(p*C2)) + I1*(- 1/(p*C1)) == 0, \
I3*(1/(p*C3) + p*L3 + 1/(p*C2)) + I2*(p*M23 - 1/(p*C2)) - I4/(p*C3) == 0, \
I4/(p*C3) + I4*RL - I3/(p*C3) == 0]
print("I1=", I1)
print("I2=", I2)
print("I3=", I3)
print("I4=", I4)
[eq.subs(param) for eq in eqn]
Output as expected:
I1= 1
I2= -I
I3= -I - 1
I4= -1
[1 == 1, 0 == 0, 0 == 0, 0 == 0]
Does anyone know why it yields false results in the first case (taking w*I instead of p) ?
Albert_ZweisteinFri, 07 Oct 2022 11:25:01 +0200https://ask.sagemath.org/question/64344/complex exponential/trigonometrichttps://ask.sagemath.org/question/7672/complex-exponentialtrigonometric/Is there any command equivalent to the Mathematica's "TrigToExp" and "ExpToTrig"?
Where should I look for general support on complex exponential/trigonometric functions?yamanSat, 25 Sep 2010 21:36:07 +0200https://ask.sagemath.org/question/7672/unique elements in complex numbers listhttps://ask.sagemath.org/question/57315/unique-elements-in-complex-numbers-list/I have a list of complex numbers, but even after converting it to a set, there are a lot of repetitions.
here is a minimal example (in my code i have 10^6 polynomials, here only 2)
from sage.rings.polynomial.complex_roots import complex_roots
x = polygen(ZZ)
r1 = complex_roots( x^3+1)
r2 = complex_roots( x^6+2*x^3+1)
s1 = set([r[0] for r in r1])
s2 = set([r[0] for r in r2])
s1 | s2
which results in
{-1,
0.500000000000000? - 0.866025403784439?*I,
0.500000000000000? - 0.866025403784439?*I,
0.500000000000000? + 0.866025403784439?*I,
0.500000000000000? + 0.866025403784439?*I}
as you can see, the roots must be only 3.
Is there an instruction to reduce a set of complex numbers differing less than a delta?alvaroFri, 28 May 2021 09:55:04 +0200https://ask.sagemath.org/question/57315/Issue with inversion of complex symbolic array initialized with numpyhttps://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/Hello!
I am trying to invert a symbolic matrix that is initialized with a combination of numpy arrays and I get an error. Below i present a simple code that gives the error.
import numpy as np
a = np.zeros((5,5) , dtype = 'complex')
np.fill_diagonal(a,1)
b = var('x')*a
c = matrix(b)
c.inverse()
Gives the error:
> ECL says: THROW: The catch
> MACSYMA-QUIT is undefined.
What I have noticed up until now is that the issue stems from the data type. When I try casting the numpy array to float before turning it into a matrix, it works. However, my actual code makes use of complex coefficients. I believe it might have to do with how the imaginary part is represented in SageMath in comparison to numpy. ( `j` vs `I` )
When I manually create the symbolic array with `I` for the imaginary part, the `.inverse()` has no issue.yorgos_sotWed, 28 Apr 2021 13:31:54 +0200https://ask.sagemath.org/question/56850/simplifying a simplicial sethttps://ask.sagemath.org/question/56477/simplifying-a-simplicial-set/The starting point to this question is a simplicial set with a large number of simplices. Sage is still computing exactly how many, but in dimension 4 there are more than 100,000 simplices. Is there a command similar to X.shrink_simplicial_complex(K), but for which the input is already a simplicial set? I am fairly certain this simplicial set description of a cell complex can be described much, much more efficiently. I would like to be able to compute homology and cohomology, but my guess is that this won't be feasible unless I can first simplify.IngridSat, 03 Apr 2021 10:21:48 +0200https://ask.sagemath.org/question/56477/Issue plotting function that is complex under certain domainhttps://ask.sagemath.org/question/34873/issue-plotting-function-that-is-complex-under-certain-domain/I have been able to recreate my issue as the following simple problem.
Say I want to take a simple cubic expression like the following and solve it for x instead of f. This will no longer be a function since there are multiple answers for the same value of f (namely 3 answers). But I could simply take the 3 expressions and overlay their graphs on the same plot to get the graph I want.
var('x,y')
solution = solve( y==x\*(x-1)\*(x+1), x )
show( solution[0].rhs() )
plot( solution[0].rhs(), (y,0,1.5) )
I have been able to get the results I'd expect using Wolfram Alpha.
Now I believe that source of my issue is coming from the fact that over a certain domain the result is complex. I have been playing around with trying to learn about symbolic rings, but I haven't been able to make heads or tails out of how to use them properly, since every time I do I get a new error. Ideally I'd like to just be able to be able to plot only where the function is only real.Cpt_BringdownTue, 20 Sep 2016 22:02:22 +0200https://ask.sagemath.org/question/34873/Implicit plot with complex functionhttps://ask.sagemath.org/question/54765/implicit-plot-with-complex-function/I have a complex function $f:\mathbb{C}\to\mathbb{C}$ and
want to draw the locus where $f$ is in the interval $[0,1]$.
I made the following (where $f$ is my complex function):
x = var('x')
assume(x, 'real')
y = var('y')
assume(y, 'real')
F = f.subs(z == x + I*y)
P = lambda x, y: F.subs(x=x, y=y).real()
Q = lambda x, y: F.subs(x=x, y=y).imag()
region_plot([P(x, y) <= 1, P(x, y) >= 0, Q(x, y) == 0], (x, 0, 5), (y, -1, 1))
That works well for small functions (low degree)
but for big functions it's too long.
For example with
$$ f(z)=-\frac{{\left(z^{4} - 6 z^{3} + 12 z^{2} - 8 \, z\right)} {\left(z - 1\right)}^{3} {\left(z - 3\right)}}{{\left(2z - 3\right)} {\left(z - 2\right)}^{3} z} $$
there is no problem (a few second), but with
$$ f(z)=-\frac{{\left(z^{8} - 16 z^{7} + 108 z^{6} - 400 z^{5} + 886 z^{4} - 1200 z^{3} + 972 z^{2} - 432 z + 81\right)} {\left(z - 2\right)}^{6} {\left(z - 4\right)} z}{{\left(6 z^{4} - 48 z^{3} + 140 z^{2} - 176 z + 81\right)} {\left(z - 1\right)}^{4} {\left(z - 3\right)}^{4}} $$
that's too long.
I think there should be something better. I've seen for example
the function `complex_plot` works well even for big functions
and seems to do something more complex.
Any idea?Gabriel SoranzoSat, 19 Dec 2020 23:44:26 +0100https://ask.sagemath.org/question/54765/Very different matrix inverse result in complex ring CC and in real ring RR while entry only contain real elementshttps://ask.sagemath.org/question/54611/very-different-matrix-inverse-result-in-complex-ring-cc-and-in-real-ring-rr-while-entry-only-contain-real-elements/I would use an dramatized example here but it was encountered during a project as well.
Suppose two matrix consisted of the same entry of only real elements. However, one was in the real domain RR, and one was in the complex domain CC.
The inverse of the complex matrix and the real matrix different significantly.
Example:
test_S=matrix(RR,5)
for ix in range(0,5):
for iy in range(0,5):
test_S[ix,iy]=sin(ix+iy)*sin(iy+1);
and
test_S_CC=matrix(CC,5)
for ix in range(0,5):
for iy in range(0,5):
test_S_CC[ix,iy]=sin(ix+iy)*sin(iy+1);
which returned
test_S
[ 0.000000000000000 0.765147401234293 0.128320060202457 -0.106799974237582 0.725716283876408]
[ 0.708073418273571 0.826821810431806 0.0199148566748170 0.572750016904307 0.919535764538226]
[ 0.765147401234293 0.128320060202457 -0.106799974237582 0.725716283876408 0.267938303940044]
[ 0.118748392158235 -0.688158561598754 -0.135323401369264 0.211462346264655 -0.630000397639817]
[ -0.636827341031836 -0.871947375471875 -0.0394311173578842 -0.497209097294248 -0.948719639025321]
test_S_CC
[ 0.000000000000000 0.765147401234293 0.128320060202457 -0.106799974237582 0.725716283876408]
[ 0.708073418273571 0.826821810431806 0.0199148566748170 0.572750016904307 0.919535764538226]
[ 0.765147401234293 0.128320060202457 -0.106799974237582 0.725716283876408 0.267938303940044]
[ 0.118748392158235 -0.688158561598754 -0.135323401369264 0.211462346264655 -0.630000397639817]
[ -0.636827341031836 -0.871947375471875 -0.0394311173578842 -0.497209097294248 -0.948719639025321]
However, their inverse
(test_S)^-1
[ 2.87043424444896e15 -9.36229301025494e15 1.13569572695379e16 -8.12065238284797e15 1.72141024999095e15]
[-1.10320559226948e15 3.19497062929143e15 5.27553667803348e15 -7.22878751888619e15 8.54302211670640e15]
[ 1.94664421070014e16 -3.02112860544102e15 0.000000000000000 1.80143985094820e16 0.000000000000000]
[ 8.29739198582516e14 9.39732200089164e15 -1.02939420054183e16 1.02939420054183e16 0.000000000000000]
[-2.15677123538261e15 -1.45141888820365e15 -7.07708512872506e15 5.95118522188244e15 -9.00719925474099e15]
(test_S_CC)^-1
[-1.00798918109020e16 1.87532050063097e15 1.07471168476593e16 -2.09935440833922e16 1.10831351979113e16]
[ 2.64758058163658e15 7.44366790239271e15 -8.04492396917722e15 9.82621377301499e15 4.42753908914828e14]
[ 2.35701470124050e16 -4.71646526274252e15 1.02939420054183e16 1.08658276723860e16 9.15017067148291e15]
[ 1.53627372015977e16 -1.32579194813951e15 -1.02939420054183e16 2.51629693465780e16 -9.15017067148291e15]
[-4.69821640166124e15 -7.20926442861200e15 5.14697100270914e15 -8.57828500451523e15 -3.43131400180609e15]
where
(test_S)^(-1)-test_S_CC^(-1)
[ 1.29503260553510e16 -1.12376135108859e16 6.09840421878568e14 1.28728917005442e16 -9.36172494792037e15]
[-3.75078617390606e15 -4.24869727310128e15 1.33204606472107e16 -1.70550012919012e16 8.10026820779157e15]
[-4.10370490540364e15 1.69533665730150e15 -1.02939420054183e16 7.14857083709603e15 -9.15017067148291e15]
[-1.45329980030152e16 1.07231139490312e16 0.000000000000000 -1.48690273411597e16 9.15017067148291e15]
[ 2.54144516627863e15 5.75784554040835e15 -1.22240561314342e16 1.45294702263977e16 -5.57588525293490e15]
Is this a bug? How to fix the issue?
ShoutOutAndCalculateThu, 10 Dec 2020 01:11:16 +0100https://ask.sagemath.org/question/54611/Argand diagram plothttps://ask.sagemath.org/question/38587/argand-diagram-plot/HI
T=[-sqrt(1/2*sqrt(13) + 1) - I*sqrt(1/2*sqrt(13) - 1),
-sqrt(1/2*sqrt(13) + 1) + I*sqrt(1/2*sqrt(13) - 1),
-sqrt(1/2*sqrt(13) + 1) +I*sqrt(1 + sqrt(13)/2),
-sqrt(1/2*sqrt(13) + 1) -I*sqrt(1 + sqrt(13)/2) ]
what is the best , fastest, way to plot Argand diagram of T ?ortolljSun, 20 Aug 2017 11:40:07 +0200https://ask.sagemath.org/question/38587/Error during integration: not of type (UNSIGNED-BYTE 15)https://ask.sagemath.org/question/52913/error-during-integration-not-of-type-unsigned-byte-15/Hi everyone, I am new to Sage and coding, and I am trying to solve a tricky integral from equation 13 in this physics paper: www.nature.com/articles/srep11876
Here is the integrand in code form (I would upload a screenshot of the full equation, but I do not have enough karma points):
sage: var ('alpha beta E_o a d_o z_o n u')
(alpha, beta, E_o, a, d_o, z_o, n, u)
sage: j = CC.0
sage: integrand = e**(-j * n * u) / (((a + d_o + 0.5 * z_o * (((16 * pi) / (alpha * beta)) ** (-3))) + 0.5 * z_o * (((16 *
....: pi) / (alpha * beta)) ** (-3)) * sin(u)) ** 3) - 1
The integral is part of a Fourier coefficient, so it is a definite integral with limits -pi and pi. I want to integrate with respect to u and then solve the full equation for beta. All of the other variables are known.
However, when I try to integrate, I get an error:
sage: integrate(integrand, u, -pi, pi)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-32-84f9192f9549> in <module>()
----> 1 integrate(integrand, u, -pi, pi)
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/misc/functional.py in integral(x, *args, **kwds)
751 """
752 if hasattr(x, 'integral'):
--> 753 return x.integral(*args, **kwds)
754 else:
755 from sage.symbolic.ring import SR
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.integral (build/cythonized/sage/symbolic/expression.cpp:64474)()
12434 R = ring.SR
12435 return R(integral(f, v, a, b, **kwds))
> 12436 return integral(self, *args, **kwds)
12437
12438 integrate = integral
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/integration/integral.py in integrate(expression, v, a, b, algorithm, hold)
929 return indefinite_integral(expression, v, hold=hold)
930 else:
--> 931 return definite_integral(expression, v, a, b, hold=hold)
932
933
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:12262)()
1136 res = self._evalf_try_(*args)
1137 if res is None:
-> 1138 res = super(BuiltinFunction, self).__call__(
1139 *args, coerce=coerce, hold=hold)
1140
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6938)()
595 for i from 0 <= i < len(args):
596 vec.push_back((<Expression>args[i])._gobj)
--> 597 res = g_function_evalv(self._serial, vec, hold)
598 elif self._nargs == 1:
599 res = g_function_eval1(self._serial,
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction._evalf_or_eval_ (build/cythonized/sage/symbolic/function.cpp:13412)()
1224 res = self._evalf_try_(*args)
1225 if res is None:
-> 1226 return self._eval0_(*args)
1227 else:
1228 return res
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/integration/integral.py in _eval_(self, f, x, a, b)
203 for integrator in self.integrators:
204 try:
--> 205 A = integrator(*args)
206 except (NotImplementedError, TypeError, AttributeError):
207 pass
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/integration/external.py in maxima_integrator(expression, v, a, b)
44 result = maxima.sr_integral(expression,v)
45 else:
---> 46 result = maxima.sr_integral(expression, v, a, b)
47 return result._sage_()
48
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/interfaces/maxima_lib.py in sr_integral(self, *args)
788 """
789 try:
--> 790 return max_to_sr(maxima_eval(([max_integrate],[sr_to_max(SR(a)) for a in args])))
791 except RuntimeError as error:
792 s = str(error)
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.EclObject.__call__ (build/cythonized/sage/libs/ecl.c:7794)()
803 """
804 lispargs = EclObject(list(args))
--> 805 return ecl_wrap(ecl_safe_apply(self.obj,(<EclObject>lispargs).obj))
806
807 def __richcmp__(left, right, int op):
/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.ecl_safe_apply (build/cythonized/sage/libs/ecl.c:5456)()
375 if ecl_nvalues > 1:
376 s = si_coerce_to_base_string(ecl_values(1))
--> 377 raise RuntimeError("ECL says: {}".format(
378 char_to_str(ecl_base_string_pointer_safe(s))))
379 else:
RuntimeError: ECL says: 1.8189894035458565e-12 is not of type (UNSIGNED-BYTE 15).
I have no idea where 1.8189894035458565e-12 came from, and I get the same error after assigning values to all of the variables except beta and u.
Has anyone encountered this error before? Is there a better way I should evaluate this?lvb884Sat, 08 Aug 2020 06:03:34 +0200https://ask.sagemath.org/question/52913/Difference between ComplexField and MPComplexField?https://ask.sagemath.org/question/52478/difference-between-complexfield-and-mpcomplexfield/Hi everyone. I'm a long-time SAGE user, a little over a decade off and on. I have a variety of interests, but most of my coding time is spent working with matrix LU decompositions, finding Taylor series of "elementary" complex functions (logarithm or algebraic roots, e.g., log(x-c) or (x-c)^(2\*pi\*I/a), where c and a are complex numbers). And for years, I've been using the RealField and ComplexField classes, with their corresponding RealNumber and ComplexNumber elements, as the backbone of my code.
Just today, I discovered the MPComplexField class. I've been writing my own cython modules, borrowing code from the sage/rings/real_mpfr.pyx file. I noticed somewhere in my boilerplate code, that I was importing the MPC library. The MPC library appears to just be a C wrapper around the MPFR library, to perform complex arithmetic using pairs of MPFR numbers.
This struck me as a bit odd, because I've casually read through the source code in sage/rings/complex_number.pyx, and it looks to me like all the complex math is coded manually. E.g.,
a+b is calculated as (a.re + b.re) + I\*(a.im + b.im)
a*b is calculated as (a.re\*b.re - a.im\*b.im) + I\*(a.re\*b.im + a.im\*b.re)
sin(a) is a horrendous calculation involving sinh, square root, sine and cosine, as well as basic operations like multiplication, squaring, etc.
(Yes, I paraphrased the mpfr_xxx calls to make them readable.)
In other words, the ComplexNumber class implements the complex calculations manually, whereas the MPComplexNumber class lets the library handle those calculations. Calculating sin(a) is just a call to mpc_sin().
Which leads to my question. Why the two different classes? Which one is preferred as the more SAGE-ish approach? Is there a long-term plan to switch to MPC as the default in tutorials and code examples?
I ran some simple tests, and it appears that the rounding is not the same, at least not with default settings. After a few simple arithmetic calculations, I get diverging results, which I can only assume is due to a difference in rounding, and possibly due to differing internal algorithms. I've been using the ComplexField for over a decade, so I'm hesitant to switch. But I'm wondering if the MPComplexNumber has either better performance or better rounding, at least for most common applications.
Sorry for the long post, but this MPComplexNumber class appears to have been in SAGE for quite some time (since 2008). It seems like it must be the red-headed stepchild, because I don't seem to see it get mentioned anywhere. But if the performance is better, I'd gladly make the switch in all my coding projects. Any insight or advice?
Edit: My asterisk got turned into *italics*, which messed up my math sections. Escaped the asterisks, hope that fixes it.jaydfoxWed, 15 Jul 2020 20:27:57 +0200https://ask.sagemath.org/question/52478/Plotting questionhttps://ask.sagemath.org/question/7843/plotting-question/I am trying to plot Sage calculation listed below. I get the following error:
verbose 0 (3989: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points.
verbose 0 (3989: plot.py, generate_plot_points) Last error message: 'unable to simplify to float approximation'
What am I doing wrong?
Thank you for your help.
Gather all for Full H off resonance
Ix=1/2*matrix(SR,2,2,[[0,1],[1,0]])
Iy=1/(2*I)*matrix(SR,2,2,[[0,1],[-1,0]])
Iz=1/2*matrix(SR,2,2,[[1,0],[0,-1]])
Al=vector(SR,[1,0])
Be=vector(SR,[0,1])
g = 26.7522128e7 # 1H Gamma rad/s/T
rf_field = 100 # B1 rf field in kHz
Brf = rf_field*1e3*2*pi/g*2 # B1 in T for Wnut/2pi= 100 kHz
Thetarf = pi/2 # angle between BO and Brf
Wnut = (1/2*g*Brf*sin(Thetarf))
t360 = 1/((1/2*g*Brf*sin(Thetarf))/(2*pi))# 360 pulse
tip = 180 # Tip angle
tp = tip/360*t360
Bo = 500e6*2*pi/g# T Spectrometer Field Strength
wref = g*Bo # spectrometer reference frequency for Protons, rad/S
Bo.n(),(wref/(2*pi)).simplify()/1e6;(Wnut/(2*pi)).simplify()/1e3,tp.simplify()*1e6
Ratio = var ('Ratio')# Wnut/O0, Ratio creates offset, 0=on res
Boo=(Wnut*Ratio+wref)/g
woo = g*Boo
O0=woo-wref
Phip = var('Phip') # Phip = phase of rf pulse
Hrfrf = (O0)*Iz + Wnut*(Ix*cos(Phip)+Iy*sin(Phip))
Pabrf=((abs(Be*exp(-i*Hrfrf(Phip=0)*tp)*Al))^2)
for Ratio in srange (-1,1,0.25,include_endpoint=True):
Pabrf(Ratio=Ratio).n(digits=5) #(pi)x on Alpha off resonance Full H
plot (Pabrf,(Ratio, -1,1))mhfreyWed, 05 Jan 2011 17:20:55 +0100https://ask.sagemath.org/question/7843/Schwarz-Christoffel transformationshttps://ask.sagemath.org/question/48332/schwarz-christoffel-transformations/I would like to visualize the Schwarz-Christoffel transformations. Do you have an idea how to implement those? For example a polygon (as a domain) in the complex plane and subsequently also a map (decomposed into several steps) into the upper half-plane of the complex plane should be drawn. I am not quite sure how to do this. The formulas for these transformations are well-known.Tintin1Mon, 14 Oct 2019 15:27:31 +0200https://ask.sagemath.org/question/48332/can't display the inverse of complex matrixhttps://ask.sagemath.org/question/47326/cant-display-the-inverse-of-complex-matrix/ Hi
it seems that Sagemath is able to calculate this matrix inverse, but it is unable to display it.
Maybe someone will explain me why ? it is not easy to find what I do wrong because there is no error message for n=8
an it's ok for n=4
#FFT
w=var('w')
n=8
FFT=matrix(SR,n,n)
#show (FFT)
for i in (0..n-1):
for j in (0..n-1):
FFT[i,j]=w^(i*j)
FFTc=1/sqrt(n) *FFT.subs(w=e^(2*I*pi/n))
#show(" FFTc : ",FFTc)
FFTc.change_ring(CC)
show(" FFTc.right_kernel().matrix() : ",FFTc.right_kernel().matrix())
show (" FFTc.determinant() for n =",str(n),": ",FFTc.determinant())
FFTcI=copy(FFTc).inverse()
#show (" FFT w :",FFT)
show (" FFTc with w=e^(2*I*pi/",str(n),") : ",FFTc)
show (" FFTc inverse with w=e^(2*I*pi/",str(n),") : ",FFTcI)
ortolljWed, 31 Jul 2019 18:19:16 +0200https://ask.sagemath.org/question/47326/Possible bug in CC needs confirmationhttps://ask.sagemath.org/question/47068/possible-bug-in-cc-needs-confirmation/Hello, Sage community. I just noticed what I believe is a bug in the implementation of `CC`, but I would like to receive confirmation.
When I try the following:
'hello' in RR
I get the natural answer: `False`. However, if I try the same with `CC`:
'hello' in CC
I get the following error message:
NameError: name 'hello' is not defined
Just in case, I am using SageMath v8.8, specifically, the binary version downloaded from the mirrors.
Can somebody confirm this is a bug? Thanks in advance!dsejasFri, 05 Jul 2019 04:05:28 +0200https://ask.sagemath.org/question/47068/desolve (how can i enable numeric?)https://ask.sagemath.org/question/45937/desolve-how-can-i-enable-numeric/ Hi, I am sorry for asking this, but I have been looking all over the internet now and I can't find any solution to this.
When i am using `desolve` i sometimes get some crazy outputs which is not relevant for me at all.
For instance i am trying to desolve this equation
y' = 0.00054 \cdot y \cdot (500 - y)
So i white
k = var('k')
x = var('x')
y = function('y')(x)
Then i define my `de`
de = diff(y,x) == 0.00054*1*y*(500-y)
And i desolve like i have done before with many other equations without any problem
desolve(de,y,[0,90])
But i get this result
-100/27*log(y(x) - 500) + 100/27*log(y(x)) == -100/27*I*pi + x - 100/27*log(410) + 100/27*log(90)
The right result should be:
y = 500/(1+4.556*e^(-0.27*x))
Please help me. I have been readning in the doc section about desolve and many other places, but there is so much stuff i dont understand. I am on the edge to install a subsitute CAS program to help me everytime sagemath fails.Martin MårtenssonFri, 29 Mar 2019 11:25:04 +0100https://ask.sagemath.org/question/45937/desolve crazy outputhttps://ask.sagemath.org/question/45938/desolve-crazy-output/ I made another one where i was plotting a function and i got a similar problem
k = var('k')
x = var('x')
T = function('T')(x)
de = diff(T,x)==1.54-0.259*(T-22)
f = desolve(de,T,[0,22])
blob = 27
f.plot(x,0,10) + plot(blob,x,0,7.098,color="red",figsize=4) + points([7.098,27],color="green",size=100) + text("[7.098,27]",(7.098,27.5),color="black")
I dont have enough points to upload the image but the graph grows and flattens out and the other graph touches with the graph at x=7.098.
Then if i want this result calculated for me I am thinking that it can be done by solving. so...
solve(f==blob,x)
And i get this output
[x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(2/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(4/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(6/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(8/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(10/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(12/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(2/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(16/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(18/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(20/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(22/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(24/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(26/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(4/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(30/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(32/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(34/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(36/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(38/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(40/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(6/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(44/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(46/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(48/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(50/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(52/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(54/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(8/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(58/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(60/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(62/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(64/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(66/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(68/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(10/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(72/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(2/7*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(76/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(78/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(80/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(82/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(12/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(86/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(88/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(90/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(92/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(94/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(96/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(14/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(100/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(102/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(104/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(106/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(108/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(110/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(16/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(114/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(116/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(118/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(120/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(122/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(124/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(18/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(128/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(130/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(132/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(134/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(136/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(138/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(20/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(142/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(144/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(146/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(4/7*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(150/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(152/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(22/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(156/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(158/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(160/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(162/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(164/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(166/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(24/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(170/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(172/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(174/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(176/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(178/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(180/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(26/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(184/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(186/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(188/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(190/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(192/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(194/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(28/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(198/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(200/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(202/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(204/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(206/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(208/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(30/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(212/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(214/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(216/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(218/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(220/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(6/7*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(32/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(226/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(228/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(230/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(232/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(234/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(236/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(34/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(240/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(242/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(244/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(246/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(248/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(250/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(36/37*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(254/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(256/259*I*pi)), x == 1000*log(1/7*44^(1/259)*7^(258/259)*e^(258/259*I*pi)), x == -258000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -256000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -254000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -36000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -250000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -248000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -246000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -244000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -242000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -240000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -34000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -236000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -234000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -232000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -230000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -228000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -226000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -32000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -6000/7*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -220000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -218000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -216000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -214000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -212000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -30000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -208000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -206000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -204000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -202000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -200000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -198000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -28000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -194000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -192000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -190000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -188000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -186000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -184000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -26000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -180000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -178000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -176000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -174000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -172000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -170000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -24000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -166000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -164000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -162000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -160000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -158000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -156000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -22000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -152000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -150000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -4000/7*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -146000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -144000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -142000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -20000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -138000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -136000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -134000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -132000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -130000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -128000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -18000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -124000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -122000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -120000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -118000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -116000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -114000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -16000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -110000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -108000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -106000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -104000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -102000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -100000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -14000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -96000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -94000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -92000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -90000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -88000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -86000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -12000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -82000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -80000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -78000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -76000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -2000/7*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -72000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -10000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -68000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -66000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -64000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -62000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -60000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -58000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -8000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -54000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -52000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -50000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -48000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -46000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -44000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -6000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -40000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -38000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -36000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -34000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -32000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -30000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -4000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -26000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -24000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -22000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -20000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -18000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -16000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -2000/37*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -12000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -10000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -8000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -6000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -4000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == -2000/259*I*pi + 1000*log(1/7*44^(1/259)*7^(258/259)), x == 1000*log(1/7*44^(1/259)*7^(258/259))]
But i can see that the simple answer should just be 7.098 so how can i get sage to sell me that?
I am not interested in the complex or imaginary numbers i just want the simple numerical or rational answer from sage.
Martin MårtenssonFri, 29 Mar 2019 11:37:29 +0100https://ask.sagemath.org/question/45938/Is there a sage command to compute complex numbers?https://ask.sagemath.org/question/43617/is-there-a-sage-command-to-compute-complex-numbers/ I trying to have sage compute the sqrt(1/z) and 1/sqrt(z) of the complex expression z = 1+ I. Is there a sage command that will execute this?nevar123Sat, 08 Sep 2018 02:16:57 +0200https://ask.sagemath.org/question/43617/inverse_laplace of a fraction whose denominator has real roots (or complex)https://ask.sagemath.org/question/42331/inverse_laplace-of-a-fraction-whose-denominator-has-real-roots-or-complex/ Hi, this is the simplest and very important method of Laplace inversion . It can be applied to any Laplace transform
(LT), by starting with a Pade approximation of the LT, then partial fractions and inversion. For some reason, a program I had written last year stopped working
def Ruin(Fx, rho, x, u):#assumes rational survival function Fx, i.e. a hyperexponential density
var('s')
L_F=laplace(Fx,x,s)
m1=L_F(s=0) #some algebraic manipulations known as Pollaczek-Khinchine formula
fe=L_F/m1
Fe=factor((1-fe)/s)
L_rui=rho*Fe/(1-rho*fe)
C=ComplexField(53);
dec=Frac(C['s'])(L_rui).partial_fraction_decomposition();
n=len(dec)
par=[inverse_laplace(dec[1][i],s,u) for i in [0..n]];
rui=sum(par)
return rui, L_rui, fe
It's pretty easy to repair the program when roots are rational -- see question "partial_fraction_decomposition" with possibly "complex roots", again. For nonrational roots, I proposed there the following test case
var('s,u')
R.<s> = PolynomialRing(QQbar)
F = R.fraction_field()
L=3/4*(19*s^2 + 156*s + 284)/(19*s^3 + 174*s^2 + 422*s + 228)
whole,LL=L.partial_fraction_decomposition()
show(LL[0])
inverse_laplace(LL[0],s,u)
Thanks in advance :)florinMon, 14 May 2018 11:04:41 +0200https://ask.sagemath.org/question/42331/"partial_fraction_decomposition" with possibly "complex roots", againhttps://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/ Hi I come back to this question, even though it's been answered before, since I am still not able to make it work. I should mention maybe that I am just trying to teach undergraduate students to invert
rational Laplace transforms (for myself I am able to afford a Mathematica licence, but it would be nice to be able to show students that such simple things may be done nowadays for free). Following a previous answer, I tried
L = 2*(s + 3)/(3*s^2 + 13*s + 10)
Ks=FractionField(PolynomialRing(CC,names='s'))
Lr=Ks(L)
Of course, with quadratic rational roots, I could do this by hand, but the purpose here is to do it when you do not know the roots. The first two commands work, but the third has error message
TypeError: ('cannot convert {!r}/{!r} to an element of {}', 2*(s + 3)/(3*s^2 + 13*s + 10), 1.00000000000000, Fraction Field of Univariate Polynomial Ring in s over Complex Field with 53 bits of precision)
I should add that a different attempt which used to work last year
C=ComplexField(53);
dec=Frac(C['s'])(Lrs).partial_fraction_decomposition();
gets now same error message. So, this is probably due to an "improvement" of sageflorinMon, 30 Apr 2018 10:58:31 +0200https://ask.sagemath.org/question/42220/"partial_fraction_decomposition" with "complex roots"https://ask.sagemath.org/question/41893/partial_fraction_decomposition-with-complex-roots/ Hi six months ago I used a hint from Zimmerman's book
var('s')
L=2*(s + 3)/(3*s^2 + 13*s + 10)
C=ComplexField(53)
dec=Frac(C['s'])(L).partial_fraction_decomposition()
but this does not work anymore (see below). How to do this? Thanks
----> 2 dec=Frac(C['s'])(Lrs).partial_fraction_decomposition();
TypeError: ('cannot convert {!r}/{!r} to an element of {}', 2*(s + 3)/(3*s^2 + 13*s + 10), 1.00000000000000, Fraction Field of Univariate Polynomial Ring in s over Complex Field with 53 bits of precision)
/opt/sagemath-8.1/src/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9641)()
937 if mor is not None:
938 if no_extra_args:
--> 939 return mor._call_(x)
940 else:
941 return mor._call_with_args(x, args, kwds)
/opt/sagemath-8.1/src/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4928)()
152 print(type(C), C)
153 print(type(C._element_constructor), C._element_constructor)
--> 154 raise
155
156 cpdef Element _call_with_args(self, x, args=(), kwds={}):
/opt/sagemath-8.1/src/sage/structure/coerce_maps.pyx in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4796)()
147 cdef Parent C = self._codomain
148 try:
--> 149 return C._element_constructor(x)
150 except Exception:
151 if print_warnings:
/opt/sagemath-8.1/local/lib/python2.7/site-packages/sage/rings/fraction_field.py in _element_constructor_(self, x, y, coerce)
616 except AttributeError:
617 raise TypeError("cannot convert {!r}/{!r} to an element of {}",
--> 618 x0, y0, self)
619 try:
620 return self._element_class(self, x, y, coerce=coerce)
florinFri, 06 Apr 2018 21:42:27 +0200https://ask.sagemath.org/question/41893/How do I plot a real function whose computation involves complex intermediate results?https://ask.sagemath.org/question/41875/how-do-i-plot-a-real-function-whose-computation-involves-complex-intermediate-results/ I have an algebraic expressing **ev** for an eigenvalue of a 4x4 matrix and I want to plot it.
The expression has square roots and cube roots, and sometimes there are complex intermediate
results although the final result is real. With float, the imaginary part will be non-zero due to numeric errors, but I can just ignore it. Still, when I try to plot this function, I get error messages.
var("u")
ev = (1/2*u + 1/2*sqrt(2/3*u^2 - 1/9*(u^4 - 18*u)/(1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3) - (1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3) + 6/sqrt((u^4 +
3*(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3)*u^2 -
18*u + 9*(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) +
1/2)^(2/3))/(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) +
1/2)^(1/3))) + 1/6*sqrt((u^4 + 3*(1/27*u^6 - u^3 + 1/6*sqrt(-32/3*u^6 +
252*u^3 + 9) + 1/2)^(1/3)*u^2 - 18*u + 9*(1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(2/3))/(1/27*u^6 - u^3 +
1/6*sqrt(-32/3*u^6 + 252*u^3 + 9) + 1/2)^(1/3)))
print ev.subs(u=2.6)
print ev.subs(u=3.0)
print ev.subs(u=3.0).real()
plot(ev.subs(u=x).real(),(x,2.6,3))
Plotting will draw only from 2.6 to 2.62 and then abort with the following error message
verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 189 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: 'math domain error'
I tried this in the sage cell server. (On my outdated office computer installation, sage 7.4, I get a different error "AssertionError".)Günter RoteWed, 04 Apr 2018 11:46:06 +0200https://ask.sagemath.org/question/41875/Sage equivalent to Mathematica Chop function?https://ask.sagemath.org/question/41518/sage-equivalent-to-mathematica-chop-function/Does anyone know a function in Sage equivalent to Mathematica _Chop()_ function ( http://reference.wolfram.com/language/ref/Chop.html )? I couldn't find any.
In fact, I would like to drop the imaginary parts that are close to zero of complex polynomial coefficients.
This shouldn't be difficult to implement but is something that I use a lot, it is strange not find an implementation.
CordiallyjoaoffTue, 13 Mar 2018 11:17:06 +0100https://ask.sagemath.org/question/41518/numerical_integral (GSL C library) does not work with complex integrands?https://ask.sagemath.org/question/41510/numerical_integral-gsl-c-library-does-not-work-with-complex-integrands/I am unable to use numerical_integral (GSL C library) with a complex integrand. I haven't found anything about this issue in the documentation. Is this correct? nintegral works correctly.
e.g.
f(x) = (exp(i*x) + exp(-i*x))/2
f(x).nintegral(x,0,1)
numerical_integral(f,0,1)
the above code outputs
(0.8414709848078965, 9.34220461887732e-15, 21, 0)
Error in lines 3-3
...
TypeError: unable to coerce to a real number
Tested with SageMath version 7.5.1 and Cocalc.
CordiallyjoaoffMon, 12 Mar 2018 16:30:59 +0100https://ask.sagemath.org/question/41510/Getting all (complex) solutions of a non polynomial equationhttps://ask.sagemath.org/question/40302/getting-all-complex-solutions-of-a-non-polynomial-equation/Hi !
I was used to solve the following equation with Mathematica.
\begin{equation}
\alpha_1 + \alpha_2x + \alpha_3x^2 + x^4 + \frac{\alpha_4}{x^2-\alpha_0} + \frac{\alpha_5 x^2}{x^2-\alpha_0}=0
\end{equation}
where $\alpha_i$ are constants.
The Mathematica function "Solve" gives me all the numerical roots of this non polynomial equation very easily. These roots can be real or complex.
I'm a very beginner at Sage. I have tried several methods to solve this equation automatically but it seems that all methods I've used work only for polynomial equations. Here they are :
alpha0 = 0.25
alpha1 = -2.5
alpha2 = 6.9282
alpha3 = -5.5
alpha4 = 0.5
alpha5 = -0.5
x = var('x')
eq = alpha1 + alpha2*x + alpha3*x**2 + x**4 + alpha4/(x**2 - alpha0) + alpha5*x**2/(x**2 - alpha0) == 0.
# test 1
# solve(eq, x, ring=CC)
# ==> [0 == 20000*x^6 - 115000*x^4 + 138564*x^3 - 32500*x^2 - 34641*x + 22500]
# test 2
# solve(eq, x, ring=CC, solution_dict=True)
# ==> [{0: 20000*x^6 - 115000*x^4 + 138564*x^3 - 32500*x^2 - 34641*x + 22500}]
# test 3
# eq.roots(x, ring=CC,multiplicities=False)
# ==> TypeError: fraction must have unit denominator
Do you have any idea of a method to get the approximated roots of the equation ?
Thanks in advance
EDIT : correction of an error in the equation ; add few testsyoda789Tue, 26 Dec 2017 15:05:07 +0100https://ask.sagemath.org/question/40302/