Tools For Mathematical Calculations
Graduate Student Colloquium Talk
February 10, 2009
Fritz Keinert

This talk gives an overview of some programs that can take the tedium out of mathematical calculations. They perform calculations such as differentiation, integration, solving nonlinear equations, ODEs and PDEs, Taylor series expansions, plotting functions, and so on. We will concentrate on four programs that are available at ISU: Mathematica, Maple, Matlab, SciLab.

This is the regular version, intended to be read on the screen. If you want to project this in front of an audience, use the large size version instead.


The Programs, And Where To Find Them

Mathematica

Mathematica is designed to provide symbolic answers, that is, to solve problems in closed form. It can also solve these problems numerically, to arbitrary precision, but gives no detail on how that is done internally. It is not well suited for extensive numerical computations. We have 10 copies running on public machines in the Math Department, all of them recent versions with good support. This is one good option if you need a symbolic calculator.

Mathematica is expensive. A couple of years ago it was $800 per license. The department pays for five public copies in Carver 449. I discovered another five copies in Carver 400, so there are 10 public copies around, with version numbers between 5.0 and 6.0. The latest version is 6.0. I used version 5.2 installed in the lab to generate the examples below.

Windows XP machines: In Carver 449, machines δ1 through δ5 (front row left, if you are standing by the whiteboard), version 5.2. In Carver 400, machines C00179 (version 5.1) and C00180 (version 5.2), and the scanner machine (version 6.0). The C00179 etc. refers to the inventory number listed on a big sticker attached to the computer. You can access Mathematica by looking for it under Start -> Programs. The machines in 449 also have a link on the desktop.

Linux machines: In Carver 400, pollux (version 5.2) and proton (version 5.0). The executables are in /usr/local/Wolfram/Mathematica/5.x/Executables/math/.

In addition, there is a copy of Mathematica 2.2.3 in /afs/iastate.edu/public/mathematica/. I am not sure which kind of machine it runs on, but it is ancient anyway.

The publisher's web site is http://www.wolfram.com/.

Maple

Maple is also a symbolic calculator. It has similar capabilities to Mathematica, and a similar price tag. In early 2007, their web site quoted me a price of $997 for a single user copy, with academic discount. The current version is 11. The examples below are from version 10.

ISU has a site license for students, for version 10 only. Details can be found at http://www.it.iastate.edu/sldb/view.php?id=18. That web page says that support will disappear in June 2010.

An ancient Version 6 of Maple is available at /afs/iastate.edu/public/maplevr6/ for Dec Alphas and for Linux PCs. Public machines of this type in Carver 400 are mercury and duluth (Alphas), and proton and pollux (Linux). I would recommend the Linux machines, they seem to be much faster than the antique Alphas. The executables are called maple (text interface) or xmaple (X windows interface). The help files are missing, so you are on your own.

The publisher's web site is http://www.mapletech.com/.

Matlab

Matlab is designed primarily to provide numerical answers. It can perform some symbolic computations by internally calling Maple, but its symbolic capabilities are limited compared to the Mathematica and Maple. The numerical algorithms used are described in detail in the documentation, and for many tasks there is a choice of different algorithms. Matlab is available on all machines in the Math department. It is used extensively in numerical analysis and other applied math courses. We have a site license that lets faculty install Matlab on their machines, even at home. There is a separate site license for students. Details are given at http://www.it.iastate.edu/sldb/view.php?id=21. The examples below were done in version 7.0.4, which is also the version installed in the labs.

The publisher's web site is http://www.mathworks.com/.

Scilab

Scilab is the open source community's answer to Matlab. Originally developed by INRIA (Institut National du Recherche en Informatique et Automatique) and ENPC (École National des Ponts et Chaussées), it is now overseen by a 23-member consortium. The program, complete with source code, is freely available. The latest version is 5.0.3. It is intended to have the same functionality as Matlab. ISU is encouraging departments to switch to free software (Microsoft Office to Open Office, for example), and they are proposing Scilab as a replacement for Matlab.

The basic syntax of Scilab is very similar to Matlab. For example, if A and B are matrices, then X = A/B means X = AB-1 (or more generally, X is the least squares solution to XB = A), while X = A./B is elementwise division. Matrices are entered as A = [1,2;3,4]. The notation a:b:c means "the numbers from a to c in steps of b", and so on. Some simple Matlab scripts will run unchanged in Scilab. However, names and syntax of many built-in routines are different, subroutine definition is different, and some of Matlab's features may not be provided. Few Matlab routines will run without changes in Scilab. Scilab does provide a tool to convert Matlab routines to Scilab. I haven't tried that yet.

My opinion: Scilab is better than Matlab in some respects. In particular, they made it a lot easier to interface to subroutines written in Fortran or C. This can be done in Matlab, but it is a programming nightmare. (I speak from experience here). The whole program is open source, so you can implement lots of things if you choose to invest the effort. I like the fact that you can set the number of displayed decimals to anything you choose; in Matlab, you either get 5 decimals or 15, nothing else. Scilab is able to read Matlab data files and other formats (Excel spreadsheets, for example).

However, Scilab also has shortcomings. For example, symbolic calculation is only offered for rational functions at this point. I also ran into inconsistent syntax and lack of explanation of the underlying algorithms (see integration example below).

If you need to do a lot of interfacing with Fortran, C or anything nonstandard for your research, Scilab is worth a look. For teaching purposes I would stick with Matlab at this point.

Scilab is not widely installed in the department, but if you want it, all you have to do is ask. Scilab can be downloaded from http://www.scilab.org/.

Other Programs

Macsyma is symbolic calculator with capabilities similar to Mathematica and Maple. We had it running in the department 15 years ago. At that time it was very clumsy to use, and obviously written by many people who didn't talk to each other: different portions of it had totally different syntax. I think they cleaned it up considerably since then, and it is probably a nice product now, but I am not aware of any copies of it on campus.

There are other symbolic calculators I have heard of: Singular and Derive. I don't know anything about their capabilities. They may be more specialized. We used to have Gap, which is a group theory calculator. It may still be around.

Octave is the GNU Project's version of Matlab: open source, freely available, intended to be fully compatible with Matlab. I looked at it once, but not in great detail. My impression was that it is equivalent to an earlier version of Matlab, before they added the fancy graphics, data types, etc.


 

  Mathematica Maple Matlab Scilab
Main Use symbolic computation symbolic computation numeric computation numeric computation
Support and Availability Availability is moderate (10 public copies).
Support is good (the copies we have are recent and come with complete manuals)
Availability for students is good. Availability for faculty is nonexistent.
Support is good (extensive online help system)
Available everywhere; widely used by faculty and in classes.
Good Support (extensive online help, many knowledgeable people around)
Available everywhere.
Support is weak (the online help is quite basic).
Output formatting text-based or typeset text-based or typeset text-based only text-based only
Save work in internal format yes yes yes yes
Save/print a transcript of your work text, TeX, HTML, MathML text, TeX, HTML text only text only
Save formulas as TeX, Fortran, C, MathML TeX, Fortran, C, Java, Matlab TeX, Fortran, C TeX (rational functions only)
Fortran (subroutines only)
Interface to other programs ?? ?? Fortran, C, Maple, Excel, lab instruments, ...
Not easy but possible
Fortran, C, whatever else you feel like programming yourself
Very extendable
Plotting capabilities yes yes yes yes
Calculate values of special functions (Bessel functions, error function, hypergeometric functions, etc.) yes yes yes yes

Problem 1: Integration

Integrate (x+1)/(x2+x+1), both as an indefinite integral and as a definite integral from 0 to 1. For the definite integral, also find the answer as a decimal with 30 digits.

Mathematica

In[2]:= Integrate[(x+1)/(x^2+x+1),x]

Out[2]=

ArcTan[(1 + 2 x)/3^(1/2)]/3^(1/2) + 1/2 Log[1 + x + x^2]

In[3]:= Integrate[(x+1)/(x^2+x+1),{x,0,1}]

Out[3]=

1/18 (3^(1/2) π + 9 Log[3])

In[7]:= N[%,30]

Out[7]=

0.851606038373091154129968994735

Maple

> int((x+1)/(x^2+x+1),x);

1/2*ln(x^2+x+1)+1/3*sqrt(3)*arctan(1/3*(2*x+1)*sqrt...

> int((x+1)/(x^2+x+1),x=0..1);

1/2*ln(3)+1/18*sqrt(3)*Pi

> evalf[30](%);

.851606038373091154129968994735

Matlab

Matlab is primarily a numerical calculator. We do the definite integral numerically first, then the symbolic integral. Note that the default accuracy is fairly low.

>> help quad
QUAD Numerically evaluate integral, adaptive Simpson quadrature.
Q = QUAD(FUN,A,B) tries to approximate the integral of function
FUN from A to B to within an error of 1.e-6 using recursive
adaptive Simpson quadrature. The function Y = FUN(X) should
accept a vector argument X and return a vector result Y, the
integrand evaluated at each element of X.
...
See also quadv, quadl, dblquad, triplequad, @, trapz.

>> f = inline('(x+1)./(x.^2+x+1)')
f =
Inline function:
f(x) = (x+1)./(x.^2+x+1)
>> quad(f,0,1)
ans =
0.85160606414482
>> quad(f,0,1,1.e-15)
ans =
0.85160603837309

>> syms x
>> int((x+1)/(x^2+x+1))
ans =
1/2*log(x^2+x+1)+1/3*3^(1/2)*atan(1/3*(2*x+1)*3^(1/2))
>> I = int((x+1)/(x^2+x+1),0,1)
I =
1/2*log(3)+1/18*3^(1/2)*pi
>> double(I)
ans =
0.85160603837309
>> vpa(I,30)
ans =
.851606038373091154129968994735

Scilab

Scilab has extremely limited symbolic tools. I wasn't able to find a symbolic integration routine. I only did the numerical integration.

-->function y = f(x), y = (x+1)/(x^2+x+1), endfunction
-->format('v',30)
-->[I,err] = intg(0,1,f)

err =
0.0000000000000094547263
I =
0.8516060383730911231837
-->integrate('(x+1)/(x^2+x+1)','x',0,1)
ans =
0.8516060383730911231837
-->ans - I
ans =
0.

Observations:

  • The program is fairly easy to use, and gives good accuracy. The documentation says that intg uses dqag0.f and dqags.f from quadpack. No algorithm is given for integrate, but it seems to be the same one.
  • I tried to set the accuracy to 30 decimals, but it refuses to display more than 22, and the last 6 are garbage. Obviously SciLab works in double precision exclusively. No extended accuracy.
  • There are two ways to invoke the integration routine. In one case, the arguments are (f,a,b), in the other case (a,b,f). They should be the same.

Problem 2: Taylor Series Expansion

We want to find two things:

  1. The Taylor series for sin(ex-1) around the point x=0 up to the x5 term
  2. The expansion of the central difference [f(x+h) - f(x-h)] /(2h) of an unspecified function f(x) in powers of h.

Mathematica

In[11]:= Series[Sin[Exp[x]-1],{x,0,5}]

Out[11]=

x + x^2/2 - (5 x^4)/24 - (23 x^5)/120 + O[x]^6

In[13]:= term1 = Series[f[x+h],{h,0,4}]

Out[13]=

f[x] + f^′[x] h + 1/2 f^′′[x] h^2 + 1/6 f^(3)[x] h^3 + 1/24 f^(4)[x] h^4 + O[h]^5

In[14]:= term2 = Series[f[x-h],{h,0,4}]

Out[14]=

f[x] - f^′[x] h + 1/2 f^′′[x] h^2 - 1/6 f^(3)[x] h^3 + 1/24 f^(4)[x] h^4 + O[h]^5

In[15]:= (term1 - term2) / (2h)

Out[15]=

f^′[x] + 1/6 f^(3)[x] h^2 + O[h]^4

 

Maple

> series(sin(exp(x)-1),x=0,5);

series(1*x+1/2*x^2-5/24*x^4+O(x^5),x,5)

> term1 := series(f(x+h),h=0,4);

term1 := series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x)*...

> term2 := series(f(x-h),h=0,4);

term2 := series(f(x)+(-D(f)(x))*h+1/2*`@@`(D,2)(f)(...

> s := (term1 - term2)/(2*h);

s := 1/2*((series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x...

> simplify(s);

-1/2*(-(series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x)*h...

> expand(s);

1/2*(series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x)*h^2+...

> combine(s);

1/2*((series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x)*h^2...

> collect(s,h);

(1/2*(series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x)*h^2...

> normal(s);

-1/2*(-(series(f(x)+D(f)(x)*h+1/2*`@@`(D,2)(f)(x)*h...

> series(s, h=0,4);

Maple equation 23

The last example shows why I tried to learn Maple two or three times, and gave up every time: I just could not convince Maple to simplify the result in the form that I wanted. The final simplification (the one that actually worked) is courtesy of Olga Ruff, a former graduate student.

Matlab

>> taylor(sin(exp(x)-1),5)
ans =
x+1/2*x^2-5/24*x^4

>> syms f x h
>> term1 = taylor('f(x+h)',h,4)

term1 =
f(x)+D(f)(x)*h+1/2*@@(D,2)(f)(x)*h^2+1/6*@@(D,3)(f)(x)*h^3
>> term2 = taylor('f(x-h)',h,4)
term2 =
f(x)-D(f)(x)*h+1/2*@@(D,2)(f)(x)*h^2-1/6*@@(D,3)(f)(x)*h^3
>> s = (term1 - term2)/(2*h)
s =
1/2*(2*D(f)(x)*h+1/3*@@(D,3)(f)(x)*h^3)/h

Scilab

I couldn't find anything in Scilab related to power series expansions.


Problem 3: Solving a Simple ODE

We want to solve the ODE y' = y+1, with or without the initial condition y(0)=1.

Mathematica

In[8]:= DSolve[y'[x] == y[x]+1, y, x]

Out[8]=

output

In[9]:= DSolve[{y'[x] == y[x]+1, y[0] == 1}, y, x]

Out[9]=

{{yFunction[{x}, -1 + ^x]}}

 

Maple

> dsolve(D(y)(x) = y(x)+1, y(x));

y(x) = -1+exp(x)*_C1

> dsolve({D(y)(x) = y(x)+1,y(0)=1}, y(x));

y(x) = -1+2*exp(x)

 

Matlab

We do the numerical solution first, then the symbolic solution.

>> yprime = inline('y+1','t','y')
yprime =
Inline function:
yprime(t,y) = y+1
>> [t,y] = ode45(yprime,[0,3],1);
>> length(t)

ans =
45
>> true_y = -1 + 2*exp(t);
>> max((y - true_y)./true_y)

ans =
1.787437079083398e-006
>> options = odeset('RelTol',1.e-10);
>> [t,y] = ode45(yprime,[0,3],1,options);
>> length(t)

ans =
89
>> true_y = -1 + 2*exp(t);
>> max((y - true_y)./true_y)

ans =
1.651705715033635e-007

>> dsolve('Dy = y+1')
ans =
-1+exp(t)*C1
>> dsolve('Dy = y+1','y(0)=1')
ans =
-1+2*exp(t)

Scilab

As before we can only do the numerical solution, not the symbolic solution.

-->function yprime = f(t,y), yprime = y + 1, endfunction
-->y0 = 1;
-->t0 = 0;
-->t = 0:0.1:3;
-->y = ode(y0,t0,t,f);
-->true_y = -1 + 2*exp(t);
-->max((y - true_y)./true_y)

ans =
0.0000002201200831737934

Comments:

  • I just used the default algorithm here (lsoda from ODEPACK), but many algorithms are available. There are also many other options you can set.
  • Matlab by default returns the solution at whatever points it was computed. SciLab insists on interpolating the solution at a sequence of points the user has to provide.

Problem 4: Plotting

We want to plot the function 1/x + ex over the interval [0.1, 3].

Mathematica

In[17]:= Plot[1/x + Exp[x], {x, 0.1, 3}]

[Graphics:matlab_talk_figures/mathematica_37.gif]

Out[17]= - Graphics -

 

Maple

> plot(1/x+exp(x),x=0.1..3);

[Maple Plot]

 

 

Matlab

>> x = linspace(0.1,3);
>> y = 1./x + exp(x);
>> plot(x,y)

ddd

Scilab

The code is virtually identical to Matlab, with one strange exception: y=1./x should be the elementwise inverse of x. Scilab instead treated it as y=1/x, which is the least squares solution of yx=1. I had to write it as y=x.\1. The reason is probably that A./B requires A and B to have the same size, and Scilab is not smart enough to expand 1 into a vector of 1s. Matlab used to have the same problem, but at least it gave an error message at that point, instead of quietly doing something different. Looks like a bug in Scilab to me.

The plot looks remarkably like Matlab's plot, including the syntax for producing it. I think the syntax for making dotted lines etc. is also the same.

-->x = linspace(0.1,3);
-->y = x.\1 + exp(x);
-->plot(x,y)

scilab plot

Comparison

These are the default figures for all programs. I am sure you can add labels, change the axis limits, etc., for all of them.

All of the programs can save the files in a variety of formats (gif, jpg, ps,...). For inclusion into TeX documents, you should save them as .eps files (encapsulated Postscript), and then convert them to PDF.


Outputting Fortran or TeX Code

We want to convert a formula to Fortran code, C code, or TeX code, for cutting and pasting into a program or TeX document.

Mathematica

In[21]:= s = 1 + (x^2+1)/Cos[x]

Out[21]=

1 + (1 + x^2) Sec[x]

In[22]:= FortranForm[s]

Out[22]//FortranForm= 1 + (1 + x**2)*Sec(x)

In[23]:= TeXForm[s]

Out[23]//TeXForm= \left(x^2+1\right) \sec (x)+1

In[24]:= CForm[s]

Out[24]//CForm= 1 + (1 + Power(x,2))*Sec(x)

In[25]:= MathXLForm[s]

Out[25]//MathXLForm=

<math>
<mrow>
<mrow>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mo>)</mo>
</mrow>
<mo>&#8290;</mo>
<mrow>
<mi>sec</mi>
<mo>&#8289;</mo>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>+</mo>
<mn>1</mn>
</mrow>
</math>

 

Maple

> s := 1 + (x^2+1)/cos(x);

s := 1+(x^2+1)/cos(x)

> with(CodeGeneration):
> Fortran(s, resultname="s");

fortran(1+(x^2+1)/cos(x))

> C(s, resultname="s");

C(1+(x^2+1)/cos(x))

> Java(s, resultname="s");

C(1+(x^2+1)/cos(x))

> Matlab(s, resultname="s");

C(1+(x^2+1)/cos(x))

> VisualBasic(s, resultname="s");

C(1+(x^2+1)/cos(x))

> latex(s);

latex

 

Matlab

>> syms x
>> s = 1 + (1+x^2)/cos(x)

s =
1+(1+x^2)/cos(x)
>> pretty(s)

                  2
1 + x
1 + ------
cos(x)

>> fortran(s)
ans =
t0 = 1+(1+x**2)/cos(x)
>> latex(s)
ans =
1+{\frac {1+{x}^{2}}{\cos \left( x \right) }}
>> ccode(s)
ans =
t0 = 1.0+(1.0+x*x)/cos(x);

Scilab

As I mentioned before, Scilab can only do rational functions in symbolic form. Those can be converted to TeX, but not Fortran or C:

-->p = poly([1,1],"x","coef")

p =
  1 + x
-->q = poly([1,0,2],"x","coef")
q =
        2
  1 + 2x
-->r = p/q
r =
  1 + x
  -----
        2
  1 + 2x
-->A = [1,r;2,3]
A =
  1 1 + x,
  - -----
          2
  1 1 + 2x
 
  2 3
  - -
  1 1
-->texprint(A)
ans =
    {\pmatrix{ 1&{ 1+x}\over{ 1+2\*x^{2}}\cr 2&  3}}

There is also a tool for converting a Scilab subroutine to Fortran. I didn't try it.

 


Last Updated: February 10, 2009