Math 21a: Multivariable Calculus, Fall 2000
Mathematica Assignment 3

Introduction and Review

The purpose of this Mathematica assignment is to introduce you to some of the differential equation solving abilitites of Mathematica. However,before starting, you might find it useful to review the introductory material from the first two Mathematica assignments about the special notation and formatting rules for Mathematica. In any event, here is a brief summary of the more ubiquitous ones:

Entering Mathematica Expressions
To create a new input cell, position the cursor between two existing cells so that Mathematica switches the cursor to a horizontal I-beam, click the mouse, and then start typing. If you type into an existing output cell, Mathematica will make a copy of the cell and allow you to edit the content. If you have two expressions in the same input cell, separate them by ; so Mathematica knows that they are two different expressions and should be evaluated separately.

Evaluating Mathematica Expressions
To evalute an expression or to have Mathematica input or read an expression, click anywhere in the input cell and press simultaneously the 'shift' and 'enter' key or the 'enter' key on the keypad.

Mathematica Syntax
Capitalize first letter of built-in  functions and variables.
	Correct:     Sin[x], E^Pi 	                                            Incorrect:  sin[x], e^pi, LOG[e]
Use () to group expressions and [] to enclose function arguments.
	Correct:    (1+2)*3, Sin[Pi], E^Log[5] 	Incorrect:  {[1+2]*3}, Sin(pi), e^ln 5
Use [], x_, and := when defining functions.
	Correct:     f[x_]:=g[x]^2+3x		                             Incorrect:   f(x)=g(x)^2+3x
Remember to evaluate function definitions by pressing  simultaneously the keys 'shift' and 'enter' or just the key 'enter' on the keypad over the cell before using them.

Implicit Multiplication
Separate expressions with a space to denote multiplication.
	Correct: (-b+Sqrt[b^2-4a c])/(2a)		Incorrect: (-b+Sqrt[b^2-4ac])/(2a)

Mathematica Notebook and Outlines
A notebook is a Mathematica document.  It is organized as an outline as indicated by the blue brackets on the right.  You can double-click  the blue brackets to collapse and expand a particular section of the outline.  Thus, if your Mathematica notebook all of a sudden "disappears," most likely, you have collapsed the entire notebook.  Just double click the right-most  blue bracket to expand the outline again.

Differential equations for functions of time

The simplest differential equations are those for a function,  x(t) of one variable t. Of these perhaps the simplest is   x'= f(x) 
where f is a function. For example, x(t)= x,  or  x'= x (1 - x). Equations of the form  x' = a x with a=constant are called "exponential  growth"  equations. To see what a solution looks like with a=2, have Mathematica evaluate the following:

New = 
        NDSolve[{x'[t] == 2 x[t], x[0] == 1}, 
          x, {t, 0, 1}];
Plot[x[t]  /.  New, {t, 0, 1}];

Explanation: The special Mathematica function NDSolve solves differential equations for functions x(t). It has the form A= NDSolve[{equation to be solved, initial value for x}, x, {t, start time, end time}]. (The equal sign == which should be read here as identical to has not to be confused with the equal sign =, which is used to assign values. For example x=5 defines x to be 5. Asking then x==5 gives the answer True.) Applying NDSolve tells Mathematica to solve the equation for x(t) from the start time to the end time with the given starting value for x. The result A is a rule pointing to an interpolating function. We get a function in the variable t from this rule with x[t] /. A[[1]] . This function is plotted with the command Plot[x[t] /. A[[1]], {t, start time, end time}];

Logistic = 
        NDSolve[{x'[t] == 2 x[t]*(1 - x[t]), 
            x[0] == 0.01}, x, {t, 0, 5}];
Plot[x[t]  /.  Logistic[[1]] , {t, 0, 5}];

To gain practice, change the starting value for the Logistics equation to x[0] = -1 and graph the result.  Also, try x[0] = 0.5 and x[0] = 2.  Do you remember enough from Harvard's Math 1b (or its equivalent)  to explain the behavior of x[t] for these various starting values? For a second example, change the equation in the NDSolve entry to 
x'[t]=Sin[x[t]] and plot solutions with x[0]=π/4, 2π/4,5π/4 and 9π/4.

Differential equations of second order

The Mathematica's NDSolve function is not limited to solutions of equations which involve only one derivative of x.  For example, it can handle equations of the form  x'' + a x' + b x = 0, where a and b are constants, and many equations which relate derivatives of even higher order in x to its lower order derivatives; for example x''' + x'' sin(x') + cos(x'' x) = 0.

However, to keep matters simple, we consider equations of the form x'' + a x' + b x = 0 which you might recall seeing in Math 1b or its  equivalent. They are called differential equations of second order. For example, to plot the solution to x'' + x' + x = 0  with initial conditions x(0) = 1 and x'(0) = 1, ask Mathematica to evaluate:

SecondOrder = 
        NDSolve[{x''[t] + x'[t] + x[t] == 0, 
            x[0] == 1, x'[0] == 1}, 
          x, {t, 0, 10}];
Plot[x[t] /. SecondOrder[[1]], {t, 0, 10}];

Comments: 
With regard to the formatting here, note that two initial conditions  must be specified for these equations with two derivatives on x. This is a fairly general rule:  For an equation of the form  n-derivatives  of x = function of derivatives of order less than n, you must specify the initial values of x and its derivatives up to order n-1 or else there won't be a unique solution. On the other hand, if you specify more than these n initial conditions, then there usually isn't any solution. By copying the preceding cell and changing the equation appropriately, have Mathematica plot solutions to x'' +  x' + x = 0, x'' + 2 x' + x = 0  and x'' + 3 x' + x = 0,  x'' - x' + x = 0  all using the same initial conditions.  
You might notice the different behavior; the initial case corresponds to a damped oscillator, the second case to a critically damped, the third to an over damped oscillator, the fourth describes a negative damping  situation. Finally, pick your own values for a and b in x''+a x'+b=0 and your own initial conditions and have Mathematica plot the solutions.

Differential Equations for vector valued functions

Mathematica's NDSolve function will also solve differential equations for vector functions of time.  A typical example which involves the first derivative only asks for vector function of time, r(t), which  solves the equation r'(t) = F(r(t)); where F is some given vector field. Of course if r(t) = (x(t), y(t)) has two components, then the graph of the solution is 3-dimensional since you have to have t, x and y pictured on the same plot.  And, if r(t) = (x(t), y(t), z(t)) has three components, then there is no way to graph the solution as a function of time with all three coordinates pictured.  (You could just graph versus time the behavior of the x and y coordinates, or just the x and z, or just the y and z.)  However, in either case, Mathematica can plot for you the curve traced out by the solution. 
For example, to plot the curve traced out by the solution to the equation r´(t) = (-2 x(t) + y(t), -2 y(t) - x(t)), with starting position (x(0) = 1, y(0) = 1), ask Mathematica to evaluate

VectorEq = 
        NDSolve[{x'[t] == -2 x[t] + y[t], y'[
                t] == -2*y[t] - x[t], x[0] == 1, 
            y[0] == 1}, {x, y}, {t, 0, 
            4}];
ParametricPlot[
        Evaluate[{x[t], y[t]}] /. VectorEq[[1]], {t, 0, 
          4}];

Comment:
With regard to the formatting of the expressions in the preceding cell, note that the NDSolve[] function has, except for the MaxSteps entry, essentially the same form as before,

NDSolve[{ equations, initial conditions},{ variable names},{ t, start time, end time}, options} .

Mathematica solve the equation by breaking the time interval into some number, N, of steps, occuring at some finite subset of times (Subscript[t, j]) and then replacing the differential equation by a difference equation. The  simplest approximation, called Euler steps, would be r((Subscript[t, j])+1) = r((Subscript[t, j])) + F(r((Subscript[t, j]))) ((Subscript[t, j + 1]) -(Subscript[t, j])), which clearly indicates how to obtain r((Subscript[t, j + 1])) from r((Subscript[t, j])). So, you must tell Mathematica the initial value, r(0),  or it won't know how to start the process. (Mathematica uses actually much more sophisticated numerical methods than the Euler step method just mentioned and even choses the method which it considers the best for the equations.)

Suggestion:
Experiment by asking Mathematica to plot the trajectories for solutions of equations of the form 
((Subscript[x, t]) = a x + b y, (Subscript[y, t]) = c x + d y) for four or five of your own choices for the constants a, b, c and d and for the starting point in the x-y plane.  For example, take the original version above and start it at (1,-1),  then (-1,-1) and finally (-1,1).  
By the way, if you choose small values for the constants a,b,c,d, you may have to increase the time interval to see more of the trajectory.

The Lorentz equation:
The Lorenz equation is a famous equation for three dimensional vector  (x(t),y(t),z(t)).  This equation is said to exhibit chaotic properties as the trajectories of solutions which start at t = 0 very close to each other evolve in drastically different ways.  A version of the Lorenz equation is 
(x'(t), y'(t), z'(t)) = (10 (y - x), - x z + 28 x - y, x + y - 8 z/3). 
To view a three dimensional plot of the trajectory of the solution starting at (0,0.4,0), have Mathematica evaluate the following cell.

Lorentz = 
        NDSolve[{x'[t] == 
              10 (y[t] - x[t]), y'[t] == -x[t] z[t] + 
                28  x[t] - y[t], z'[t] == x[t]*y[t] - 8  z[t]/3, 
            x[0] == z[0] == 0, y[0] ==  .4}, {x, y, z}, {t, 
            0, 40}, MaxSteps -> 5000];
ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. Lorentz], {t, 0, 40}, 
        PlotPoints -> 5000];

(With regards to the formatting: The MaxStep -> 5000 entry tells Mathematica to increase the default maximal number of step (which is 1000) to 5000. In the ParametricPlot[ ] entry, the PlotPoints -> 5000 entry tells Mathematica to use 5000 points to plot the trajectory.) Now have Mathematica plot the solutions for various other intial conditinos, say (1,0,0), (0,1,0), (0,0,1), (1,1,1) or four of your own choosing. Notice that the trajectories all seem to converge upon the same, convoluted subset of space. This part of space is called the `Lorentz attractor' as it attracts trajectories.

Just for fun, we ask Mathematica to show the endpoint (x(40),y(40),z(40)) at the time t=40 as a function of the starting position (a,b,c)=(0,y(0),0), where we vary b=y[0]:

EndPoint[a_, b_, c_] := 
        Module[{L, f}, 
          L = NDSolve[{x'[t] == 
                  10 (y[t] - x[t]), y'[t] == -x[t] z[t] + 
                    28  x[t] - y[t], z'[t] == 
                  x[t]*y[t] - 8  z[t]/3, x[0] == a, y[0] == b, 
                z[0] == c}, {x, y, z}, {t, 0, 40}, 
              MaxSteps -> 5000]; 
f[
              t_] := Evaluate[{x[t], y[t], z[t]} /. L][[1]];  
          f[40]];

(With regards to the formatting: Module[{variables},  commands  ] is used to make a procedural unit from a set of  commands. The variables in the first bracket are local variables and not accessible from outside.)

You can observer that even a tiny change of initial conditions of (10^(-16)) (which in  centimeters is 1 million'th of the size of an atom) leads to a completely different end point:

EndPoint[0, 0.4, 0]
EndPoint[0, 0.4000000000000001, 0]

This sensitive dependence on initial conditions is also called "chaos". Because of its fractal shape of the orbit and this sensitivity, Chaologists (scientists studing chaos) call the Lorentz attractor a strange attractor. It is used to illustrate the "butterfly effect": A small change of initial conditions of the atmosphere like the disturbance of a wing of a butterfly could produce a storm in a couple of days. Actually, the Lorentz equations were derived originally from a physical model for the atmosphere.

The advection equation

Mathematica  can also plot solutions to partial differential equations for functions u(t, x) which depend on time t, and a spacial coordinate x. One of the simplest such equation is the advection equation (Subscript[u, t])=-c (Subscript[u, x]), where c is a constant. Solutions to this equation describe travelling waves. Moreover, given any initial function g(x), there is a unique solution to the advection equation, u(t,x), with the property that u(0, x) = g(x).  To see an example solution, have Mathematica  evaluate the following:

Adv  = 
        NDSolve[{D[u[x, t], t] == -2 D[u[x, t], x], 
            u[x, 0] == 1/(1 + x^2), u[-20, t] == 1/401, 
            u[20, t] == 1/401}, 
          u, {x, -20, 20}, {t, 0, 9}];
Plot3D[
        Evaluate[u[x, t] /. Adv [[1]]], {x, -20, 20}, {t, 0, 9}, 
        AxesLabel -> {"x", "t", "u[t,x]"}, 
        PlotRange -> {0, 1}, PlotPoints -> 100, Mesh -> False, 
        Boxed -> False];

Comments:
Here are some comments about the formatting of the expressions in the preceding cell: First off, notice that the NDSolve function is still being used, even for partial differential equations. In this case, its formatting has the following schematic form:
NDSolve[{equation for function, function at the start time, constraints on boundary values}, name of the function, {x, lower bound for x, upper bound for x}, {t, start time, end time}]. You should also notice how time and space derivatives are indicated in Mathematica.  In particular, the operation D[ ] has entries D[function,derivative direction],  so D[u[t,x], t] tells Mathematica to take the partial derivative of u[t,x] in the t-direction, while D[u[t,x],x] tells Mathematica to take the partial derivative in the x-direction. You should also notice that the constant c in this example of the advection equation is set equal to 2. Notice as well that the input for the NDSolve[ ] command requires not just the initial value, u[0,x], but the value of u at an end point of its range.  Without this input, there wouldn't be a unique solution and Mathematica would chose its own boundary condition. By the way, you can set this  boundary condition for u[t,-20]  to be any function of t which is consistent with the choice of u[0,-20]. We choose the constant function u[t,-20] = u[0,-20]  for simplicity. 

Notice the entries for the Plot3D[  ] command. First is  Evaluate[x[t] /. Adv[[1]]] which tells Mathematica what function to plot. Next in the Plot3D  brackets comes 
{x, lower bound for x, upper bound for x}, {t, start time, end time}
which should be compatible with the values computed in the interpolation function Adv. You could use different lower and upper bounds for x, and different start and end times for t in the Plot3D command as compared to those used in the NDSolve command, but those used in the Plot3D command had better be covered by those used in NDSolve, or the interpolated function Adv computed by Mathematica would not give good values for the desired t and x values. The other options in the Plot3D function are of cosmetic nature and could be left away:  AxisLabel -> {names for the axis on your plot} tells Mathematica how to label the axis on the graph.  Following this is  PlotRange -> {0,1} which tells Mathematica  the range for the u[t,x] axis on the plot.  The PlotPoints->100 command tells Mathematica to make a plot with 100 points.  Finally, the Mesh -> False and Boxed -> False options tell Mathematica not to draw the Mesh on the surface and not to draw a box around the graph. 

To get a feel for how this all works for the advection equation, tell Mathematica to graph the solution to  (Subscript[u, t])=-2 (Subscript[u, x]) with the same initial condition, u[0, x] = 1/(1+(x^2)).  Then, use the single advection equation, (Subscript[u, t])=- (Subscript[u, x]) with different initial conditions.  
For example, try the initial conditions u[0,x] = Cos[x], or u[0,x] = x/(1+(x^2)). Then, make up some of your own initial conditions for Mathematica to plot. You might also try (Subscript[u, t]) = (Subscript[u, x]); thus the sign of c has been changed.

The heat equation

HeatEq = 
        NDSolve[{D[u[x, t], t] == D[u[x, t], {x, 2}], 
            u[x, 0] == Exp[-x^2], u[-3, t] == 1/E^9, 
            u[3, t] == 1/E^9}, 
          u, {x, -3, 3}, {t, 0, 3}];
Plot3D[
        Evaluate[u[x, t] /. HeatEq[[1]]], {x, -3, 3}, {t, 0, 3}, 
        AxesLabel -> {"x", "t", "u[t,x]"}, 
        PlotPoints -> 100, Mesh -> False, 
        Boxed -> False];

The wave equation

The NDSolve command can also be used to solve equations which involve more than one time derivative and more than one space derivative. An important example is the wave equation, (Subscript[u, tt])= c (Subscript[u, xx]), where the square root of c has the interpretation as the velocity of the wave. (If c<0, this equation behaves very differently and is called Laplace's equation.) Now, the wave equation takes two derivatives of u, and so one might expect (by analogy with the equation (Subscript[u, tt]) + a (Subscript[u, t]) + b u), that the t=0 values for u and for (Subscript[u, t ]) must be supplied in order to find a unique solution. This is indeed the case. Thus, to have Mathematica plot a solution to   (Subscript[u, tt])= c (Subscript[u, xx]), it is necessary to enter functions for u[0,x] and (Subscript[u, t ])[0,x]. In the following cell, these are taken to by Sin[4 x] and 0:

Waves = 
        NDSolve[{D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}], 
            u[x, 0] == 
              Sin[4*Pi*x], Derivative[0, 1][u][x, 0] == 0, 
            u[0, t] == u[1, t]}, 
          u, {x, 0, 1}, {t, 0, 1}];
Plot3D[Evaluate[u[x, t] /. Waves[[1]]], {x, 0, 1}, {t, 0, 1}, 
        AxesLabel -> {"x", "t", "u[t,x]"}, 
        PlotRange -> {-2, 2}, PlotPoints -> 100, 
        Mesh -> False, Boxed -> False];

There are three formatting changes here that you should note. First, the entry to the D[] operation for the time derivative has D[u[t, x],{t, 2}] which, as you probably guessed, tells Mathematica to take two time derivatives now, not one. Second, since there are now second order partial derivatives with respect to t, I have to specify not just the initial  t=0 condition for u, but also for its time derivative. The condition for u at t = 0 is specified by writing u[0, x] = function of x. We took Sin[4 π x] to be the function of x here.  Meanwhile, the initial condition for the time derivative of u is specified by writing Derivative[u][0, x] = function of x. For this example, we took the constant function 0. Third, we have changed the boundary conditions for u[t,0] and u[t,1]. Instead of setting both to a constant, we told Mathematica  to set u[t,-1] to equal u[t,1]. We could also have set both to zero, or each to a different function of t as long as these functions are compatible with both the initial condition for u(0,x) and  (Subscript[u, t])(0,x). By the way, we could also have set u[t, -1] = u[t, 1] for the heat equation boundary condition. We did so here and not there just to add something new to the story. To get some experience with Mathematica's solutions to the wave equation, have Mathematica plot the solutions to (Subscript[(u ), tt]) = (Subscript[u, xx]) for with 
u[0,x] = Cos[2 π x] and (Subscript[u, t])[0, x] = 0.  
Then have Mathematica plot the solution with 
u[0, x] = Cos[4 π x] and u_t[0, x] = 0. 
Below, we let now Mathematica plot the solution with 
u[0, x] = Exp[ -(x^2)] and (Subscript[u, t])[0,x] = 0.

Waves = 
        NDSolve[{D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}], 
            u[x, 0] == 
              Exp[-x^2], Derivative[0, 1][u][x, 0] == 0, 
            u[-10, t] == u[10, t]}, 
          u, {x, -10, 10}, {t, 0, 10}];
Plot3D[
        Evaluate[u[x, t] /. Waves[[1]]], {x, -10, 10}, {t, 0, 10}, 
        AxesLabel -> {"x", "t", "u[t,x]"}, 
        PlotRange -> {-2, 2}, PlotPoints -> 100, 
        Mesh -> False, Boxed -> False];

Next, have Mathematica plot the solutions for the wave equations (Subscript[u, tt])= 4 (Subscript[u, xx]), with the same initial conditions u[0,x] = Exp[-(x^2)] and (Subscript[u, t])[0, x] = 0. 
By the way, solutions to the wave equation (Subscript[u, tt]) = (Subscript[u, xx ])can also be multiplied by constants and and added together to get new solutions. For example, have Mathematica  plot the solutions to the wave equation with the initial conditions 
u[0, x] = 2 Sin[π x], (Subscript[u, t])[0, x] = 0 
and then have it plot the solution with 
u[0, x] = 4 Sin[ π x], (Subscript[u, t])[0, x] = 0.  
Finally, have it plot the solution with the initial condition 
u[0, x] = Sin[π x] + Sin[4 π x], (Subscript[u, t])[0, x] = 0.

Credit for this assignment

To get credit for this assignment, have Mathematica plot a solution for a wave equation with the schematic form (Subscript[u, tt])=(Subscript[u, xx])+ f(u), where you choose  the function f and the initial conditions. For example, you can take  f = 0 as we did above, or you can take something fancy. In general, take a function which is small because otherwise, you will experience blowup's. You may use the boundary conditions u(x,-10) = u(x,10), or change these as well. Turn in a printout of the plot with the following information:  Your name, the equation used, the initial conditions used and the boundary conditions used.  This should all be placed in your section TF's Mathematics Department mail box on the 3'rd floor of the Science Center next to Room 325 by Wednesday, December 20 in order to receive credit for this assignment. The most elegant plots will be posted on the website.