Tools For Mathematical Calculations
Graduate Student Colloquium Talk
February 20, 2008
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.
If you are looking at this on the web: the images are very large because
this is intended to be projected in front of an audience. I can crank
up the size of text, but the pictures stay fixed.
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. Details can be found at http://www.it.iastate.edu/sldb/view.php?id=18.
The instructions there are still for Maple version 10. Version 11 came
out in mid-2007. It is not clear whether the site-licensed software is
version 10 or 11.
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/.
Side remarks: I just discovered last week that some lab computer are
running a versioin called R2006. That brings the number of different
version labeling schemes up to three. My own verson 7.0.4 is also known
as release 12, and the R2006 is probably also called R13 or R14. I don't
know why it is so hard to stick with a single scheme.
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 4.1.2. 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]=
In[3]:= Integrate[(x+1)/(x^2+x+1),{x,0,1}]
Out[3]=
In[7]:= N[%,30]
Out[7]=

|
Maple
> int((x+1)/(x^2+x+1),x);
> int((x+1)/(x^2+x+1),x=0..1);
> evalf[30](%);
|
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)
ans =
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, anyway. 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:
- The Taylor series for sin(ex-1) around the point x=0
up to the x5 term
- 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]=
In[13]:= term1 = Series[f[x+h],{h,0,4}]
Out[13]=
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](tools_figures/mathematica_06.png)
In[15]:= (term1 - term2) / (2h)
Out[15]=
|
Maple
> series(sin(exp(x)-1),x=0,5);

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

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

> simplify(s);
> expand(s);
> combine(s);

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

|
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.
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]=
In[9]:= DSolve[{y'[x] == y[x]+1, y[0] ==
1}, y, x]
Out[9]=
|
Maple
> dsolve(D(y)(x) = y(x)+1, y(x));
> dsolve({D(y)(x) = y(x)+1,y(0)=1}, y(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]](tools_figures/mathematica_figure.gif)
Out[17]= - Graphics - |
Maple
> plot(1/x+exp(x),x=0.1..3);
![[Maple Plot]](tools_figures/maple_figure.gif) |
Matlab
>> x = linspace(0.1,3);
>> y = 1./x + exp(x);
>> plot(x,y)

|
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)

|
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).
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]=
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>⁢</mo> <mrow> <mi>sec</mi> <mo>⁡</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);
> with(CodeGeneration):
> Fortran(s, resultname="s");
> C(s, resultname="s");

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

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

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

> latex(s);
 |
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. |