Math 21a: Multivariable Calculus, Fall 2004
Mathematica Assignment 2

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 Mathematica assignment about the special notation and formatting rules for Mathematica. In any event, here is a brief summary of some of the more useful 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 that involve a function, f(t), of a single variable t. Of these "Ordinary Differential Equations" one of the simplest is just f'(t) = g(f) where g is a function. For example, f' = f,  or  f' = f (1 - f). Equations of the form  f' = k f with k a constant are called exponential growth equations.  To see what a solution looks like with k = 2, have Mathematica evaluate the following:

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

Explanation: The special Mathematica function NDSolve solves differential equations for functions f(t). It has the form "Name = NDSolve[{equation to be solved, initial value for f}, f, {t, start time, end time}]".  The notation is a bit tricky - here the double equal sign == means "identical to" and it shouldn't be confused with the equal sign =, which is used to assign values. For example "t=5" will set the variable t to be 5. Then if someone later evaluates the statement "t==5" then Mathematica should respond with the answer "True" as Mathematica has simply checked to see if the variable t is now equal to 5. 
Applying NDSolve tells Mathematica to solve the equation for f(t) from the start time to the end time with the given starting value for f. The result "Name" (in the above example Name is simply equal to the label "New") is a rule that refers to a function derived by interpolation. We get a full fledged function of the variable t from this rule by writing "f[t] /. Name[[1]]" . This function is plotted with the command Plot[f[t] /. Name[[1]], {t, start time, end time}];

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

To gain practice, change the starting value for the Logistic equation to f[0] = -1 and graph the result.  Also, try f[0] = 0.5 and f[0] = 2.  Do you remember enough from Harvard's Math 1b (or an equivalent introduction to differential equations) to explain the behavior of f[t] for these various starting values? If not, don't be too concerned - the main results for the 21a class will come in the study of the PDEs in the following sections.  For a second example, change the equation in the NDSolve entry to 
f'[t]=Sin[f[t]] and plot solutions with f[0]=π/4, 2π/4,5π/4 and 9π/4.

Differential equations of second order

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

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

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

Comments: 
With regard to the formatting here, note that for this second order differential equation two initial conditions must be specified to nail down a specific solution function. This is a pretty general rule:  For an nth order differential equation you must specify the initial values of f 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 often isn't any solution that fits all the given constraints. By copying the preceding cell and changing the equation appropriately, have Mathematica plot solutions to f'' +  f' + f = 0, f'' + 2 f' + f = 0, f'' + 3 f' + f = 0, and f'' - f' + f = 0 all using the same initial conditions.  
You might notice the different behavior; the first 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 f'' + a f' + b = 0 and your own initial conditions and have Mathematica plot the solutions.

The Advection Equation - a first order PDE

Mathematica can also plot solutions to partial differential equations (PDEs) for two variable functions u(x, t) which depend a spacial coordinate x and on time t. One of the simplest such PDE is the so-called "Advection Equation" (Subscript[u, t]) = -c (Subscript[u, x]), where c is a constant. Solutions to this equation describe travelling waves (as opposed to solutions to the Wave Equation PDE that describes standing waves, such as those of a vibrating string). Given any initial function g(x), there is a unique solution to the advection equation, u(x,t), with the property that u(x,0) = g(x).  To see an example solution, have Mathematica  evaluate the following (don't sweat any warning messages that come out - Mathematica can be pretty fickle!):

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[x,t]"}, 
        PlotRange -> {0, 1}, PlotPoints -> 100, Mesh -> False, 
        Boxed -> False];

Comments:
First off, try to get a feel for how to "read" this 3D graph.  The front cross section of the graph is the plot of the wave at time t = 0, just a single high peaked wave located with maximum height at x = 0.  As time passes by, i.e. as t goes from 0 up past 8 seconds (if t is taken to represent seconds), the wave slides, i.e. travels, off to the right (in the positive x direction).  In the solution to the advection equation note that the wave doesn't lose height as it travels, whereas in reality a wave might slowly decrease in height as it spreads out due to other physical forces.
Now, next here are some comments about the formatting of the expressions in the preceding cell.  There are a lot of comments, so bear with us as you read through this - you can always come back and refer to this later if you get confused in the PDE examples that follow this one.  First, 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, variable for the partial derivative], so D[u[x,t], t] tells Mathematica to take the partial derivative of u[x,t] with respect to t, while D[u[x,t],x] tells Mathematica to take the partial derivative with respect to x. 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[x,0], 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[-20,t] to be any function of t which is consistent with the choice of u[-20,0]. We choose the constant function u[-20,t] = u[-20,0] for simplicity. 

Notice the entries for the Plot3D[  ] command. First is  Evaluate[u[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 x and t values. The other options in the Plot3D function are of cosmetic nature and could be left away:  AxesLabel -> {names for the axes on your plot} tells Mathematica how to label the axes on the graph.  Following this is PlotRange -> {0,1} which tells Mathematica  the range for the u[x,t] 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[x,0] = 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[x,0] = Cos[x], or u[x,0] = 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[x,t]"}, 
        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^2 (Subscript[u, xx]), where in the solution function the constant c is related to the speed of the wave. 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 both the function u[x,0] and for the partial derivative (Subscript[u, t ])[x, 0] must be supplied in order to determine 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[x,0] and (Subscript[u, t ])[x,0]. In the following cell, these are taken to by Sin[4 x] and 0 (again, feel free to laugh at any warning signs that come out as the solution is being graphed):

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[x,t]"}, 
        PlotRange -> {-2, 2}, PlotPoints -> 100, 
        Mesh -> False, Boxed -> False];

Again, the interpretation of this graph is that the front cross section corresponds to the vibrating string at the start of the time being observed.  As time passes, look at the cross sections going further back on the graph - each cross section corresponds to the status of the vibrating string at a particular time.  Note that for each cross section, the end points of the string end up constant, with  u[x,t] equaling 0 at the two endpoints, x = 0 and x = 1 (corresponding to the back left and front right of the graph), as if the string is tied down at these two points.
There are three formatting changes here that you should note. First, the entry to the D[] operation for the time derivative has D[u[x,t],{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[x,0] = 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][x,0] = function of x. For this example, we took the constant function 0. Third, we have changed the boundary conditions for u[0,t] and u[1,t]. Instead of setting both to a constant, we told Mathematica to set u[0,t] to equal u[1,t]. 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(x,0) and  (Subscript[u, t])(x,0). 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[x,0] = Cos[2 π x] and (Subscript[u, t])[x,0] = 0.  
Then have Mathematica plot the solution with 
u[x,0] = Cos[4 π x] and (Subscript[u, t])[x,0] = 0.

Next, have Mathematica plot the solutions for the wave equations (Subscript[u, tt])= 4 (Subscript[u, xx]), with the same initial conditions u[x,0] = Exp[-(x^2)] and (Subscript[u, t])[x,0] = 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[x,0] = 2 Sin[π x], (Subscript[u, t])[x,0] = 0 
and then have it plot the solution with 
u[x,0] = 4 Sin[ π x], (Subscript[u, t])[x,0] = 0.  
Finally, have it plot the solution with the initial condition 
u[x,0] = Sin[π x] + Sin[4 π x], (Subscript[u, t])[x,0] = 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 just take f = 0 as we did above, or you can take something fancy. In general, take a function which is small as otherwise, you will experience blowups as the PDE is being resolved. You may use the same boundary conditions as in the example above 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 be dropped off in your CA's mail box (these are all located in the black filing cabinet on the 3rd floor of the Science Center to the left of room 310) by 5pm on Friday, November 19th in order to receive credit for this assignment.