| 1 | r""" |
|---|
| 2 | Interface to Maxima |
|---|
| 3 | |
|---|
| 4 | Maxima is a free GPL'd general purpose computer algebra system whose |
|---|
| 5 | development started in 1968 at MIT. It contains symbolic manipulation |
|---|
| 6 | algorithms, as well as implementations of special functions, including |
|---|
| 7 | elliptic functions and generalized hypergeometric functions. Moreover, |
|---|
| 8 | Maxima has implementations of many functions relating to the invariant |
|---|
| 9 | theory of the symmetric group $S_n$. (However, the commands for group |
|---|
| 10 | invariants, and the corresponding Maxima documenation, are in French.) |
|---|
| 11 | For many links to Maxima documentation see |
|---|
| 12 | \url{http://maxima.sourceforge.net/docs.shtml/}. |
|---|
| 13 | |
|---|
| 14 | AUTHORS OF THIS MODULE: |
|---|
| 15 | - William Stein (2005-12): Initial version |
|---|
| 16 | - David Joyner: Improved documentation |
|---|
| 17 | - William Stein (2005-01-08): Fixed bug in parsing |
|---|
| 18 | |
|---|
| 19 | If the string "error" (case insensitive) occurs in the output of |
|---|
| 20 | anything from maxima, a RuntimeError exception is raised. |
|---|
| 21 | |
|---|
| 22 | EXAMPLES: |
|---|
| 23 | We evaluate a very simple expression in maxima. |
|---|
| 24 | sage: maxima('3 * 5') |
|---|
| 25 | 15 |
|---|
| 26 | |
|---|
| 27 | We factor $x^5 - y^5$ in Maxima in several different ways. |
|---|
| 28 | The first way yields a Maxima object. |
|---|
| 29 | sage: F = maxima.factor('x^5 - y^5') |
|---|
| 30 | sage: F |
|---|
| 31 | -(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4) |
|---|
| 32 | sage: type(F) |
|---|
| 33 | <class 'sage.interfaces.maxima.MaximaElement'> |
|---|
| 34 | |
|---|
| 35 | Note that Maxima objects can also be displayed using ``ASCII art''; |
|---|
| 36 | to see a normal linear representation of any Maxima object x, |
|---|
| 37 | use \code{str(x)}. |
|---|
| 38 | sage: F.display2d() |
|---|
| 39 | 4 3 2 2 3 4 |
|---|
| 40 | - (y - x) (y + x y + x y + x y + x ) |
|---|
| 41 | |
|---|
| 42 | We can make this the default: |
|---|
| 43 | sage: maxima.display2d(True) |
|---|
| 44 | sage: F |
|---|
| 45 | 4 3 2 2 3 4 |
|---|
| 46 | - (y - x) (y + x y + x y + x y + x ) |
|---|
| 47 | |
|---|
| 48 | You can always use \code{x.str()} to obtain the linear representation |
|---|
| 49 | of an object, even without changing the display2d flag. This can |
|---|
| 50 | be useful for moving maxima data to other systems. |
|---|
| 51 | sage: F.str() |
|---|
| 52 | '-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)' |
|---|
| 53 | |
|---|
| 54 | sage: maxima.display2d(False) |
|---|
| 55 | sage: F |
|---|
| 56 | -(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4) |
|---|
| 57 | |
|---|
| 58 | |
|---|
| 59 | The \code{maxima.eval} command evaluates an expression in maxima |
|---|
| 60 | and returns the result as a string. |
|---|
| 61 | |
|---|
| 62 | sage: print maxima.eval('factor(x^5 - y^5)') |
|---|
| 63 | -(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4) |
|---|
| 64 | |
|---|
| 65 | We can create the polynomial $f$ as a Maxima polynomial, then call |
|---|
| 66 | the factor method on it. Notice that the notation \code{f.factor()} |
|---|
| 67 | is consistent with how the rest of \sage works. |
|---|
| 68 | sage: f = maxima('x^5 - y^5') |
|---|
| 69 | sage: f^2 |
|---|
| 70 | (x^5 - y^5)^2 |
|---|
| 71 | sage: f.factor() |
|---|
| 72 | -(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4) |
|---|
| 73 | |
|---|
| 74 | Control-C interruption works well with the maxima interface, |
|---|
| 75 | because of the excellent implementation of maxima. For example, |
|---|
| 76 | try the following sum but with a much bigger range, and hit |
|---|
| 77 | control-C. |
|---|
| 78 | sage: maxima('sum(1/x^2, x, 1, 10)') |
|---|
| 79 | 1968329/1270080 |
|---|
| 80 | |
|---|
| 81 | \subsection{Tutorial} |
|---|
| 82 | We follow the tutorial at |
|---|
| 83 | \url{http://maxima.sourceforge.net/docs/intromax/}. |
|---|
| 84 | |
|---|
| 85 | sage: maxima('1/100 + 1/101') |
|---|
| 86 | 201/10100 |
|---|
| 87 | |
|---|
| 88 | sage: a = maxima('(1 + sqrt(2))^5'); a |
|---|
| 89 | (sqrt(2) + 1)^5 |
|---|
| 90 | sage: a.expand() |
|---|
| 91 | 29*sqrt(2) + 41 |
|---|
| 92 | |
|---|
| 93 | sage: a = maxima('(1 + sqrt(2))^5') |
|---|
| 94 | sage: float(a) |
|---|
| 95 | 82.012193308819747 |
|---|
| 96 | sage: a.numer() |
|---|
| 97 | 82.01219330881975 |
|---|
| 98 | |
|---|
| 99 | sage: maxima.eval('fpprec : 100') |
|---|
| 100 | '100' |
|---|
| 101 | sage: a.bfloat() |
|---|
| 102 | 8.20121933088197564152489730020812442785204843859314941221237124017312418754011041266612384955016056B1 |
|---|
| 103 | |
|---|
| 104 | sage: maxima('100!') |
|---|
| 105 | 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 |
|---|
| 106 | |
|---|
| 107 | sage: f = maxima('(x + 3*y + x^2*y)^3') |
|---|
| 108 | sage: f.expand() |
|---|
| 109 | x^6*y^3 + 9*x^4*y^3 + 27*x^2*y^3 + 27*y^3 + 3*x^5*y^2 + 18*x^3*y^2 + 27*x*y^2 + 3*x^4*y + 9*x^2*y + x^3 |
|---|
| 110 | sage: f.subst('x=5/z') |
|---|
| 111 | (5/z + 25*y/z^2 + 3*y)^3 |
|---|
| 112 | sage: g = f.subst('x=5/z') |
|---|
| 113 | sage: h = g.ratsimp(); h |
|---|
| 114 | (27*y^3*z^6 + 135*y^2*z^5 + (675*y^3 + 225*y)*z^4 + (2250*y^2 + 125)*z^3 + (5625*y^3 + 1875*y)*z^2 + 9375*y^2*z + 15625*y^3)/z^6 |
|---|
| 115 | sage: h.factor() |
|---|
| 116 | (3*y*z^2 + 5*z + 25*y)^3/z^6 |
|---|
| 117 | |
|---|
| 118 | sage: eqn = maxima(['a+b*c=1', 'b-a*c=0', 'a+b=5']) |
|---|
| 119 | sage: s = eqn.solve('[a,b,c]'); s |
|---|
| 120 | [[a = (25*sqrt(79)*%i + 25)/(6*sqrt(79)*%i - 34),b = (5*sqrt(79)*%i + 5)/(sqrt(79)*%i + 11),c = (sqrt(79)*%i + 1)/10],[a = (25*sqrt(79)*%i - 25)/(6*sqrt(79)*%i + 34),b = (5*sqrt(79)*%i - 5)/(sqrt(79)*%i - 11),c = - (sqrt(79)*%i - 1)/10]] |
|---|
| 121 | |
|---|
| 122 | Here is an example of solving an algebraic equation: |
|---|
| 123 | sage: maxima('x^2+y^2=1').solve('y') |
|---|
| 124 | [y = - sqrt(1 - x^2),y = sqrt(1 - x^2)] |
|---|
| 125 | sage: maxima('x^2 + y^2 = (x^2 - y^2)/sqrt(x^2 + y^2)').solve('y') |
|---|
| 126 | [y = - sqrt(( - y^2 - x^2)*sqrt(y^2 + x^2) + x^2),y = sqrt(( - y^2 - x^2)*sqrt(y^2 + x^2) + x^2)] |
|---|
| 127 | |
|---|
| 128 | You can even nicely typeset the solution in latex: |
|---|
| 129 | sage: print latex(s) |
|---|
| 130 | \left[\left[a = \frac{25\sqrt{79}i + 25}{6\sqrt{79}i - 34},b = \frac{5\sqrt{79}i + 5}{\sqrt{79}i + 11},c = \frac{\sqrt{79}i + 1}{10}\right],\left[a = \frac{25\sqrt{79}i - 25}{6\sqrt{79}i + 34},b = \frac{5\sqrt{79}i - 5}{\sqrt{79}i - 11},c = - \frac{\sqrt{79}i - 1}{10}\right]\right] |
|---|
| 131 | |
|---|
| 132 | To have the above appear onscreen via \code{xdvi}, type \code{view(s)}. |
|---|
| 133 | (TODO: For OS X should create pdf output and use preview instead?) |
|---|
| 134 | |
|---|
| 135 | sage: e = maxima('sin(u + v) * cos(u)^3'); e |
|---|
| 136 | cos(u)^3*sin(v + u) |
|---|
| 137 | sage: f = e.trigexpand(); f |
|---|
| 138 | cos(u)^3*(cos(u)*sin(v) + sin(u)*cos(v)) |
|---|
| 139 | sage: f.trigreduce() |
|---|
| 140 | (sin(v + 4*u) + sin(v - 2*u))/8 + (3*sin(v + 2*u) + 3*sin(v))/8 |
|---|
| 141 | sage: w = maxima('3 + k*%i') |
|---|
| 142 | sage: f = w^2 + maxima('%e')^w |
|---|
| 143 | sage: f.realpart() |
|---|
| 144 | %e^3*cos(k) - k^2 + 9 |
|---|
| 145 | |
|---|
| 146 | sage: f = maxima('x^3 * %e^(k*x) * sin(w*x)'); f |
|---|
| 147 | x^3*%e^(k*x)*sin(w*x) |
|---|
| 148 | sage: f.diff('x') |
|---|
| 149 | k*x^3*%e^(k*x)*sin(w*x) + 3*x^2*%e^(k*x)*sin(w*x) + w*x^3*%e^(k*x)*cos(w*x) |
|---|
| 150 | sage: f.integrate('x') |
|---|
| 151 | (((k*w^6 + 3*k^3*w^4 + 3*k^5*w^2 + k^7)*x^3 + (3*w^6 + 3*k^2*w^4 - 3*k^4*w^2 - 3*k^6)*x^2 + ( - 18*k*w^4 - 12*k^3*w^2 + 6*k^5)*x - 6*w^4 + 36*k^2*w^2 - 6*k^4)*%e^(k*x)*sin(w*x) + (( - w^7 - 3*k^2*w^5 - 3*k^4*w^3 - k^6*w)*x^3 + (6*k*w^5 + 12*k^3*w^3 + 6*k^5*w)*x^2 + (6*w^5 - 12*k^2*w^3 - 18*k^4*w)*x - 24*k*w^3 + 24*k^3*w)*%e^(k*x)*cos(w*x))/(w^8 + 4*k^2*w^6 + 6*k^4*w^4 + 4*k^6*w^2 + k^8) |
|---|
| 152 | |
|---|
| 153 | sage: f = maxima('1/x^2') |
|---|
| 154 | sage: f.integrate('x', 1, 'inf') |
|---|
| 155 | 1 |
|---|
| 156 | sage: g = maxima('f/sinh(k*x)^4') |
|---|
| 157 | sage: g.taylor('x', 0, 3) |
|---|
| 158 | f/(k^4*x^4) - 2*f/(3*k^2*x^2) + 11*f/45 - 62*k^2*f*x^2/945 |
|---|
| 159 | |
|---|
| 160 | \subsection{Examples involving matrices} |
|---|
| 161 | We illustrate computing with the matrix whose $i,j$ entry |
|---|
| 162 | is $i/j$, for $i,j=1,\ldots,4$. |
|---|
| 163 | |
|---|
| 164 | sage: f = maxima.eval('f[i,j] := i/j') |
|---|
| 165 | sage: A = maxima('genmatrix(f,4,4)'); A |
|---|
| 166 | matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1]) |
|---|
| 167 | sage: A.determinant() |
|---|
| 168 | 0 |
|---|
| 169 | sage: A.echelon() |
|---|
| 170 | matrix([1,1/2,1/3,1/4],[0,0,0,0],[0,0,0,0],[0,0,0,0]) |
|---|
| 171 | sage: A.eigenvalues() |
|---|
| 172 | [[0,4],[3,1]] |
|---|
| 173 | sage: A.eigenvectors() |
|---|
| 174 | [[[0,4],[3,1]],[1,0,0, - 4],[0,1,0, - 2],[0,0,1, - 4/3],[1,2,3,4]] |
|---|
| 175 | |
|---|
| 176 | We can also compute the echelon form in \sage: |
|---|
| 177 | sage: B = matrix(A, QQ) |
|---|
| 178 | sage: B.echelon_form() |
|---|
| 179 | [ 1 1/2 1/3 1/4] |
|---|
| 180 | [ 0 0 0 0] |
|---|
| 181 | [ 0 0 0 0] |
|---|
| 182 | [ 0 0 0 0] |
|---|
| 183 | sage: B.charpoly().factor() |
|---|
| 184 | (x - 4) * x^3 |
|---|
| 185 | |
|---|
| 186 | \subsection{Laplace Transforms} |
|---|
| 187 | We illustrate Laplace transforms: |
|---|
| 188 | sage: _ = maxima.eval("f(t) := t*sin(t)") |
|---|
| 189 | sage: maxima("laplace(f(t),t,s)") |
|---|
| 190 | 2*s/(s^2 + 1)^2 |
|---|
| 191 | |
|---|
| 192 | sage: maxima("laplace(delta(t-3),t,s)") #Dirac delta function |
|---|
| 193 | %e^ - (3*s) |
|---|
| 194 | |
|---|
| 195 | sage: _ = maxima.eval("f(t) := exp(t)*sin(t)") |
|---|
| 196 | sage: maxima("laplace(f(t),t,s)") |
|---|
| 197 | 1/(s^2 - 2*s + 2) |
|---|
| 198 | |
|---|
| 199 | sage: _ = maxima.eval("f(t) := t^5*exp(t)*sin(t)") |
|---|
| 200 | sage: maxima("laplace(f(t),t,s)") |
|---|
| 201 | 360*(2*s - 2)/(s^2 - 2*s + 2)^4 - 480*(2*s - 2)^3/(s^2 - 2*s + 2)^5 + 120*(2*s - 2)^5/(s^2 - 2*s + 2)^6 |
|---|
| 202 | sage: maxima("laplace(f(t),t,s)").display2d() |
|---|
| 203 | 3 5 |
|---|
| 204 | 360 (2 s - 2) 480 (2 s - 2) 120 (2 s - 2) |
|---|
| 205 | --------------- - --------------- + --------------- |
|---|
| 206 | 2 4 2 5 2 6 |
|---|
| 207 | (s - 2 s + 2) (s - 2 s + 2) (s - 2 s + 2) |
|---|
| 208 | |
|---|
| 209 | sage: maxima("laplace(diff(x(t),t),t,s)") |
|---|
| 210 | s*laplace(x(t),t,s) - x(0) |
|---|
| 211 | |
|---|
| 212 | sage: maxima("laplace(diff(x(t),t,2),t,s)") |
|---|
| 213 | -at('diff(x(t),t,1),t = 0) + s^2*laplace(x(t),t,s) - x(0)*s |
|---|
| 214 | |
|---|
| 215 | It is difficult to read some of these without the 2d representation: |
|---|
| 216 | sage.: maxima("laplace(diff(x(t),t,2),t,s)").display2d() |
|---|
| 217 | ! |
|---|
| 218 | d ! 2 |
|---|
| 219 | - -- (x(t))! + s laplace(x(t), t, s) - x(0) s |
|---|
| 220 | dt ! |
|---|
| 221 | !t = 0 |
|---|
| 222 | |
|---|
| 223 | Even better, use \code{view(maxima("laplace(diff(x(t),t,2),t,s)"))} to see |
|---|
| 224 | a typeset version. |
|---|
| 225 | |
|---|
| 226 | \subsection{Continued Fractions} |
|---|
| 227 | |
|---|
| 228 | A continued fraction $a + 1/(b + 1/(c + \cdots))$ is |
|---|
| 229 | represented in maxima by the list $[a, b, c, \ldots]$. |
|---|
| 230 | |
|---|
| 231 | sage: maxima("cf((1 + sqrt(5))/2)") |
|---|
| 232 | [1,1,1,1,2] |
|---|
| 233 | sage: maxima("cf ((1 + sqrt(341))/2)") |
|---|
| 234 | [9,1,2,1,2,1,17,1,2,1,2,1,17,1,2,1,2,1,17,2] |
|---|
| 235 | |
|---|
| 236 | \subsection{Special examples} |
|---|
| 237 | |
|---|
| 238 | In this section we illustrate calculations that would be awkward to do |
|---|
| 239 | (as far as I know) in non-symbolic computer algebra systems like MAGMA |
|---|
| 240 | or GAP. |
|---|
| 241 | |
|---|
| 242 | We compute the gcd of $2x^{n+4} - x^{n+2}$ and $4x^{n+1} + 3x^n$ |
|---|
| 243 | for arbitrary $n$. |
|---|
| 244 | |
|---|
| 245 | sage: f = maxima('2*x^(n+4) - x^(n+2)') |
|---|
| 246 | sage: g = maxima('4*x^(n+1) + 3*x^n') |
|---|
| 247 | sage: f.gcd(g) |
|---|
| 248 | x^n |
|---|
| 249 | |
|---|
| 250 | You can plot 3d graphs (via gnuplot): |
|---|
| 251 | |
|---|
| 252 | sage.: maxima('plot3d(x^2-y^2, [x,-2,2], [y,-2,2], [grid,12,12])') |
|---|
| 253 | [displays a 3 dimensional graph] |
|---|
| 254 | |
|---|
| 255 | You can formally evaluate sums (note the \code{nusum} command): |
|---|
| 256 | |
|---|
| 257 | sage: S = maxima('nusum(exp(1+2*i/n),i,1,n)') |
|---|
| 258 | sage.: S.display2d() |
|---|
| 259 | 2/n + 3 2/n + 1 |
|---|
| 260 | %e %e |
|---|
| 261 | ----------------------- - ----------------------- |
|---|
| 262 | 1/n 1/n 1/n 1/n |
|---|
| 263 | (%e - 1) (%e + 1) (%e - 1) (%e + 1) |
|---|
| 264 | |
|---|
| 265 | We formally compute the limit as $n\to\infty$ of $2S/n$ as follows: |
|---|
| 266 | |
|---|
| 267 | sage: T = S*maxima('2/n') |
|---|
| 268 | sage: T.tlimit('n','inf') |
|---|
| 269 | %e^3 - %e |
|---|
| 270 | |
|---|
| 271 | \subsection{Miscellaneous} |
|---|
| 272 | Obtaining digits of $\pi$: |
|---|
| 273 | sage: maxima.eval('fpprec : 100') |
|---|
| 274 | '100' |
|---|
| 275 | sage: maxima(pi).bfloat() |
|---|
| 276 | 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068B0 |
|---|
| 277 | |
|---|
| 278 | Defining functions in maxima: |
|---|
| 279 | sage: maxima.eval('fun[a] := a^2') |
|---|
| 280 | 'fun[a] := a^2' |
|---|
| 281 | sage: maxima('fun[10]') |
|---|
| 282 | 100 |
|---|
| 283 | |
|---|
| 284 | \subsection{Interactivity} |
|---|
| 285 | Unfortunately maxima doesn't seem to have a non-interactive mode, |
|---|
| 286 | which is needed for the \sage interface. If any \sage call leads |
|---|
| 287 | to maxima interactively answering questions, then the questions |
|---|
| 288 | can't be answered and the maxima session may hang. |
|---|
| 289 | See the discussion at \url{http://www.ma.utexas.edu/pipermail/maxima/2005/011061.html} for some ideas about how to fix this problem. An |
|---|
| 290 | example that illustrates this problem is |
|---|
| 291 | \code{maxima.eval('integrate (exp(a*x), x, 0, inf)')}. |
|---|
| 292 | |
|---|
| 293 | \subsection{Latex Output} |
|---|
| 294 | The latex output of Maxima is not perfect. E.g., |
|---|
| 295 | |
|---|
| 296 | sage: maxima.eval('tex(sin(u) + sinh(v^2))') |
|---|
| 297 | '$$\\sinhv^2 + \\sinu$$false' |
|---|
| 298 | |
|---|
| 299 | Notice the lack of space after the sin macro, which is a latex syntax |
|---|
| 300 | error. In \sage this is automatically fixed via a substition for |
|---|
| 301 | trig functions, which may have potentially bad side effects: |
|---|
| 302 | |
|---|
| 303 | sage: latex(maxima('sin(u) + sinh(v^2)')) |
|---|
| 304 | '\\sin{}hv^2 + \\sin{}u' |
|---|
| 305 | |
|---|
| 306 | It would be nice if somebody would fix this problem. One way would |
|---|
| 307 | be to improve Maxima by making the fix to Maxima and giving this back |
|---|
| 308 | to the Maxima people. |
|---|
| 309 | |
|---|
| 310 | \subsection{Long Input} |
|---|
| 311 | The MAXIMA interface reads in even very long input (using files) in a |
|---|
| 312 | robust manner, as long as you are creating a new object. |
|---|
| 313 | \note{Using \code{maxima.eval} for long input |
|---|
| 314 | is much less robust, and is not recommended.} |
|---|
| 315 | |
|---|
| 316 | sage: t = '"%s"'%10^10000 # ten thousand character string. |
|---|
| 317 | sage: a = maxima(t) |
|---|
| 318 | """ |
|---|
| 319 | |
|---|
| 320 | #***************************************************************************** |
|---|
| 321 | # Copyright (C) 2005 William Stein <wstein@ucsd.edu> |
|---|
| 322 | # |
|---|
| 323 | # Distributed under the terms of the GNU General Public License (GPL) |
|---|
| 324 | # |
|---|
| 325 | # This code is distributed in the hope that it will be useful, |
|---|
| 326 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |
|---|
| 327 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|---|
| 328 | # General Public License for more details. |
|---|
| 329 | # |
|---|
| 330 | # The full text of the GPL is available at: |
|---|
| 331 | # |
|---|
| 332 | # http://www.gnu.org/licenses/ |
|---|
| 333 | #***************************************************************************** |
|---|
| 334 | |
|---|
| 335 | import os, re |
|---|
| 336 | |
|---|
| 337 | from expect import Expect, ExpectElement, tmp |
|---|
| 338 | |
|---|
| 339 | from sage.misc.misc import verbose |
|---|
| 340 | |
|---|
| 341 | SAGE_START = '_s_start_' |
|---|
| 342 | SAGE_END = '__s_stop_' |
|---|
| 343 | cnt = 0 |
|---|
| 344 | |
|---|
| 345 | # The Maxima "apropos" command, e.g., apropos(det) gives a list |
|---|
| 346 | # of all identifiers that begin in a certain way. This could |
|---|
| 347 | # maybe be useful somehow... (?) Also maxima has a lot for getting |
|---|
| 348 | # documentation from the system -- this could also be useful. |
|---|
| 349 | |
|---|
| 350 | class Maxima(Expect): |
|---|
| 351 | """ |
|---|
| 352 | Interface to the Maxima interpreter. |
|---|
| 353 | """ |
|---|
| 354 | def __call__(self, x): |
|---|
| 355 | import sage.rings.all |
|---|
| 356 | if sage.rings.all.is_Infinity(x): |
|---|
| 357 | return Expect.__call__(self, 'inf') |
|---|
| 358 | else: |
|---|
| 359 | return Expect.__call__(self, x) |
|---|
| 360 | |
|---|
| 361 | def __init__(self, script_subdirectory=None, logfile=None, server=None): |
|---|
| 362 | """ |
|---|
| 363 | Create an instance of the Maxima interpreter. |
|---|
| 364 | """ |
|---|
| 365 | # TODO: Input and output prompts in maxima can be changed by |
|---|
| 366 | # setting inchar and outchar.. |
|---|
| 367 | eval_using_file_cutoff = 200 |
|---|
| 368 | self.__eval_using_file_cutoff = eval_using_file_cutoff |
|---|
| 369 | Expect.__init__(self, |
|---|
| 370 | name = 'maxima', |
|---|
| 371 | prompt = '\(\%i[0-9]+\)', |
|---|
| 372 | command = "maxima", |
|---|
| 373 | maxread = 1, # CRUCIAL to use less buffering for maxima (or get all kinds of hangs on OS X and 64-bit machines, etc! |
|---|
| 374 | script_subdirectory = script_subdirectory, |
|---|
| 375 | restart_on_ctrlc = False, |
|---|
| 376 | verbose_start = False, |
|---|
| 377 | init_code = ['display2d : false', # no ascii art output |
|---|
| 378 | 'load("mactex-utilities")' # latex instead of plain tex from tex command |
|---|
| 379 | ], |
|---|
| 380 | logfile = logfile, |
|---|
| 381 | eval_using_file_cutoff=eval_using_file_cutoff ) |
|---|
| 382 | self._display2d = False |
|---|
| 383 | |
|---|
| 384 | def _start(self): |
|---|
| 385 | # For some reason sending a single input line at startup avoids |
|---|
| 386 | # lots of weird timing issues when doing doctests. |
|---|
| 387 | Expect._start(self) |
|---|
| 388 | self(1) |
|---|
| 389 | |
|---|
| 390 | # this doesn't work. |
|---|
| 391 | #def x_start(self): |
|---|
| 392 | # Expect._start(self) |
|---|
| 393 | # self._expect.sendline('inchar:"__SAGE__";') |
|---|
| 394 | # self._change_prompt('__SAGE__[0-9]+\)') |
|---|
| 395 | # self.expect().expect('__SAGE__[0-9]+\)') |
|---|
| 396 | |
|---|
| 397 | def _eval_line_using_file(self, line, tmp): |
|---|
| 398 | F = open(tmp, 'w') |
|---|
| 399 | F.write(line) |
|---|
| 400 | F.close() |
|---|
| 401 | if self._expect is None: |
|---|
| 402 | self._start() |
|---|
| 403 | self._expect.sendline('batchload("%s");'%tmp) |
|---|
| 404 | self._expect.expect(self._prompt) |
|---|
| 405 | return '' |
|---|
| 406 | |
|---|
| 407 | def __reduce__(self): |
|---|
| 408 | return reduce_load_Maxima, tuple([]) |
|---|
| 409 | |
|---|
| 410 | def _quit_string(self): |
|---|
| 411 | return 'quit();' |
|---|
| 412 | |
|---|
| 413 | def _eval_line(self, line, reformat=True, allow_use_file=False): |
|---|
| 414 | line = line.rstrip().rstrip(';') |
|---|
| 415 | if line == '': |
|---|
| 416 | return '' |
|---|
| 417 | line = '%s; %s; %s;'%(SAGE_START, line, SAGE_END) |
|---|
| 418 | if self._expect is None: |
|---|
| 419 | self._start() |
|---|
| 420 | if allow_use_file and self.__eval_using_file_cutoff and \ |
|---|
| 421 | len(line) > self.__eval_using_file_cutoff: |
|---|
| 422 | return self._eval_line_using_file(line, tmp) |
|---|
| 423 | try: |
|---|
| 424 | E = self._expect |
|---|
| 425 | #print "in = '%s'"%line |
|---|
| 426 | E.sendline(line) |
|---|
| 427 | self._expect.expect(SAGE_END) |
|---|
| 428 | self._expect.expect(SAGE_END) |
|---|
| 429 | out = self._expect.before |
|---|
| 430 | #print "out = '%s'"%out |
|---|
| 431 | self._expect.expect(self._prompt) |
|---|
| 432 | out += self._expect.before |
|---|
| 433 | |
|---|
| 434 | except KeyboardInterrupt: |
|---|
| 435 | self._keyboard_interrupt() |
|---|
| 436 | |
|---|
| 437 | i = out.rfind(SAGE_START) |
|---|
| 438 | j = out.rfind(SAGE_END) |
|---|
| 439 | out = out[i+len(SAGE_START):j] |
|---|
| 440 | if not reformat: |
|---|
| 441 | return out |
|---|
| 442 | if out.find('error') != -1: |
|---|
| 443 | return out |
|---|
| 444 | out = out.lstrip() |
|---|
| 445 | #i = out.rfind('(') |
|---|
| 446 | #out = out[:i].strip() |
|---|
| 447 | i = out.find('(%o') |
|---|
| 448 | out0 = out[:i].strip() |
|---|
| 449 | i += out[i:].find(')') |
|---|
| 450 | out1 = out[i+1:].strip() |
|---|
| 451 | out = out0 + out1 |
|---|
| 452 | out = ''.join(out.split()) # no whitespace |
|---|
| 453 | out = out.replace('-',' - ').replace('+',' + ').replace('=',' = ').replace(': =',' :=') |
|---|
| 454 | if out[:3] == ' - ': |
|---|
| 455 | out = '-' + out[3:] |
|---|
| 456 | i = out.rfind('(%o') |
|---|
| 457 | if i != -1: |
|---|
| 458 | out = out[:i] |
|---|
| 459 | return out |
|---|
| 460 | |
|---|
| 461 | def _object_class(self): |
|---|
| 462 | return MaximaElement |
|---|
| 463 | |
|---|
| 464 | def set(self, var, value): |
|---|
| 465 | """ |
|---|
| 466 | Set the variable var to the given value. |
|---|
| 467 | """ |
|---|
| 468 | cmd = '%s : %s;'%(var, value) |
|---|
| 469 | #out = self._eval_line(cmd, reformat=False) |
|---|
| 470 | out = self._eval_line(cmd, reformat=False, allow_use_file=True) |
|---|
| 471 | |
|---|
| 472 | if out.find("error") != -1: |
|---|
| 473 | raise TypeError, "Error executing code in Maxima\nCODE:\n\t%s\nMaxima ERROR:\n\t%s"%(cmd, out) |
|---|
| 474 | |
|---|
| 475 | |
|---|
| 476 | def get(self, var): |
|---|
| 477 | """ |
|---|
| 478 | Get the string value of the variable var. |
|---|
| 479 | """ |
|---|
| 480 | s = self._eval_line('%s'%var) |
|---|
| 481 | return s |
|---|
| 482 | |
|---|
| 483 | #def clear(self, var): |
|---|
| 484 | # """ |
|---|
| 485 | # Clear the variable named var. |
|---|
| 486 | # """ |
|---|
| 487 | # if self._expect is None: |
|---|
| 488 | # return |
|---|
| 489 | # self._expect.sendline('kill(%s);'%var) |
|---|
| 490 | # self._expect.expect(self._prompt) |
|---|
| 491 | |
|---|
| 492 | def console(self): |
|---|
| 493 | maxima_console() |
|---|
| 494 | |
|---|
| 495 | def version(self): |
|---|
| 496 | return maxima_version() |
|---|
| 497 | |
|---|
| 498 | def display2d(self, flag=True): |
|---|
| 499 | """ |
|---|
| 500 | Set the flag that determines whether Maxima objects are |
|---|
| 501 | printed using their 2-d ASCII art representation. When the |
|---|
| 502 | maxima interface starts the default is that objects are not |
|---|
| 503 | represented in 2-d. |
|---|
| 504 | |
|---|
| 505 | INPUT: |
|---|
| 506 | flag -- bool (default: True) |
|---|
| 507 | |
|---|
| 508 | EXAMPLES |
|---|
| 509 | sage: maxima('1/2') |
|---|
| 510 | 1/2 |
|---|
| 511 | sage: maxima.display2d(True) |
|---|
| 512 | sage: maxima('1/2') |
|---|
| 513 | 1 |
|---|
| 514 | - |
|---|
| 515 | 2 |
|---|
| 516 | sage: maxima.display2d(False) |
|---|
| 517 | """ |
|---|
| 518 | self._display2d = bool(flag) |
|---|
| 519 | |
|---|
| 520 | def plot2d(self, *args): |
|---|
| 521 | r""" |
|---|
| 522 | Plot a 2d graph using Maxima / gnuplot. |
|---|
| 523 | |
|---|
| 524 | maxima.plot2d(f, '[var, min, max]', options) |
|---|
| 525 | |
|---|
| 526 | INPUT: |
|---|
| 527 | f -- a string representing a function (such as f="sin(x)") |
|---|
| 528 | [var, xmin, xmax] |
|---|
| 529 | options -- an optional string representing plot2d options in gnuplot format |
|---|
| 530 | |
|---|
| 531 | EXAMPLES: |
|---|
| 532 | sage.: maxima.plot2d('sin(x)','[x,-5,5]') |
|---|
| 533 | sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "sin-plot.eps"]' |
|---|
| 534 | sage.: maxima.plot2d('sin(x)','[x,-5,5]',opts) |
|---|
| 535 | |
|---|
| 536 | The eps file is saved in the current directory. |
|---|
| 537 | """ |
|---|
| 538 | self('plot2d(%s)'%(','.join([str(x) for x in args]))) |
|---|
| 539 | |
|---|
| 540 | def plot2d_parametric(self, r, var, trange, nticks=50, options=None): |
|---|
| 541 | r""" |
|---|
| 542 | Plots r = [x(t), y(t)] for t = tmin...tmax using gnuplot with options |
|---|
| 543 | |
|---|
| 544 | INPUT: |
|---|
| 545 | r -- a string representing a function (such as r="[x(t),y(t)]") |
|---|
| 546 | var -- a string representing the variable (such as var = "t") |
|---|
| 547 | trange -- [tmin, tmax] are numbers with tmin<tmax |
|---|
| 548 | nticks -- int (default: 50) |
|---|
| 549 | options -- an optional string representing plot2d options in gnuplot format |
|---|
| 550 | |
|---|
| 551 | EXAMPLES: |
|---|
| 552 | sage.: maxima.plot2d_parametric(["sin(t)","cos(t)"], "t",[-3.1,3.1]) |
|---|
| 553 | |
|---|
| 554 | sage.: opts = '[gnuplot_preamble, "set nokey"], [gnuplot_term, ps], [gnuplot_out_file, "circle-plot.eps"]' |
|---|
| 555 | sage.: maxima.plot2d_parametric(["sin(t)","cos(t)"], "t", [-3.1,3.1], options=opts) |
|---|
| 556 | |
|---|
| 557 | The eps file is saved to the current working directory. |
|---|
| 558 | |
|---|
| 559 | Here is another fun plot: |
|---|
| 560 | sage.: maxima.plot2d_parametric(["sin(5*t)","cos(11*t)"], "t", [0,2*pi()], nticks=400) |
|---|
| 561 | """ |
|---|
| 562 | tmin = trange[0] |
|---|
| 563 | tmax = trange[1] |
|---|
| 564 | cmd = "plot2d([parametric, %s, %s, [%s, %s, %s], [nticks, %s]]"%( \ |
|---|
| 565 | r[0], r[1], var, tmin, tmax, nticks) |
|---|
| 566 | if options is None: |
|---|
| 567 | cmd += ")" |
|---|
| 568 | else: |
|---|
| 569 | cmd += ", %s)"%options |
|---|
| 570 | self(cmd) |
|---|
| 571 | |
|---|
| 572 | def plot3d(self, *args): |
|---|
| 573 | r""" |
|---|
| 574 | Plot a 3d graph using Maxima / gnuplot. |
|---|
| 575 | |
|---|
| 576 | maxima.plot3d(f, '[x, xmin, xmax]', '[y, ymin, ymax]', '[grid, nx, ny]', options) |
|---|
| 577 | |
|---|
| 578 | INPUT: |
|---|
| 579 | f -- a string representing a function (such as f="sin(x)") |
|---|
| 580 | [var, min, max] |
|---|
| 581 | |
|---|
| 582 | EXAMPLES: |
|---|
| 583 | sage.: maxima.plot3d('1 + x^3 - y^2', '[x,-2,2]', '[y,-2,2]', '[grid,12,12]') |
|---|
| 584 | sage.: maxima.plot3d('sin(x)*cos(y)', '[x,-2,2]', '[y,-2,2]', '[grid,30,30]') |
|---|
| 585 | sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "sin-plot.eps"]' |
|---|
| 586 | sage.: maxima.plot3d('sin(x+y)', '[x,-5,5]', '[y,-1,1]', opts) |
|---|
| 587 | |
|---|
| 588 | The eps file is saved in the current working directory. |
|---|
| 589 | """ |
|---|
| 590 | self('plot3d(%s)'%(','.join([str(x) for x in args]))) |
|---|
| 591 | |
|---|
| 592 | def plot3d_parametric(self, r, vars, urange, vrange, options=None): |
|---|
| 593 | r""" |
|---|
| 594 | Plot a 3d parametric graph with r=(x,y,z), x = x(u,v), y = y(u,v), z = z(u,v), |
|---|
| 595 | for u = umin...umax, v = vmin...vmax using gnuplot with options. |
|---|
| 596 | |
|---|
| 597 | INPUT: |
|---|
| 598 | x, y, z -- a string representing a function (such as x="u^2+v^2", ...) |
|---|
| 599 | vars is a list or two strings representing variables (such as vars = ["u","v"]) |
|---|
| 600 | urange -- [umin, umax] |
|---|
| 601 | vrange -- [vmin, vmax] are lists of numbers with |
|---|
| 602 | umin < umax, vmin < vmax |
|---|
| 603 | options -- optional string representing plot2d options in gnuplot format |
|---|
| 604 | |
|---|
| 605 | OUTPUT: |
|---|
| 606 | displays a plot on screen or saves to a file |
|---|
| 607 | |
|---|
| 608 | EXAMPLES: |
|---|
| 609 | sage.: maxima.plot3d_parametric(["v*sin(u)","v*cos(u)","v"], ["u","v"],[-3.2,3.2],[0,3]) |
|---|
| 610 | sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "sin-cos-plot.eps"]' |
|---|
| 611 | sage.: maxima.plot3d_parametric(["v*sin(u)","v*cos(u)","v"], ["u","v"],[-3.2,3.2],[0,3],opts) |
|---|
| 612 | |
|---|
| 613 | The eps file is saved in the current working directory. |
|---|
| 614 | |
|---|
| 615 | Here is a torus: |
|---|
| 616 | |
|---|
| 617 | sage.: _ = maxima.eval("expr_1: cos(y)*(10.0+6*cos(x)); expr_2: sin(y)*(10.0+6*cos(x)); expr_3: -6*sin(x);") # optional |
|---|
| 618 | sage.: maxima.plot3d_parametric(["expr_1","expr_2","expr_3"], ["x","y"],[0,6],[0,6]) |
|---|
| 619 | |
|---|
| 620 | Here is a Mobius strip: |
|---|
| 621 | sage.: x = "cos(u)*(3 + v*cos(u/2))" |
|---|
| 622 | sage.: y = "sin(u)*(3 + v*cos(u/2))" |
|---|
| 623 | sage.: z = "v*sin(u/2)" |
|---|
| 624 | sage.: maxima.plot3d_parametric([x,y,z],["u","v"],[-3.1,3.2],[-1/10,1/10]) |
|---|
| 625 | """ |
|---|
| 626 | umin = urange[0] |
|---|
| 627 | umax = urange[1] |
|---|
| 628 | vmin = vrange[0] |
|---|
| 629 | vmax = vrange[1] |
|---|
| 630 | cmd = 'plot3d([%s, %s, %s], [%s, %s, %s], [%s, %s, %s]'%( |
|---|
| 631 | r[0], r[1], r[2], vars[0], umin, umax, vars[1], vmin, vmax) |
|---|
| 632 | if options is None: |
|---|
| 633 | cmd += ')' |
|---|
| 634 | else: |
|---|
| 635 | cmd += ', %s)'%options |
|---|
| 636 | maxima(cmd) |
|---|
| 637 | |
|---|
| 638 | def de_solve(maxima, de, vars, ics=None): |
|---|
| 639 | """ |
|---|
| 640 | Solves a 1st or 2nd order ordinary differential equation (ODE) |
|---|
| 641 | in two variables, possibly with initial conditions. |
|---|
| 642 | |
|---|
| 643 | INPUT: |
|---|
| 644 | de -- a string representing the ODE |
|---|
| 645 | vars -- a list of strings representing the two variables. |
|---|
| 646 | ics -- a triple of numbers [a,b1,b2] representing |
|---|
| 647 | y(a)=b1, y'(a)=b2 |
|---|
| 648 | |
|---|
| 649 | EXAMPLES: |
|---|
| 650 | sage.: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'], [1,1,1]) |
|---|
| 651 | y = 3*x - 2*%e^(x - 1) |
|---|
| 652 | sage.: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y']) |
|---|
| 653 | y = %k1*%e^x + %k2*%e^ - x + 3*x |
|---|
| 654 | sage.: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y']) |
|---|
| 655 | y = (%c - 3*( - x - 1)*%e^ - x)*%e^x |
|---|
| 656 | sage.: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'],[1,1]) |
|---|
| 657 | y = - %e^ - 1*(5*%e^x - 3*%e*x - 3*%e) |
|---|
| 658 | """ |
|---|
| 659 | if not isinstance(vars, str): |
|---|
| 660 | str_vars = '%s, %s'%(vars[1], vars[0]) |
|---|
| 661 | else: |
|---|
| 662 | str_vars = vars |
|---|
| 663 | maxima.eval('depends(%s)'%str_vars) |
|---|
| 664 | m = maxima(de) |
|---|
| 665 | a = 'ode2(%s, %s)'%(m.name(), str_vars) |
|---|
| 666 | if ics != None: |
|---|
| 667 | if len(ics) == 3: |
|---|
| 668 | cmd = "ic2("+a+",%s=%s,%s=%s,diff(%s,%s)=%s);"%(vars[0],ics[0], vars[1],ics[1], vars[1], vars[0], ics[2]) |
|---|
| 669 | return maxima(cmd) |
|---|
| 670 | if len(ics) == 2: |
|---|
| 671 | return maxima("ic1("+a+",%s=%s,%s=%s);"%(vars[0],ics[0], vars[1],ics[1])) |
|---|
| 672 | return maxima(a+";") |
|---|
| 673 | |
|---|
| 674 | def de_solve_laplace(self, de, vars, ics=None): |
|---|
| 675 | """ |
|---|
| 676 | Solves an ordinary differential equation (ODE) using Laplace transforms. |
|---|
| 677 | |
|---|
| 678 | INPUT: |
|---|
| 679 | de -- a string representing the ODE |
|---|
| 680 | (e.g., de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)") |
|---|
| 681 | vars -- a list of strings representing the variables |
|---|
| 682 | (e.g., vars = ["x","f"]) |
|---|
| 683 | ics -- a list of numbers representing initial conditions, |
|---|
| 684 | with symbols allowed which are represented by strings |
|---|
| 685 | (eg, f(0)=1, f'(0)=2 is ics = [0,1,2]) |
|---|
| 686 | |
|---|
| 687 | EXAMPLES: |
|---|
| 688 | sage.: maxima.clear('x'); maxima.clear('f') |
|---|
| 689 | sage.: maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)", ["x","f"], [0,1,2]) |
|---|
| 690 | f(x) = x*%e^x + %e^x |
|---|
| 691 | |
|---|
| 692 | sage.: maxima.clear('x'); maxima.clear('f') |
|---|
| 693 | sage.: f = maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)", ["x","f"]) |
|---|
| 694 | sage.: f |
|---|
| 695 | f(x) = x*%e^x*(at('diff(f(x),x,1),x = 0)) - f(0)*x*%e^x + f(0)*%e^x |
|---|
| 696 | sage.: f.display2d() |
|---|
| 697 | ! |
|---|
| 698 | x d ! x x |
|---|
| 699 | f(x) = x %e (-- (f(x))! ) - f(0) x %e + f(0) %e |
|---|
| 700 | dx ! |
|---|
| 701 | !x = 0 |
|---|
| 702 | |
|---|
| 703 | |
|---|
| 704 | \note{The second equation sets the values of $f(0)$ and |
|---|
| 705 | $f'(0)$ in maxima, so subsequent ODEs involving these |
|---|
| 706 | variables will have these initial conditions automatically |
|---|
| 707 | imposed.} |
|---|
| 708 | """ |
|---|
| 709 | if not (ics is None): |
|---|
| 710 | d = len(ics) |
|---|
| 711 | for i in range(0,d-1): |
|---|
| 712 | ic = 'atvalue(diff(%s(%s), %s, %s), %s = %s, %s)'%( |
|---|
| 713 | vars[1], vars[0], vars[0], i, vars[0], ics[0], ics[1+i]) |
|---|
| 714 | maxima.eval(ic) |
|---|
| 715 | return maxima('desolve(%s, %s(%s))'%(de, vars[1], vars[0])) |
|---|
| 716 | |
|---|
| 717 | ## def de_solve_laplace_plot(self, de,vars,ics,xrange,yrange,options=None): |
|---|
| 718 | ## """ |
|---|
| 719 | ## Plots the solution to an ODE using laplace transforms. |
|---|
| 720 | ## INPUT: de is a string representing the ODE |
|---|
| 721 | ## (eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)") |
|---|
| 722 | ## vars is a list of strings representing the variables |
|---|
| 723 | ## (eg, vars = ["x","f"]) |
|---|
| 724 | ## ics is a list of numbers representing initial conditions, |
|---|
| 725 | ## with symbols allowed which are represented by strings |
|---|
| 726 | ## (eg, f(0)=1, f'(0)=2 is ics = [0,1,2]) |
|---|
| 727 | |
|---|
| 728 | ## EXAMPLES: |
|---|
| 729 | ## sage: self.de_solve_laplace_plot("diff(f(x),x,2)=2*diff(f(x),x)-f(x)",["x","f"],[0,1,2],[-1,1],[-5,5]) |
|---|
| 730 | |
|---|
| 731 | ## Warning: The second equation sets the values of f(0) and f'(0) in maxima, so |
|---|
| 732 | ## subsequent ODEs involving these variables will have these initial conditions |
|---|
| 733 | ## automatically imposed. |
|---|
| 734 | ## """ |
|---|
| 735 | ## raise NotImplementedError |
|---|
| 736 | |
|---|
| 737 | ## def de_plot(self, de,vars,ic,xrange,yrange,options=None): |
|---|
| 738 | ## r""" |
|---|
| 739 | ## Plots solution to a 2nd order ODE. |
|---|
| 740 | |
|---|
| 741 | ## INPUT: |
|---|
| 742 | ## de is a string representing the ODE |
|---|
| 743 | ## (eg, de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)") |
|---|
| 744 | ## vars is a list or two strings representing variables (such as vars = ["x","y"]) |
|---|
| 745 | ## ics is a list of numbers representing initial conditions, |
|---|
| 746 | ## with symbols allowed which are represented by strings |
|---|
| 747 | ## (eg, y(0)=1, y'(0)=2 is ic = [0,1,2]) |
|---|
| 748 | ## xrange = [xmin, xmax], yrange = [ymin, ymax] are lists ofnumbers with xmin<xmax, ymin<ymax |
|---|
| 749 | ## options is an optional string representing plot2d options in gnuplot format |
|---|
| 750 | |
|---|
| 751 | ## EXAMPLES: |
|---|
| 752 | ## sage.: de = "diff(y,x,2) = 2*(1+x)" |
|---|
| 753 | ## sage.: de_plot(de,["x","y"],[1,2,3],[-4,4],[-10,10]) |
|---|
| 754 | ## sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "de_plot.eps"]' |
|---|
| 755 | ## sage.: de_plot(de,["x","y"],[1,2,3],[-4,4],[-10,10],opts) |
|---|
| 756 | |
|---|
| 757 | ## The eps file is saved in the current working directory. |
|---|
| 758 | ## """ |
|---|
| 759 | ## y = vars[1] |
|---|
| 760 | ## x = vars[0] |
|---|
| 761 | ## x0 = ic[0] |
|---|
| 762 | ## y0 = ic[1] |
|---|
| 763 | ## y1 = ic[2] |
|---|
| 764 | ## xmin = xrange[0] |
|---|
| 765 | ## xmax = xrange[1] |
|---|
| 766 | ## ymin = yrange[0] |
|---|
| 767 | ## ymax = yrange[1] |
|---|
| 768 | ## cmd1 = "(soln:ode2('"+de+","+y+","+x+"), tmp:IC2(soln,"+x+"="+str(x0)+","+y+"="+str(y0)+",'diff("+y+","+x+")="+str(y1)+"));" |
|---|
| 769 | ## #print cmd1 |
|---|
| 770 | ## print self(cmd1) |
|---|
| 771 | ## if options==None: |
|---|
| 772 | ## cmd2 = "plot2d(sublis(solve(tmp,"+y+"),"+y+"),["+x+","+str(xmin)+","+str(xmax)+"],["+y+","+str(ymin)+","+str(ymax)+"]);" |
|---|
| 773 | ## #print cmd2 |
|---|
| 774 | ## self(cmd2) |
|---|
| 775 | ## if options!=None: |
|---|
| 776 | ## cmd2 = "plot2d(sublis(solve(tmp,"+y+"),"+y+"),["+x+","+str(xmin)+","+str(xmax)+"],["+y+","+str(ymin)+","+str(ymax)+"],"+options+");" |
|---|
| 777 | ## #print cmd2 |
|---|
| 778 | ## self(cmd2) |
|---|
| 779 | |
|---|
| 780 | def solve_linear(self, eqns,vars): |
|---|
| 781 | """ |
|---|
| 782 | Wraps maxima's linsolve. |
|---|
| 783 | |
|---|
| 784 | INPUT: |
|---|
| 785 | eqns is a list of m strings, each rperesenting a linear question |
|---|
| 786 | in m <= n variables |
|---|
| 787 | vars is a list of n strings, each representing a variable |
|---|
| 788 | |
|---|
| 789 | EXAMPLES: |
|---|
| 790 | sage: eqns = ["x + z = y","2*a*x - y = 2*a^2","y - 2*z = 2"] |
|---|
| 791 | sage: vars = ["x","y","z"] |
|---|
| 792 | sage: maxima.solve_linear(eqns, vars) |
|---|
| 793 | [x = a + 1,y = 2*a,z = a - 1] |
|---|
| 794 | """ |
|---|
| 795 | eqs = "[" |
|---|
| 796 | for i in range(len(eqns)): |
|---|
| 797 | if i<len(eqns)-1: |
|---|
| 798 | eqs = eqs + eqns[i]+"," |
|---|
| 799 | if i==len(eqns)-1: |
|---|
| 800 | eqs = eqs + eqns[i]+"]" |
|---|
| 801 | vrs = "[" |
|---|
| 802 | for i in range(len(vars)): |
|---|
| 803 | if i<len(vars)-1: |
|---|
| 804 | vrs = vrs + vars[i]+"," |
|---|
| 805 | if i==len(vars)-1: |
|---|
| 806 | vrs = vrs + vars[i]+"]" |
|---|
| 807 | return self('linsolve(%s, %s)'%(eqs, vrs)) |
|---|
| 808 | |
|---|
| 809 | def unit_quadratic_integer(self, n): |
|---|
| 810 | r""" |
|---|
| 811 | Finds a unit of the ring of integers of the quadratic number |
|---|
| 812 | field $\Q(\sqrt{n})$, $n>1$, using the qunit maxima command. |
|---|
| 813 | |
|---|
| 814 | EXAMPLE: |
|---|
| 815 | sage: u = maxima.unit_quadratic_integer(101) |
|---|
| 816 | sage: u.parent() |
|---|
| 817 | Number Field in a with defining polynomial x^2 - 101 |
|---|
| 818 | sage: u |
|---|
| 819 | a + 10 |
|---|
| 820 | sage: u = maxima.unit_quadratic_integer(13) |
|---|
| 821 | sage: u |
|---|
| 822 | 5*a + 18 |
|---|
| 823 | sage: u.parent() |
|---|
| 824 | Number Field in a with defining polynomial x^2 - 13 |
|---|
| 825 | """ |
|---|
| 826 | from sage.rings.all import QuadraticField, Integer |
|---|
| 827 | # Take square-free part so sqrt(n) doesn't get simplified further by maxima |
|---|
| 828 | # (The original version of this function would yield wrong answers if |
|---|
| 829 | # n is not squarefree.) |
|---|
| 830 | n = Integer(n).square_free_part() |
|---|
| 831 | if n < 1: |
|---|
| 832 | raise ValueError, "n (=%s) must be >= 1"%n |
|---|
| 833 | s = str(self('qunit(%s)'%n)).lower() |
|---|
| 834 | r = re.compile('sqrt\(.*\)') |
|---|
| 835 | s = r.sub('a', s) |
|---|
| 836 | a = QuadraticField(n, 'a').gen() |
|---|
| 837 | return eval(s) |
|---|
| 838 | |
|---|
| 839 | def plot_list(self, ptsx, ptsy, options=None): |
|---|
| 840 | r""" |
|---|
| 841 | Plots a curve determined by a sequence of points. |
|---|
| 842 | |
|---|
| 843 | INPUT: |
|---|
| 844 | ptsx -- [x1,...,xn], where the xi and yi are real, |
|---|
| 845 | ptsy -- [y1,...,yn] |
|---|
| 846 | options -- a string representing maxima plot2d options. |
|---|
| 847 | |
|---|
| 848 | The points are (x1,y1), (x2,y2), etc. |
|---|
| 849 | |
|---|
| 850 | This function requires maxima 5.9.2 or newer. |
|---|
| 851 | |
|---|
| 852 | \note{More that 150 points can sometimes lead to the program |
|---|
| 853 | hanging. Why?} |
|---|
| 854 | |
|---|
| 855 | EXAMPLES: |
|---|
| 856 | sage.: zeta_ptsx = [ (pari(1/2 + i*I/10).zeta().real()).precision(1) for i in range (70,150)] |
|---|
| 857 | sage.: zeta_ptsy = [ (pari(1/2 + i*I/10).zeta().imag()).precision(1) for i in range (70,150)] |
|---|
| 858 | sage.: maxima.plot_list(zeta_ptsx, zeta_ptsy) |
|---|
| 859 | sage.: opts='[gnuplot_preamble, "set nokey"], [gnuplot_term, ps], [gnuplot_out_file, "zeta.eps"]' |
|---|
| 860 | sage.: maxima.plot_list(zeta_ptsx, zeta_ptsy, opts) |
|---|
| 861 | """ |
|---|
| 862 | cmd = 'plot2d([discrete,%s, %s]'%(ptsx, ptsy) |
|---|
| 863 | if options is None: |
|---|
| 864 | cmd += ')' |
|---|
| 865 | else: |
|---|
| 866 | cmd += ', %s)'%options |
|---|
| 867 | self(cmd) |
|---|
| 868 | |
|---|
| 869 | |
|---|
| 870 | def plot_multilist(self, pts_list, options=None): |
|---|
| 871 | r""" |
|---|
| 872 | Plots a list of list of points pts_list=[pts1,pts2,...,ptsn], |
|---|
| 873 | where each ptsi is of the form [[x1,y1],...,[xn,yn]] |
|---|
| 874 | x's must be integers and y's reals |
|---|
| 875 | options is a string representing maxima plot2d options. |
|---|
| 876 | |
|---|
| 877 | Requires maxima 5.9.2 at least. |
|---|
| 878 | \note{More that 150 points can sometimes lead to the |
|---|
| 879 | program hanging.} |
|---|
| 880 | |
|---|
| 881 | EXAMPLES: |
|---|
| 882 | sage.: xx = [ i/10.0 for i in range (-10,10)] |
|---|
| 883 | sage.: yy = [ i/10.0 for i in range (-10,10)] |
|---|
| 884 | sage.: x0 = [ 0 for i in range (-10,10)] |
|---|
| 885 | sage.: y0 = [ 0 for i in range (-10,10)] |
|---|
| 886 | sage.: zeta_ptsx1 = [ (pari(1/2+i*I/10).zeta().real()).precision(1) for i in range (10)] |
|---|
| 887 | sage.: zeta_ptsy1 = [ (pari(1/2+i*I/10).zeta().imag()).precision(1) for i in range (10)] |
|---|
| 888 | sage.: maxima.plot_multilist([[zeta_ptsx1,zeta_ptsy1],[xx,y0],[x0,yy]]) |
|---|
| 889 | sage.: zeta_ptsx1 = [ (pari(1/2+i*I/10).zeta().real()).precision(1) for i in range (10,150)] |
|---|
| 890 | sage.: zeta_ptsy1 = [ (pari(1/2+i*I/10).zeta().imag()).precision(1) for i in range (10,150)] |
|---|
| 891 | sage.: maxima.plot_multilist([[zeta_ptsx1,zeta_ptsy1],[xx,y0],[x0,yy]]) |
|---|
| 892 | sage.: opts='[gnuplot_preamble, "set nokey"]' |
|---|
| 893 | sage.: maxima.plot_multilist([[zeta_ptsx1,zeta_ptsy1],[xx,y0],[x0,yy]],opts) |
|---|
| 894 | """ |
|---|
| 895 | n = len(pts_list) |
|---|
| 896 | cmd = '[' |
|---|
| 897 | for i in range(n): |
|---|
| 898 | if i < n-1: |
|---|
| 899 | cmd = cmd+'[discrete,'+str(pts_list[i][0])+','+str(pts_list[i][1])+'],' |
|---|
| 900 | if i==n-1: |
|---|
| 901 | cmd = cmd+'[discrete,'+str(pts_list[i][0])+','+str(pts_list[i][1])+']]' |
|---|
| 902 | #print cmd |
|---|
| 903 | if options is None: |
|---|
| 904 | self('plot2d('+cmd+')') |
|---|
| 905 | else: |
|---|
| 906 | self('plot2d('+cmd+','+options+')') |
|---|
| 907 | |
|---|
| 908 | |
|---|
| 909 | class MaximaElement(ExpectElement): |
|---|
| 910 | def numer(self): |
|---|
| 911 | return self.comma('numer') |
|---|
| 912 | |
|---|
| 913 | def real(self): |
|---|
| 914 | return self.realpart() |
|---|
| 915 | |
|---|
| 916 | def imag(self): |
|---|
| 917 | return self.imagpart() |
|---|
| 918 | |
|---|
| 919 | def str(self): |
|---|
| 920 | self._check_valid() |
|---|
| 921 | P = self.parent() |
|---|
| 922 | return P.get(self._name) |
|---|
| 923 | |
|---|
| 924 | def __repr__(self): |
|---|
| 925 | self._check_valid() |
|---|
| 926 | P = self.parent() |
|---|
| 927 | if P._display2d: |
|---|
| 928 | return self.display2d(onscreen=False) |
|---|
| 929 | else: |
|---|
| 930 | return P.get(self._name) |
|---|
| 931 | |
|---|
| 932 | def display2d(self, onscreen=True): |
|---|
| 933 | """ |
|---|
| 934 | EXAMPLES: |
|---|
| 935 | sage: F = maxima('x^5 - y^5').factor() |
|---|
| 936 | sage: F.display2d () |
|---|
| 937 | 4 3 2 2 3 4 |
|---|
| 938 | - (y - x) (y + x y + x y + x y + x ) |
|---|
| 939 | """ |
|---|
| 940 | self._check_valid() |
|---|
| 941 | P = self.parent() |
|---|
| 942 | s = P._eval_line('display2d : true; %s'%self.name(), reformat=False) |
|---|
| 943 | P._eval_line('display2d : false', reformat=False) |
|---|
| 944 | i = s.find('true') |
|---|
| 945 | i += s[i:].find('\n') |
|---|
| 946 | #j = s.rfind('(%o') |
|---|
| 947 | #s = s[:j] |
|---|
| 948 | j = s.rfind('(%o') |
|---|
| 949 | s = s[i:j-2] |
|---|
| 950 | i = s.find('(%o') |
|---|
| 951 | j = i + s[i:].find(')') |
|---|
| 952 | s = s[:i] + ' '*(j-i+1) + s[j+1:] |
|---|
| 953 | s = s.lstrip('\n') |
|---|
| 954 | if onscreen: |
|---|
| 955 | print s |
|---|
| 956 | else: |
|---|
| 957 | return s |
|---|
| 958 | |
|---|
| 959 | def diff(self, var='x', n=1): |
|---|
| 960 | """ |
|---|
| 961 | Return the n-th derivative of self. |
|---|
| 962 | |
|---|
| 963 | INPUT: |
|---|
| 964 | var -- variable (default: 'x') |
|---|
| 965 | n -- integer (default: 1) |
|---|
| 966 | |
|---|
| 967 | OUTPUT: |
|---|
| 968 | n-th derivative of self with respect to the variable var |
|---|
| 969 | |
|---|
| 970 | EXAMPLES: |
|---|
| 971 | sage: f = maxima('x^2') |
|---|
| 972 | sage: f.diff() |
|---|
| 973 | 2*x |
|---|
| 974 | sage: f.diff('x') |
|---|
| 975 | 2*x |
|---|
| 976 | sage: f.diff('x', 2) |
|---|
| 977 | 2 |
|---|
| 978 | sage: maxima('sin(x^2)').diff('x',4) |
|---|
| 979 | 16*x^4*sin(x^2) - 12*sin(x^2) - 48*x^2*cos(x^2) |
|---|
| 980 | |
|---|
| 981 | sage: f = maxima('x^2 + 17*y^2') |
|---|
| 982 | sage: f.diff('x') |
|---|
| 983 | 2*x |
|---|
| 984 | sage: f.diff('y') |
|---|
| 985 | 34*y |
|---|
| 986 | """ |
|---|
| 987 | return ExpectElement.__getattr__(self, 'diff')(var, n) |
|---|
| 988 | |
|---|
| 989 | derivative = diff |
|---|
| 990 | |
|---|
| 991 | def integral(self, var='x', min=None, max=None): |
|---|
| 992 | r""" |
|---|
| 993 | Return the integral of self with respect to the variable x. |
|---|
| 994 | |
|---|
| 995 | INPUT: |
|---|
| 996 | var -- variable |
|---|
| 997 | min -- default: None |
|---|
| 998 | max -- default: None |
|---|
| 999 | |
|---|
| 1000 | Returns the definite integral if xmin is not None, |
|---|
| 1001 | otherwise returns an indefinite integral. |
|---|
| 1002 | |
|---|
| 1003 | EXAMPLES: |
|---|
| 1004 | sage: maxima('x^2+1').integral() |
|---|
| 1005 | x^3/3 + x |
|---|
| 1006 | sage: maxima('x^2+ 1 + y^2').integral('y') |
|---|
| 1007 | y^3/3 + x^2*y + y |
|---|
| 1008 | sage: maxima('x / (x^2+1)').integral() |
|---|
| 1009 | log(x^2 + 1)/2 |
|---|
| 1010 | sage: maxima('1/(x^2+1)').integral() |
|---|
| 1011 | atan(x) |
|---|
| 1012 | sage.: maxima('1/(x^2+1)').integral('x', 0, infinity) |
|---|
| 1013 | %pi/2 |
|---|
| 1014 | sage: maxima('x/(x^2+1)').integral('x', -1, 1) |
|---|
| 1015 | 0 |
|---|
| 1016 | |
|---|
| 1017 | sage: f = maxima('exp(x^2)').integral('x',0,1); f |
|---|
| 1018 | -sqrt(%pi)*%i*erf(%i)/2 |
|---|
| 1019 | sage: f.numer() # I wonder how to get a real number (~1.463)?? |
|---|
| 1020 | -.8862269254527579*%i*erf(%i) |
|---|
| 1021 | """ |
|---|
| 1022 | I = ExpectElement.__getattr__(self, 'integrate') |
|---|
| 1023 | if min is None: |
|---|
| 1024 | return I(var) |
|---|
| 1025 | else: |
|---|
| 1026 | if max is None: |
|---|
| 1027 | raise ValueError, "neither or both of min/max must be specified." |
|---|
| 1028 | return I(var, min, max) |
|---|
| 1029 | |
|---|
| 1030 | integrate = integral |
|---|
| 1031 | |
|---|
| 1032 | |
|---|
| 1033 | |
|---|
| 1034 | |
|---|
| 1035 | def __float__(self): |
|---|
| 1036 | return float(str(self.numer())) |
|---|
| 1037 | |
|---|
| 1038 | def __len__(self): |
|---|
| 1039 | """ |
|---|
| 1040 | Return the length of a list. |
|---|
| 1041 | |
|---|
| 1042 | EXAMPLES: |
|---|
| 1043 | sage: v = maxima('create_list(x^i,i,0,5)') |
|---|
| 1044 | sage: len(v) |
|---|
| 1045 | 6 |
|---|
| 1046 | """ |
|---|
| 1047 | self._check_valid() |
|---|
| 1048 | return int(self.parent().eval('length(%s)'%self.name())) |
|---|
| 1049 | |
|---|
| 1050 | def __getitem__(self, n): |
|---|
| 1051 | r""" |
|---|
| 1052 | Return the n-th element of this list. |
|---|
| 1053 | |
|---|
| 1054 | \note{Lists are 0-based when accessed via the \sage interface, |
|---|
| 1055 | not 1-based as they are in the Maxima interpreter.} |
|---|
| 1056 | |
|---|
| 1057 | EXAMPLES: |
|---|
| 1058 | sage: v = maxima('create_list(i*x^i,i,0,5)'); v |
|---|
| 1059 | [0,x,2*x^2,3*x^3,4*x^4,5*x^5] |
|---|
| 1060 | sage: v[3] |
|---|
| 1061 | 3*x^3 |
|---|
| 1062 | sage: v[0] |
|---|
| 1063 | 0 |
|---|
| 1064 | sage: v[10] |
|---|
| 1065 | Traceback (most recent call last): |
|---|
| 1066 | ... |
|---|
| 1067 | IndexError: n = (10) must be between 0 and 5 |
|---|
| 1068 | """ |
|---|
| 1069 | n = int(n) |
|---|
| 1070 | if n < 0 or n >= len(self): |
|---|
| 1071 | raise IndexError, "n = (%s) must be between %s and %s"%(n, 0, len(self)-1) |
|---|
| 1072 | return ExpectElement.__getitem__(self, n+1) |
|---|
| 1073 | |
|---|
| 1074 | def subst(self, val): |
|---|
| 1075 | return self.comma(val) |
|---|
| 1076 | |
|---|
| 1077 | def comma(self, args): |
|---|
| 1078 | self._check_valid() |
|---|
| 1079 | P = self.parent() |
|---|
| 1080 | return P('%s, %s'%(self.name(), args)) |
|---|
| 1081 | |
|---|
| 1082 | def _latex_(self): |
|---|
| 1083 | self._check_valid() |
|---|
| 1084 | P = self.parent() |
|---|
| 1085 | s = maxima.eval('tex(%s)'%self.name()) |
|---|
| 1086 | s = s[2:-7] |
|---|
| 1087 | # Actually trying the latex on some examples |
|---|
| 1088 | # quickly reveals serious bugs in it. The |
|---|
| 1089 | # following are some attempts to program around |
|---|
| 1090 | # these. |
|---|
| 1091 | s = s.replace('\\sin', '\\sin{}') |
|---|
| 1092 | s = s.replace('\\cos', '\\cos{}') |
|---|
| 1093 | s = s.replace('\\tan', '\\tan{}') |
|---|
| 1094 | s = s.replace('\\arcsin', '\\sin^{-1}{}') |
|---|
| 1095 | s = s.replace('\\arccos', '\\cos^{-1}{}') |
|---|
| 1096 | s = s.replace('\\arctan', '\\tan^{-1}{}') |
|---|
| 1097 | # TODO: What to do about this, which won't work!? |
|---|
| 1098 | #s = s.replace('\\sinh', '\\sinh{}') |
|---|
| 1099 | #s = s.replace('\\cosh', '\\cosh{}') |
|---|
| 1100 | #s = s.replace('\\tanh', '\\tanh{}') |
|---|
| 1101 | return s |
|---|
| 1102 | |
|---|
| 1103 | def _matrix_(self, R): |
|---|
| 1104 | r""" |
|---|
| 1105 | If self is a Maxima matrix, return the corresponding \sage |
|---|
| 1106 | matrix over the \sage ring $R$. |
|---|
| 1107 | |
|---|
| 1108 | This may or may not work depending in how complicated the |
|---|
| 1109 | entries of self are! It only works if the entries of self |
|---|
| 1110 | can be coerced as strings to produce meaningful elements |
|---|
| 1111 | of $R$. |
|---|
| 1112 | |
|---|
| 1113 | EXAMPLES: |
|---|
| 1114 | sage: _ = maxima.eval("f[i,j] := i/j") |
|---|
| 1115 | sage: A = maxima('genmatrix(f,4,4)'); A |
|---|
| 1116 | matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1]) |
|---|
| 1117 | sage: A._matrix_(QQ) |
|---|
| 1118 | [ 1 1/2 1/3 1/4] |
|---|
| 1119 | [ 2 1 2/3 1/2] |
|---|
| 1120 | [ 3 3/2 1 3/4] |
|---|
| 1121 | [ 4 2 4/3 1] |
|---|
| 1122 | |
|---|
| 1123 | You can also use the \code{matrix} command (which is defined |
|---|
| 1124 | in \code{sage.misc.functional}): |
|---|
| 1125 | sage: matrix(A, QQ) |
|---|
| 1126 | [ 1 1/2 1/3 1/4] |
|---|
| 1127 | [ 2 1 2/3 1/2] |
|---|
| 1128 | [ 3 3/2 1 3/4] |
|---|
| 1129 | [ 4 2 4/3 1] |
|---|
| 1130 | """ |
|---|
| 1131 | from sage.matrix.all import MatrixSpace |
|---|
| 1132 | self._check_valid() |
|---|
| 1133 | P = self.parent() |
|---|
| 1134 | nrows = int(P.eval('length(%s)'%self.name())) |
|---|
| 1135 | if nrows == 0: |
|---|
| 1136 | return MatrixSpace(R, 0, 0)(0) |
|---|
| 1137 | ncols = int(P.eval('length(%s[1])'%self.name())) |
|---|
| 1138 | M = MatrixSpace(R, nrows, ncols) |
|---|
| 1139 | s = self.str().replace('matrix','').replace(',',"','").\ |
|---|
| 1140 | replace("]','[","','").replace('([',"['").replace('])',"']") |
|---|
| 1141 | s = eval(s) |
|---|
| 1142 | return M([R(x) for x in s]) |
|---|
| 1143 | |
|---|
| 1144 | def partial_fraction_decomposition(self, var='x'): |
|---|
| 1145 | """ |
|---|
| 1146 | Return the partial fraction decomposition of self with respect to |
|---|
| 1147 | the variable var. |
|---|
| 1148 | |
|---|
| 1149 | EXAMPLES: |
|---|
| 1150 | sage: f = maxima('1/((1+x)*(x-1))') |
|---|
| 1151 | sage: f.partial_fraction_decomposition('x') |
|---|
| 1152 | 1/(2*(x - 1)) - 1/(2*(x + 1)) |
|---|
| 1153 | sage: f.partial_fraction_decomposition('x').display2d() |
|---|
| 1154 | 1 1 |
|---|
| 1155 | --------- - --------- |
|---|
| 1156 | 2 (x - 1) 2 (x + 1) |
|---|
| 1157 | """ |
|---|
| 1158 | return self.partfrac(var) |
|---|
| 1159 | |
|---|
| 1160 | |
|---|
| 1161 | def is_MaximaElement(x): |
|---|
| 1162 | return isinstance(x, MaximaElement) |
|---|
| 1163 | |
|---|
| 1164 | # An instance |
|---|
| 1165 | maxima = Maxima(script_subdirectory=None) |
|---|
| 1166 | |
|---|
| 1167 | def reduce_load_Maxima(): |
|---|
| 1168 | return maxima |
|---|
| 1169 | |
|---|
| 1170 | import os |
|---|
| 1171 | def maxima_console(): |
|---|
| 1172 | os.system('maxima') |
|---|
| 1173 | |
|---|
| 1174 | def maxima_version(): |
|---|
| 1175 | return os.popen('maxima --version').read().split()[1] |
|---|
| 1176 | |
|---|
| 1177 | def __doctest_cleanup(): |
|---|
| 1178 | import sage.interfaces.quit |
|---|
| 1179 | sage.interfaces.quit.expect_quitall() |
|---|
