Excercise 2

Plotting

Simple plots

Plot the following functions
1. sin(x)
2. cos(x)
3. (e^x)
4. Random[]
5. sin(cos(x)) and (e^(sin(x)))in the same figure.
6. (Subscript[J, 0])(x),  (Subscript[J, 1])(x),  (Subscript[J, 2])(x), and (Subscript[J, 3])(x) where (Subscript[J, n])(x)is the Bessel function of order n, given in Mathematica by the function BesselJ[n,x].
7. Try plotting the factorial function, n!.

Solution

Plot[
      Sin[x], {x, -Pi, Pi}]

Plot[
      Cos[x], {x, -Pi, Pi}]

Plot[
      Exp[x], {x, -1, 1}]

The function Random[] requires no argument (allthough you can give it some). For each point to be plotted, Mathematica computes a random number. When increasing the value of option PlotPoints, the graph of Random[] approaches a black box.

Plot[
      Random[], {x, 0, 1}]

Plot[{Sin[Cos[x]], 
        Exp[Sin[x]]}, {x, -Pi, Pi}]

In the following, I will plot four Bessel functions in the same figure and distinguish them by color, using the option PlotStyle. The graphics directives, such as RGBColor, will be applied cyclically to the list of functions to be plotted.

Plot[{BesselJ[0, x], 
        BesselJ[1, x], BesselJ[2, x], BesselJ[3, x]}, {x, 0, 5}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
          RGBColor[0, 0, 1], RGBColor[0, 0, 0]}]

Plot[n!, {n, -5, 
        5}]

For non-integer arguments of the factorial function, Mathematica automatically uses the generalized factorial function, i.e. the gamma function instead. See the help file for Factorial and Gamma for more information.

Parametric plots

1. Plot the curve x(t)  =  (t^2) + 10  t + 1,  
      y(t) = (1/((t^2) + 10))in x,  y -coordinates.

2. Plot the curve (x(t),  y(t))  =  (cos(2 t), 
          sin(π t))in x, y -coordinates.

3*. Plot the function z = (e^(i t)),  
      0  <  t  <  2  π on the complex z-plane, that is, in (Re z,  Im z) coordinates. The functions Re[...] and Im[...] give the real and imaginary parts of a variable

Solution

ParametricPlot[{t*t + 10  t + 1, 
        1/(t*t + 10)}, {t, -15, 7}]

ParametricPlot[{Cos[2  t], 
        Sin[Pi t]}, {t, -Pi, Pi}]

This is an example of a chaotic system. Try increasing the limits of t. You will see that the curve slowly fills the entire unit square.

ParametricPlot[{Re[Exp[I t]], 
        Im[Exp[I t]]}, {t, -Pi, Pi}]

A complex number can be represented as z=x+yi,and graphically in the complex plane,where the real part of z is its x-coordinate and the imaginary part is the y-coordinate.Furthermore, the exponent function for a complex argument is (e^z) = 
        cos x + i sin x, so the above parametric graph represents the unit circle in xy-plane.

3D and contour plots

1. Use Plot3D[...] and ContourPlot[...] to plot the following functions:
	a. (1/((x^2) + (y^2))), -5 < x < 5,  -5 < y < 5
	
	b. sin(x)cos(y)

Solution

Plot3D[
      1/(x*x + y*y), {x, -2, 2}, {y, -2, 2}]

    ContourPlot[
      1/(x*x + y*y), {x, -2, 2}, {y, -2, 2}]

Plot3D[
      Sin[x] Cos[y], {x, -Pi, Pi}, {y, -Pi, Pi}]

    ContourPlot[
      Sin[x] Cos[y], {x, -Pi, Pi}, {y, -Pi, Pi}]

Plot options

1. Using Mathematicas help system, find out how you can make coloured plots. Try it by plotting sin(x) and cos(x) to same figure with different colours. Can you also find a way to change the thickness of the lines?

2. Plot cos(2000 sin(x)), where 0 < x < π. Find out how you can increase the number of points used in the plot to produce a better plot.

Solution

Options[Plot]

?? PlotStyle

See the help file for PlotStyle. Options for Plot are used like this: (try also Hue instead of RGBColor)

Plot[{Sin[x], 
        Cos[x]}, {x, -Pi, Pi}, 
      PlotStyle -> {{RGBColor[1, 0, 0], 
            Thickness[0.001]}, {RGBColor[0, 1, 0], 
            Thickness[0.01]}}]

Plot[
      Cos[2000  Sin[x]], {x, 0, Pi}]

Pull the graph larger to see the effect of a small number of plotting points.

Plot[Cos[2000  Sin[x]], {x, 0, Pi}, PlotPoints -> 300]

Additional excersises

1. Mathematica knows many special functions. Try plotting some of the following:
	LegendreP[n,x], try with some integer values of n.
	HermiteH[n,x]
	LaguerreL[n,x]
	MathieuC[a,q,x], try for example  a = 3,  q = 2
	Zeta[s,a]
If you're interested, you can use the Mathematica help system to get more information about these functions.

2. Parametric 3D plot. Surfaces in three dimensions can be plotted by parametrizing the surface with two variables. As an example, try plotting the Möbius strip, given by the equations
	x(t, s)  =  [R + s cos((1/2) t)] 
          cos t
y(t, s) =  [R + s cos((1/2) t)] 
          sin t
z(t, s)  =  s sin((1/2) t),	0 ≤ t < 2  π,  -w ≤ 
        s < w,
where R is the radius of the strip and w its width. Use Mathematicas built-in function ParametricPlot3D[{x, y, z}, {s, (Subscript[s, start]),  (Subscript[s, end])}, {t, (Subscript[t, start]),(Subscript[t, end])}], where x,  y,  and z are the coordinates and s and t are the parameters. Try different values of R and w.

Solution

Try different values of n. The polynomial LegendreP[n,x] is an n:th degree polynomial with all its zeros on the interval [-1,1]. You will encounter these functions in quantum mechanics, where they constitute the solution for the angular part of the Hydrogen atom Schrödinger equation.

Plot[LegendreP[9, x], {x, -1, 1}]

Hermite polynomials, which form the solution for the wave function of the harmonic oscillator in quantum mechanics, are also n:th degree.

Plot[HermiteH[4, x], {x, -3, 3}]

Laguerre polynomials solve the radial part of the Hydrogen atom Schrödinger equation.

Plot[LaguerreL[6, x], {x, -1, 10}]

Mathieu functions are needed in the theory of oscillators. Maybe you'll avoid these.

Plot[MathieuC[3, 2, x], {x, 0, 20}]

The Riemann zeta-function is an important function in mathematics, but rarely needed in physics.

Plot[Zeta[x, 4], {x, -5, 5}]

R = 5;  
    w = 2;

    ParametricPlot3D[{(R + s Cos[t/2]) Cos[t], (R + s Cos[t/2]) 
          Sin[t], s Sin[t/2]}, {t, 0, 2  Pi}, {s, -w, w}]

Defining functions

1. Define the function f(x) = cos(10 sin(x))and plot it over x ∈  [0, 2  π]
2. Define the function a(x, y) that calculates the average of two numbers x and y.
3. Define functions f(x) = (x^2), g(x) = exp(-x) and h(x)  =  g(f(x)). Plot the function over -5 < x < 5.
4. Define a function that computes (e^(-((x^2) + (y^2)))) and plot it as a three dimensional plot.
5. Test the difference between delayed and direct assignment.
	i. Set a = 1
	ii. Define the function f1(x) = exp(a x) using direct assignment ('=' operator)
	iii. Define the function f2(x) = exp(a x) using delayed assignment (':=' operator)
	iv. Plot f1(x) and f2(x) for 0 < x < 2.
	v. Set a = 1/2
	vi. Plot f1(x) and f2(x) for 0 < x < 2.

Solution

f[x_] = 
      Cos[10  Sin[x]]

    Plot[f[x], {x, 0, 2  Pi}]

a[x_, y_] = (x + y)/
        2

f[x_] = 
        x*x;

    g[x_] = Exp[-x];

    h[x_] = g[f[x]]

    Plot[h[x], {x, -2, 2}]

f[x_, y_] = 
      Exp[-(x*x + y*y)]

    Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, 
      PlotRange -> {0, 1}]

The option PlotRange determines the range for the z-axis, or y-axis in two-dimensional plots. In the next excercise, the option Hue is used to set the color of the plot. See its help file and try different arguments.

a = 1;

    f1[x_] = Exp[a x];

    f2[x_] := Exp[a x]

    Plot[{f1[x], f2[x]}, {x, 0, 2}, 
      PlotStyle -> {Hue[0.6], Hue[0.1]}]

    a = 0.5;

    Plot[{f1[x], f2[x]}, {x, 0, 2}, 
      PlotStyle -> {Hue[0.6], Hue[0.1]}]

Equations

Polynomial equations

1. Solve the following sets of equations symbolically, and find the numerical values of the solutions:
a.	(x^3) + x + 1 = 0

b. 	(x^2) + (y^2) = 1, 
y = 
        2  x + 1

2. Use Eliminate[...] to eliminate the variable z from the following set of equations
	x + y + z = 0
x y z =  1
(x/z) = 2

3. Use the function Reduce[...] to find the roots of the following equations:
a.	(x^3) + p (x^2) + q (x^2) + r = 0

b.	x y  =  -2
x + y = 1

c.	x y  =  -2
x + y = 1
x  >  0

d.	(x^3) - a (x^2) + x - a = 0
x > 0

4. Find the general solution of a (x^4) +  b (x^3) + c (x^2) + d x + e = 0. The result is really complicated!

5. Find the general solution of (x^2) + p x + 1 = 0 and plot the first solution as a function of p.

6*. Find the general solution of (x^2) + p x + q = 0, make a function f(p, q) that gives the first solution, and finally plot the function as a three dimensional plot, for -15 < p < 15,  -15 < q < 15.

7*. Try what SolveAlways[...] does. Use it like Solve[...] on the second order polynomial equation a (x^2) + b x + c = 0.

Solution

Solve[
      x*x*x + x + 1 == 0, x]

    % // N

Solve[{x*x + y*y == 1, 
        y == 2  x + 1}, {x, y}]

    % // N

Eliminate[{x + y + z == 0, 
        x y z == 1, x/z == 2}, z]

See the help file for Eliminate for more info.

Reduce[
      x^3 + p x^2 + q x + r == 0, x]

    ToRadicals[%]

Here Reduce produces Root-objects, which are roots of a pure function, with argument #1. The Root-objects can be converted into the usual form with the function ToRadicals. This works always for polynomials of order 4 or less.

Reduce[{x y == -2, 
        x + y == 1}, {x, y}]

Reduce[{x y == -2, 
        x + y == 1, x > 0}, {x, y}]

The same as 3a with a further restriction on x.

Reduce[{x^3 - a x^2 + x - a == 
          0, x > 0}, x]

    Reduce[{x^3 - a x^2 + x - a == 0, a > 0}, x]

Note that the symbols are treated as complex numbers unless they appear explicitly in an inequality, or the assumption Element[x,Reals] is made.

Clear[a, b, c, d, e, 
      x]

    Solve[a x*x*x*x + b x*x*x + c x*x + d x + e == 0, 
      x]

Solve[x*x + p x + 1 == 0, 
      x]

    Plot[x /. %[[1]], {p, 2, 6}]

Remember how to apply rules, and extract an element from a list? The function Solve gives a list of rules for x. Then, when we apply the rule to x, we get an expression featuring p and q. Note that the discriminant of the equation is negative for p < 2. Plot will complain if you try to plot it. There's a similar situation in the next excercise.

Solve[x*x + p x + q == 0, 
      x]

    f[p_, q_] = x /. %[[1]]

    Plot3D[f[p, q], {p, -15, 15}, {q, -15, 15}, 
      AxesLabel -> {p, q, f}]

SolveAlways[
      a x^2 + b x + c == 0, x]

SolveAlways tries to find the values of all free parameters appearing in an equation, here a, b and c, for which the equation holds for all values of the variable, here x.

Algebraic equations

1. Use Solve[...] to find the solutions to the following equations
a.	(e^x) = 2

b	(Sqrt[(x^2) - 1]) = (1/2)

c.	(x^y) = 1
x + y = 2

d.	x (e^x) = 1

e.	x exp(exp(x))  =  1

2. Use Reduce[...] to find the solutions to the following equations
a.	((sin(x))^2) + cos(x) = 0

b.	(e^(2  x)) + 2  (e^x) - 
          3 = 0

3. Find the roots of the following equations. To get an initial guess for the solution, plot the right- and left-hand sides of the equations (or right-hand side minus the left-hand side alone) and make an estimate for the position of the root.
a.	cos(x)  =  x

b.	ln x = 1/(x^2)

c.	(Subscript[J, 0])(x) = 0,  5 < x < 6 
	(Subscript[J, 0])(x)is the Bessel J function, BesselJ[0,x] in Mathematica.
	
d. 	Define the function g(x) =  exp(x sin(x)) - x, plot it, and find its root.

Solution

Solve can not in general solve non-polynomial equations. It might give a false answer, or at least leave some solutions out. It will give a warning message when that might be the case.

Solve[Exp[x] == 2, 
      x]

    Solve[Sqrt[x*x - 1] == 1/2, x]

    Solve[{x^y == 1, x + y == 2}, {x, y}]

    Solve[x Exp[x] == 1, x]

    Solve[x Exp[Exp[x]] == 1, x]

Note that the last equation couldn't be solved at all. See the help file for ProductLog if you wish.

Reduce[
      Sin[x]^2 + Cos[x] == 0, x]

    Reduce[Exp[2  x] + 2  Exp[x] - 3 == 0, x]

Reduce never makes approximations. The output corresponds with logical exactness to the input.

Plot[
      Cos[x] - x, {x, -8, 8}]

    FindRoot[Cos[x] == x, {x, 0.5}]

Plot[
      Log[x] - 1/x^2, {x, 0.5, 3}]

    FindRoot[Log[x] == 1/x^2, {x, 1.5}]

In 3.c), since the interval is allready given, we wouldn't have to plot the function to find the initial guess, but let's plot it anyway.

Plot[
      BesselJ[0, x], {x, -10, 10}]

    FindRoot[BesselJ[0, x] == 0, {x, 5, 6}]

The following function has infinitely many zeros, which you can see by increasing the range of x. Let's find the smallest root.

Clear[x]

    g[x_] = Exp[x Sin[x]] - x;

    Plot[g[x], {x, -8, 8}, 
      PlotRange -> {-8, 12}]

    FindRoot[g[x] == 0, {x, 2.5}]

Linear equations

Solve the following linear equations:
a.	(Subscript[x, 1]) + (Subscript[x, 2]) = 1
(Subscript[x, 1]) - (Subscript[x, 2]) = -1

b.	(Subscript[x, 1]) + 2  (Subscript[x, 2]) + 10  (Subscript[x, 3]) = 
        10, 
10  (Subscript[x, 1]) + 2  (Subscript[x, 2]) + 10  (Subscript[x, 3]) = 
        1, 
(Subscript[x, 1]) + 2  (Subscript[x, 2]) + 3  (Subscript[x, 3]) = 0

c.	(({{1, 2, 3}, {0, 1, 2}, {0, 0, 1}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}}))=(({{1}, {1}, {1}}))

d.	(({{1, 1, 1}, {1, 0, -1}, {2, 1, 0}}))(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {(Subscript[y, 3])}}))=(({{0}, {1}, {2}}))

Solution

LinearSolve[{{1, 
          1}, {1, -1}}, {1, -1}]

LinearSolve[{{1, 2, 10}, {10, 2, 
          10}, {1, 2, 3}}, {10, 1, 0}]

LinearSolve[{{1, 2, 3}, {0, 1, 
          2}, {0, 0, 1}}, {1, 1, 1}]

M = {{1, 1, 1}, {1, 
            0, -1}, {2, 1, 0}};

    LinearSolve[M, {0, 1, 2}]

    Det[M]

The last set of equations has no solutions because the determinant of the matrix of coefficients is zero.

Misc. equations

Solve the following equations. Choose an appropriate function (Solve, NSolve, FindRoot, LinearSolve) 
a.	a b  +  c = 0
c + a = 1
a b c = 1, 	a,  b,  c  =  ?

b.	s + 2  t + 3  u = 10
t + u = 1
u - t = 0,	s,  t,  u  =  ?

c.	x log(x) = 10, 	x = ?

d.	sin(2  x) + 3 cos(x) + 1 = 0,	x = ?

Solution

In (a) the equations are all polynomial, so lets use Solve. For comparison, use also Reduce.

Clear[a, b, c]

    Solve[{a b + c == 0, c + a == 1, a b c == 1}, {a, b, 
        c}]

    Reduce[{a b + c == 0, c + a == 1, a b c == 1}, {a, 
        b, c}]

The equations in (b) are all linear.

LinearSolve[{{1, 2, 3}, {0, 1, 
          1}, {0, -1, 1}}, {10, 1, 0}]

In (c), it seems that we'll be needing FindRoot. However, let's try Solve. It gives a solution in terms of ProductLog. And, I'll be damned, the Solve solution is correct!

Solve[x Log[x] == 10, 
      x]

    % // N

    Plot[x Log[x] - 10, {x,  .1, 8}]

    FindRoot[x Log[x] == 10, {x, 6}]