Matlab and Its Relatives
Graduate Student Colloquium Talk
February 8, 2006
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 three programs that are available at ISU: Matlab, Mathematica, and Maple.

The Programs

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. It is not well suited for extensive numerical computations.

Mathematica is available on four of the machines in the lab in Carver 449 (δ1 through δ4; from the front of the room, they are the leftmost machines in the front row). I also found an old copy on machine ca400w1 (version 2, as opposed to version 5 in the lab).

Maple

Maple is also mainly a symbolic calculator. It is very similar to Mathematica, but poorly supported at ISU. It is available on Project Vincent in directory /afs/iastate.edu/public/maplevr6/bin as either maple (text interface) or xmaple (X windows interface). From machine isua1, the command "add maple" had the desired effect. From machine ca400w1 the "add" did not work, but I was able to start maple from its afs directory. There are also directories maplevr4 and maplevr5 (presumably versions 4 and 5, and the current version is 6). The help files are missing, so you are on your own.

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 other two. The numerical algorithms used are described in detail in the documentation, and for many operations there is a choice of different algorithms. Matlab is available on all machines in the department.

Matlab is also available on Project Vincent, but they let their site license lapse several years ago. You can get to version 6, installed in 2000, by attach matlab6. The current version available in the Math Department is version 7.

 

  Mathematica Maple Matlab
Main Application symbolic symbolic numeric
Support and Availability few copies around, but up-to-date; public copies are for Windows XP only available for several types of workstations (including Linux PC); out of date, missing help files available everywhere; widely used by faculty and in classes; up to date
Output formatting mathematical typesetting text-based or mathematical typesetting (depending on whether you invoke maple or xmaple) text-based only
Save work in internal format yes yes yes
Save/print a transcript of your work text, TeX, HTML, MathML text, TeX, HTML text only
Save formulas as TeX, Fortran, C TeX (supposedly also Fortran and C, but I couldn't get it to work) TeX, Fortran, C
Interface to other programs ?? ?? Fortran, C, Maple, Excel, lab instruments, ...
Plot yes yes yes
Calculate values of special functions (Bessel functions, error function, hypergeometric functions, etc.) yes yes yes

Problem 1: Integration

We want to integrate (x+1)/(x^2+x+1), either as an indefinite integral or as a definite integral from 0 to 1. For the definite integral, also state 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.

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

 

Problem 2: Taylor Series Expansion

We want to find the Taylor series for sin(e^x-1) around the point x=0 up to the x^5 term, and also the expansion of the central finite difference approximation to the first derivative.

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

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.

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

 

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

(I accidentally used the initial condition y(0)=0 here)

In[8]:=

DSolve[y '[x]  y[x] + 1, y, x]

Out[8]=

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

In[9]:=

DSolve[{y '[x]  y[x] + 1, y[0]  0}, 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)

 

Problem 4: Plotting

We want to plot the function 1/x + e^x 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)

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

Mathematica and Maple saved the figures as gif files by default. I am sure you can save them in a variety of other formats as well. For inclusion into TeX documents, you should save them as .eps files (encapsulated Postscript). Matlab needs to be told the format explicitly (this one is a .jpg file), but it can handle a lot of different formats, including .eps and .epsc (color Postscript).

Outputting Fortran or TeX Code

You can convert a single 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)

 

Maple

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

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

> fortran(s);

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

> latex(s); 

1+{\frac {{x}^{2}+1}{\cos(x)}}

> C(s);

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

The Fortran output form is listed in the manual, but it did not work on the version of Maple I tried. The syntax for the C form was just a guess on my part, which also did not work.

 

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

Last Updated: March 7, 2006