On Some Unsolved Problems of Mathematica

In this paper, some unsolved problems of the Mathematica software package are documented. At first, it is shown, using a number of examples, that the processing (simplification) of rational-fractional expressions involving powers in the general form, has been implemented in Maple more carefully than in Mathematica. Then, an error in Mathematica is demonstrated, leading to incorrect results at a change of a function’s body. For each example of symbolic or symbolicnumeric computations, alternate routes for solving the emerged problem are proposed (where possible). An added problem, related to Mathematica’s computing speed is documented, employing examples related to two-dimensional gas dynamics problems. It is shown that Mathematica computes rather slowly (about one thousand times slower) compared to a Fortran code. Keywords-Mathematica; analytic; numerical; computations; simplification; expressions; Fortran


INTRODUCTION
Modern theoretical mathematics and the development of information technologies are strongly connected to powerful software packages, widely applied in the areas of science and engineering, such as Mathematica [1], which is one of the most powerful tools supporting scientific research, applied design, and teaching of specialists, postgraduates and students.Mathematica was designed as software mainly aiming to automate the work of scientists and mathematicians.The software contains a wide spectrum of mathematical methods and numerical algorithms for conducting theoretical and applied research as well as a powerful system for the visualization of obtained results.Starting from version 2, the subsequent versions of Mathematica are the world leaders in the area of symbolic mathematics for PCs, providing the execution of numerical computations, employing elegant graphical output forms, and also the execution of especially laborious analytic computations and transformations.Mathematica has been adapted to the most recent operating systems and hardware, making it available to most platforms including Microsoft Windows, Apple Mac and Linux.
Mathematica enables the performance of laborious analytic computations.However, the possibility of errors in the produced results due to internal software design is a rather important issue, considering that the source of such problems is rarely detected by the user.Even after the error has been localized, the end-user will need further time in order to replace the piece of his code where some built-in operators or Mathematica functions work incorrectly.Even though, it is usually possible to implement such alternate approaches, due to the large number of available functions, such internal errors and/or imperfections of Mathematica reduce its efficiency.Most importantly, in some cases, the erroneous results produced may be externally believable and thus not identified as such.Therefore, when there are no verifying mathematical relations, it is necessary to duplicate the analytic computations by using a different software system, for example, Maple, REDUCE, etc.
The capabilities of seven general-purpose systems for analytic computations, namely Axiom, Derive, Macsyma, Maple, Mathematica, MuPAD, and REDUCE were analyzed previously [2] on 542 short problems.Mathematica produced completely wrong results in 13 of the 542 examples.Some added examples, not included in the above set, are documented in the present paper.In these cases, Mathematica either outputs an incorrect result or outputs a correct result having a bulky form, and one fails in its simplification with the aid of the builtin functions of the system.For each such example, alternate solving routes are proposed where possible.

FUNCTIONS
Consider the basic problem of mathematical analysis of the calculation of the sum:

www.etasr.com Vorozhtsov: On Some Unsolved Problems of Mathematica
The Mathematica command:

S=Sum[j*(j-1)*x^(j-2),{j,2,n}],
where Sum[…] is the built-in function of the system, gives the following solution: ( 1 ) It is clear that this expression can be simplified by dividing its numerator and denominator by x.Instead of doing this by hand using the pen and paper one can try to use the built-in function Simplify[…], which is intended for the simplification of expressions.However, instead of the expected simplification of (1), the function produces the same result (1).Another built-in function, PowerExpand[…], enables the simplification of expressions involving power functions.But the application of this function to (1) also leaves this expression without change.The Mathematica function Expand[S] expands a complex fraction S into the sum of elementary fractions: It can be seen from ( 2) that Mathematica has performed in each elementary fraction the divisions of both the numerators and denominators by the factor x. Therefore, one can now hope that the subsequent application of the function Simplify[…] to the sum of fractions (2) would yield the desired result.However, at the simplification of (2), Mathematica has multiplied both the numerator and denominator of the final expression by x, with the output being again the original expression (1).Mathematica also contains the function Together[…] which reverts the sum of fractions to a single fraction.Since there is no factor x in the denominators of (2), one may hope that this function will produce the result in the desired simplified form.However, the application of this function to (2) leads again to expression (1).
It should be noted that the Maple software produces the desired result for this task, using its simplify command, as shown below: S:=simplify(sum(j*(j-1)*x^(j-2),j=2..n)); . ) 1 ( The problem of the ideal incompressible fluid flow around a motionless sphere is a relatively simple example of a threedimensional fluid dynamics problem.The fluid flow is assumed to be directed at infinity along the z axis of a , respectively, and let p be the dimensional pressure of the fluid.The nondimensionalization of independent and dependent variables in the problem under study is given by: where R is the dimensional sphere radius, U ∞ is the freestream velocity and ρ is the fluid density.It is assumed for definiteness that the sphere center lies in the coordinate origin (0,0,0) of the coordinate system Oxyz.The exact solution of the problem under consideration is then obtained with the aid of the velocity potential theory so that, for example, the velocity component u has the form: The following notations are used here: fi is the velocity potential, V0 = U ∞ , x1, y1 and z1 are the dimensionless spatial coordinates and u1a is the dimensionless velocity component along the x-axis.The following expression is obtained for u1a as a result: It is shown from (4) that the function PowerExpand[u1/V0] has "missed" the presence of the common factor R 2•5/2 = R 5 in the denominator.
And, finally, It is seen from ( 5) that the combination of the functions Simplify[…] and PowerExpand[…] has led to the correct result.
It is interesting to solve the same task with Maple and to see how many steps are needed to obtain ( 5 Let now z = cos 2 x + i tg(x)/x, where x is a real variable.The application of function Re[z] in this case gives the answer Re(z) = Re(cos 2 (x)) -Im(tg(x)/x) because it is assumed in Mathematica that each of the functions cos 2 (x) and tg(x)/x may be complex.The function ComplexExpand[z] performs the simplification of the complex number z = a + ib under the assumption that both expressions a and b are real. Example: The answer is obtained in the form Re(z) = cos 2 x.However, the operation zim1 = ComplexExpand[Im[z]] leads to: instead of the expected: Using the well-known trigonometric relations it is easy to transform (6) to (7).The right-hand side of (6), however, requires 20 symbols, and the right-hand side of (7) requires 7 symbols, if one counts the horizontal bar separating the numerator from the denominator as a single symbol.The builtin function ByteCount[expr] available in Mathematica counts the number of bytes used internally by the system to store the expression expr.The call ByteCount[zim1] gives for the right-hand side of (6) the answer 488 bytes.Let us now try to simplify expression (6) with the aid of the function Simplify[…]: the command zim2 = Simplify[zim1] gives (7).ByteCount[zim2] gives 160 bytes.This means that the internal representation of (6) in Mathematica is three times longer compared to (7).
The danger of the appearance of bulky expressions like (6) in the process of symbolic manipulations consists of the fact that if (6) is used at the next stages of the analytic procedure then this may lead to a multiple "swell" of the final expressions.Further, bulky expressions deteriorate the visibility and understanding of final results.It is, therefore, desirable to execute the multi-stage analytic procedures in the interactive regime.As soon as a bulky expression of type (6) suddenly appears at a stage, one can try to simplify it with the aid of the following functions: Let us now see how the Maple system solves the task of finding the imaginary part of the complex number z = cos 2 x + i tg(x)/x.The function evalc(…)in the Maple system is an analog of the function ComplexExpand[…] of Mathematica.
Let us return to Mathematica.In a search for a simple alternative way of obtaining Im(z) in the optimal form (7) one can use the following definition of the imaginary part of the number z: Im(z) is the coefficient affecting i.

www.etasr.com Vorozhtsov: On Some Unsolved Problems of Mathematica
To find the coefficients affecting some variable, the function Coefficient[…] is provided.
Let us use it: In this case, Mathematica prints the following diagnostic message: "Coefficient:: ivar : i is not a valid variable.>>" and does not output any result.One can circumvent this situation by replacing I in (8) with a different letter, say j: Coefficient[z/.I -> j, j].
The answer is then obtained in the optimal form (7).
Summarizing the contents of this and foregoing sections one can draw the following conclusion: the simplification of analytic expressions has been implemented more carefully in Maple than in Mathematica.

IV. AN EXAMPLE OF INCORRECT WORK OF MATHEMATICA AT A CHANGE OF A FUNCTION'S BODY
As a preliminary discussion, we assume that there is a variable parameter M in the user's Mathematica code, which must be specified at the beginning of the code.To conserve the information about already computed variants it is important to store in the program the information about already considered values of the parameter M.This is also convenient from the viewpoint that it will be easy to return again to one of the computed runs, which proved to be the best among all other variants.It is very simple to store the information about already used values of M in the following way: M = 5; M = 10; M = 15; M = 20.Here the last assignment operator cancels any of the foregoing assignment operators of the form M = …, so that one can be sure that the Mathematica program of the user will execute the next run with the value M = 20.This is always the case in practice.This time, Mathematica has produced the correct result: expr2 = (y 2 -4z) 2 + 2z.

Let us now consider the following
V. SLOW NUMERICAL COMPUTATIONS Many engineering problems are described by systems of partial differential equations.Common examples include problems of fluid dynamics, magnetohydrodynamics, elasticity theory etc. Problems whose solutions depend on two or three spatial variables and the time variable t are considered the most laborious computationally.Such problems are commonly computed with the use of advanced supercomputers and parallel computations.
The ability to perform an on-line analysis of the results obtained through the simulation of multidimensional problems increases significantly with the use of graphical representations.ANSYS Fluent (www.ansys.com) is a well known software that integrates not only the programs for the numerical solution of gas dynamics problems, structural mechanics, and electromagnetics but also a number of computer graphics procedures for the visualization of the results of multidimensional computations.
Since the Mathematica software had already integrated initially the means for symbolic manipulations, numerical computations and graphical visualization of the results with the aid of numerous built-in functions, it appeared to be attractive for the specialists concerned with numerical modeling of various applied problems.Examples of the numerical solutions of some relatively simple one-, two-, and three-dimensional problems of continuum mechanics were presented in [3], which may be computed with acceptable CPU time expenses within the framework of Mathematica at a fairly moderate total number of spatial computing mesh nodes.An attempt to apply Mathematica for the numerical solution of two-dimensional problems of gas dynamics on curvilinear spatial computing meshes was implemented successfully in [4].But it turned out that the developed Mathematica code had a computational speed, which was many times smaller than the speed of the corresponding Fortran code.The developers of the considered system recommended the use of the built-in function Compile[…] to accelerate the numerical computations in Mathematica.However, its application reduced the CPU time needed for solving the two-dimensional problems of gas dynamics only by the factor of 5/3 = 1.67 [4] rather than by the factor of 20, as it was declared in [1].
The Euler equations, which represent one of the basic mathematical models for gas dynamics problems, have the following form for the two-dimensional flow of an inviscid, non-heat-conducting compressible gas: ∂w/∂t + ∂f(w)/∂x + ∂g(w)/∂y = 0, ( where x and y are the Cartesian spatial coordinates, t is the time, and Here p,ρ,u, and v denote the pressure, density, and velocity components; , ε is the specific internal energy, H = E + (p/ρ).The superscript T denotes the transposition operation.The ideal gas equation of state p = (γ -1)ρε, γ = const > 1 is used.An explicit TVD method as applied to equations ( 9) and (10) was described in detail in [5], and therefore, the computational formulas of this method are not presented here.The method has either the second or third order of accuracy in spatial variables and in the subregions of smooth flow depending on the numerical value of the available weight parameter φ.
The method uses the flux limiters as well as the Roe averaging.As in other TVD methods, a switching from the higher approximation order to the first approximation order is carried out in the regions of shock waves and contact discontinuities for the purpose of obtaining monotonous profiles of numerical solution [5].
The problem of the oblique shock reflection from a wall is frequently used for verification of new numerical methods for solving (9) and (10).The exact solution of this test problem represents a piecewise constant function and is found numerically with machine accuracy with the aid of the theory of the oblique shock waves [3].The value φ = π/6 for the angle between the incident shock and the x axis was used in the computational examples presented below, as shown in Figure 1.
Numerical solution was found by the pseudo-unsteady method.At the initial moment of time, the entire flow field was specified equal to the values of the undisturbed supersonic free stream, which were specified for subregion 1, as shown in Figure 1, that is the initial gas flow was parallel with the xaxis.
The criterion for the difference solution w n convergence to the limiting stationary solution was taken in the form Res(n) < ε, where ε is the user-specified small positive number; where , ˆ± ± are the difference approximations of fluxes f(w) and g(w) in ( 9) by the method of [5], h 1 and h 2 are the steps of the uniform rectangular computational grid in the (x,y) plane.In both runs -on the 160x40 grid and on the 80x20 grid -the information about the solution residual was printed every 50 time steps.On the 160x40 grid, 2600 time steps were executed, and on the 80x20 grid, a much smaller number of 1600 time steps were done.
One can assume with an insignificant error that the twodimensional gas dynamics problems are solved in Mathematica 4.1.0one thousand times slower compared to Fortran.In the case of Mathematica 8.0.4.0, the corresponding slow-down factor is higher and is equal to about 1270.Such unfavorable slow-down factors make the considered versions of Mathematica inapplicable for numerical solution of complex multidimensional problems of fluid dynamics.
It should be noted that Mathematica has been written using the C language.This language was developed at later times compared to the Fortran language.Therefore, the C compilers have absorbed the later developments in the area of the programming theory.In this connection, it is reasonable to expect that C compilers should generate exe files which should Mathematica has a large number of the built-in functions of computer graphics, which enable the obtaining of highquality graphs of curves, pictures of surfaces, contour plots, etc.In order to use these possibilities for the purpose of a rapid on-line visualization of the results of numerical solution of complex multidimensional problems of fluid dynamics and mathematical physics the present author has implemented a simple interface between Fortran and Mathematica.The realization of this interface reduces to the following: upon termination of the basic computation by Fortran, an output of the computed results to the external file or several external files in Mathematica format was implemented.And a small program is then written in Mathematica which ensures the graphical visualization of computed results in a form suitable for the user.The exit from the Fortran shell and entry into the Mathematica program needs about 3 seconds of the user's time, and then the Mathematica program is started by pressing the keys Shift/Enter.It reads the arrays of numerical results from the external files created by the Fortran code, and the user can see in a few seconds the computed results in a graphical form on the monitor screen.
Figure 2 shows the results of numerical computation, by Fortran, of the problem of the oblique shock wave reflection from a solid wall, which were obtained by the steadying computation using the method of [5] on a uniform grid of 160x40 cells; the parameter of the method φ = 1/3, which ensures the third-order accuracy of the method in spatial variables in smooth flow subregions.All these graphical results were obtained with the aid of the above-described small Mathematica 8 program, which reads the arrays of numerical results from external files created by Fortran.

2 1 ;
Mathematica program: f[x_]:=x*x/;x>1; expr1=y^2-f[x]*z/.x->Here the first command defines the function in the case when x>1.The system yields a correct result of the symbolic computation in the second command: expr1 =y 2 -4z.Let us now assume that it is necessary to replace the body of the function f[x_] at some stage of the analytic computation with a different body, for example, with the body f[x_]:=-x/;x ≠ It is obvious that this new definition of the function should cancel the foregoing definition of this function.Let us now assume that it is necessary to calculate the expression: expr2=expr1^2-f[x]*z/.x->2; in the next line of the Mathematica code.Mathematica outputs the following answer: expr2 = (y 2 -4z) 2 -4z.This answer is obviously incorrect because f[2]=-2 in accordance with the new definition of function f[x_].That is, Mathematica has ignored the new definition of the function f[x_]:=-x/;x ≠ 1 and continues using the old definition f[x_]:=x*x/;x>1.Below we present the solution of the given problem within the Mathematica framework: For the numerical solution of the given problem, two codes were developed: one in Mathematica and the other in Fortran 90.Both results have coincided.All the computations were carried out on a PC with an 3GHz Intel processor.The Mathematica code was run both in Mathematica 4.1.0and in Mathematica 8.0.4.0.To execute 2600 time steps the Mathematica 4.1.0and Mathematica 8.0.4.0 needed 22457 and 27312 seconds, respectively.Fortran needed only 21.45 seconds for the same computation.Thus, Mathematica 4.1.0computes this two-dimensional task 1047 times slower than the Fortran code.In the case of Mathematica 8.0.4.0, the corresponding slow-down factor is equal to 27312/21.45= 1273.29.The computations of the same problem were also done in Mathematica 4.1.0and in Fortran 90 on a coarser grid of 80x20 cells, and 1600 time steps were done.In this case, Mathematica 4.1.0needed 3455 sec. of CPU time, and the Fortran code needed only 3.36 sec. of CPU time.The corresponding slow-down ratio is 3455/3.36= 1028.28.The discrepancy in the slow-down factors of 1047 and 1028 is explained by the well-known fact that a larger number of the intermediate printouts leads to larger CPU time expenses in Mathematica.

Fig. 1 .Fig. 2 .
Fig. 1.Spatial region in the problem of the oblique shock wave reflection.

Figure 2
Figure 2(b) is the exact solution, and the dashed line is the result of the computation by the method of [5]. Figure 2(c) shows the residual (11) in the case of the method of [5].It can be seen that 2400 time steps are needed to ensure a drop of the solution residual norm to the level of machine round-off errors 10 -10 -10 -12 .