Tutorial 09 - Symbolic and Numerical Solutions to Differential Equations and their Finite Difference Counterparts

Copyright © 2003 by Bill Titus, Carleton College,
Department of Physics and Astronomy, Northfield, MN  55057-4025
March 1, 2003

Part 1 - Getting Started

Step 1.  Start up Mathematica and open the notebook  Tutorial09.nb.  Remember, this notebook is located in the workshop folder.

Step 2.  OPTIONAL.  Since you're going to edit this notebook, you need to make a copy of this notebook and work with that copy.  To do so, select  File / SaveAs  from the Mathematica menu and save the notebook in the workshop folder under the name  myTutorial09.nb.

Part 2 - Symbolic Solutions to Differential Equations

Topics and Skills

1.  The Mathematica commands  DSolve[]  and  Derivative[]
2.  Replacement rules

Step 1.  Read over the following problem.

Problem 1

In Tutorial 03, you looked at an example where the rate of accumulation,  p,  of carbonate sediments on a reef as a function of the depth  z  below sea level was given by the equation  p[z] = (Subscript[p, 0])(e^(-z/(Subscript[z, 0]))),  where  (Subscript[p, 0])  and  (Subscript[z, 0])  are positive constants.  On a more fundamental level that perhaps speaks directly to the physics and geology behind this model, the expression for  p[z]  is equivalent to saying that the rate of change of  p[z]  with respect to  z  is proportional to  p[z],  where the constant of proportionality is  -(1/(Subscript[z, 0]))  and  p[z]   is  (Subscript[p, 0])  at  z = 0,  i.e.,

				dp[z]/dz  = -p[z]/(Subscript[z, 0]),      p[0] = (Subscript[p, 0]) .
			
Suppose you were given the above differential equation and initial condition as your knowledge about this sedimentation process.  How would you go about solving this differential equation in order to derive the analytic expression for  p[z]?

Step 2.  The Mathematica command that symbolically solves differential equations subject to initial conditions is  DSolve[].   Get Mathematica help on this function.

Comment 2

1.  Here  eqn  is the list of the equations that you want to solve, using the double equal sign notation discussed in Tutorial 06.

Step 3.  To see how to use  DSolve[]  to find the solution to the differential equation in Problem 1, look over the code in the cell below and then activate the cell.

Clear[soln,  p, z, z0, p0]
soln = 
DSolve[
{
p'[
              z] == -p[z]/z0, 
p[0] == 
            p0
}, p[z], z]

Comment 3

1.  Note the form of the output:  it's a list of all the solutions possible, with each solution expressed as a list of replacement rules.  However, in your case, there's only one solution.

Step 4.  You might have guessed that  f'[z]  represents the derivative of  f  with respect to  z.  However if you try to get help on the symbol  ',  Mathematica claims not to know about it.  That's because  f'[z]  is shorthand for the Mathematica function  Derivative[1][f][z].  Get Mathematica help on  Derivative[]  and verify that this is indeed the case.  Make sure you understand the arguments and structure for this function.

Step 5.  Suppose you want to create a function corresponding to the solution obtained in Step 3.  Since the solution is so simple, you'd probably just type in the following code directly:

Clear[p, z, p0, z0]
p[z_, p0_, z0_] = p0*Exp[-z/z0]

However, if the solution were more complicated, it would be a waste of your time to retype an expression that's already typed out on the screen.  Look over the code in the cell below that produces the desired function without explicitly retyping the function.  Activate the cell and then use  ?p  to verify that Mathematica has stored the function correctly.

Clear[p, z, p0, z0]
p[z_, p0_, z0_] = p[z] /. soln[[1]];

Comment 5

1.  Remember, you can interpret the symbol  /.  as  "such that".
2.  The double bracket notation in  soln[[1]]  projects out the desired replacement rule.  Go back to the code in Step 5 and see what happens when the double bracket notation is missing.  After you're seen what happens, restore the code to its origin form.

Step 6.  As practice, create a frame plot for  p[z]  for  (Subscript[p, 0]) = 3.0  and  (Subscript[z, 0])= 20.0.  These values correspond to the parameters given for the sediment problem as stated in Tutorial 03.

Step 7.  Read over the following problem.

Problem 2

Suppose, based on some sound geologic evidence, you felt that a better model for the sediment accumulation problem is given by the differential equation
				
				dp[z]/dz = -(1/(Subscript[z, 0]))cos[2π z/(Subscript[z, 1])]p[z],     p[0] = (Subscript[p, 0]),  
				
where  (Subscript[p, 0]),  (Subscript[z, 0]),  and  (Subscript[z, 1])  are all positive constants.  What's the analytic solution to this problem?  What does   p[z]  look like when the parameters are  (Subscript[p, 0]) = 3.0,  (Subscript[z, 0]) = 20.0,  and  (Subscript[z, 1]) = 30.0.  Remember, copying, pasting and editing can save you time.

Step 8.  Based on Steps 1-6, find the analytic solution for  p[z]  for Problem 2.  Also create a frame plot of  p[z]  using the parameters stated in the problem.

Step 9.  Read over the next problem.

Problem 3

Suppose you have two chemical species, with concentrations  (Subscript[c, 1])[t]   and  (Subscript[c, 2])[t],  that are interacting with one another as a function of the time  t.  A very simple model would be to assume that 

				(Subscript[dc, 1])[t]/dt  = a (Subscript[c, 2])[t]     and   (Subscript[dc, 2])[t]/dt  = b (Subscript[c, 1])[t] ,
		
where  a  and  b  are constants.  Assume that at  t = 0,   (Subscript[c, 1])[0] = 1  and  (Subscript[c, 2])[0] = 0.    Determine analytical expressions for  (Subscript[c, 1])[t]   and   then examine their graphical behavior for different values of  a  and  b.

Step 10.  The code in the cell below determines the analytical solutions to the coupled differential equations given in Problem 3.  Look over the code carefully and then activate the cell.

Clear[soln, c1, c2, t, a, b]
soln = 
DSolve[
{
c1'[
              t] == a*c2[t], 
c2'[t] == 
            b*c1[t], 
c1[0] == 
            1, 
c2[0] == 0
}, {c1[
            t], c2[t]}, t]

Step 11.  The form of the expressions for  (Subscript[c, 1])[t]  and  (Subscript[c, 2])[t]  looks a bit ugly because Mathematica didn't combine  (Sqrt[a]) (Sqrt[b])   into  (Sqrt[a b]).  However, you can tell Mathematica to do so by attaching the replacement rule

				/.(Sqrt[a]) (Sqrt[b]) -> (Sqrt[a*b])
			
at the end of the code in Step 10.  Go back to Step 10 and do this.

Comment 11

1.  Keep in mind that replacement rules allow you to easily transform expressions as well as data.

Step 12.  Use the result for  soln  in Step 10 to define two functions,  c1[t,a,b]  and  c2[t,a,b].  To do so, you'll need to generalize what you did in Step 5 above.  Get Mathematica help on your functions as a quick check to see that you've defined the functions correctly.

Comment 12

1.  If  (Subscript[c, 1])[t]   and  (Subscript[c, 2])[t]  are really concentrations, there's something not physical about your plot.  Do you know what it is?

Step 13.  Graph  (Subscript[c, 1])[t]  and  (Subscript[c, 2])[t]  on the same plot for  a  = 1.0  and  b  = 1.0.  Adjust the range of  t  so you get a good sense for what the concentrations are doing as a function of time.  Plot each curve in a different color so you can tell the curves apart.

Step 14.  Use  ParametricPlot[]  to create a graph of  (Subscript[c, 1])[t]  versus  (Subscript[c, 2])[t]  for  a = 1.0  and  b  = 1.0.  This type of plot is also useful in determining how the system behaves.

Comment 14

1.  I used the option  AspectRatio->Automatic  in order to make the scale for  (Subscript[c, 1])[t]  and  (Subscript[c, 2])[t]  visually the same.

Step 15.  To get a good sense of the behavior of these concentrations as  a  and  b  change, repeat Steps 13 and 14 for a number of different values of  a  and  b.  Remember, there's no reason that  a  and/or  b  couldn't be negative or even zero.

Part 3 - Numerical Solutions to Differential Equations

Topics and Skills

1.  The Mathematica commands  NDSolve[]  and  InterpolatingFunction[]
2.  The Mathematica commands  GraphicsArray[]  and  Partition[]

Step 1.  Read over the following problem.

Problem 4

Suppose you're studying the concentration of a certain chemical species that is subject to interactions from both inside and outside the system.  You have a sense that you can model the species using

  				             dc[t]/dt  =  a c[t] (1 - c[t] / b),     c[0] = (Subscript[c, 0]).  
   
Here  a,  b,  and  (Subscript[c, 0])  are positive constants but constrained so that  c[t],  the concentration, always lie between  0  and  1.

Step 2.  Look over the code below.  This code uses  DSolve[]  to find the analytic solution to the differential equation in Problem 4 and then uses that result to define the function  cExact[t,a,b,c0].   Activate the cell and see what the solution looks like.

Clear[soln,  c, t, a, b, c0, cExact]
soln = 
DSolve[
{
c'[
              t] == 
            a*c[t] (1 - c[t]/b), 
c[0] == 
            c0
}, c[t], 
        t]

cExact[t_, a_, b_, c0_] = c[t] /. soln[[1]]

Comment 2

1.  The warning message is just telling you that Mathematica isn't quite certain that it has found all the solutions.  However, in this case it actually did.

Step 3.  To get a sense for what the solutions look like for  0  <  t  <  20,  create a table of plots for  (Subscript[c, 0]) = 1/10  with  b = 3/4  and  a  going from  0  to  5  in steps of  1/2  by activating the code below.  Be sure you understand what this code is doing, especially the options in  Plot[].

Clear[c0, a, b, tMax, da, plot]
c0 = 1/10; b = 3/4; tMax = 20; da = 1/2;
Table[
        plot[a] = 
          Plot[cExact[t, a, b, c0], {t, 0, tMax}, 
            PlotStyle -> RGBColor[1, 0, 0], Frame -> True, 
            PlotRange -> {0, 1.02}, 
            PlotLabel -> 
              StringForm["a =``", PaddedForm[N[a], {2, 2}]]], {a, 0, 5, 
          da}];

Comment 3

1.  Each plot has been given its own name so you can refer to the plots in the next step.

Step 4.  It's a bit hard to page through all the solutions to see how the plot behaves as a function of  a.  One way to compare the plots is to create a single frame that contains all of them.  You can do so by using two new Mathematica functions,  GraphicsArray[]  and  Partition[].  Get Mathematica help on these functions.  Then look over the code below and activate the code, but only after you understand what it's doing.

Show[GraphicsArray[Partition[Table[plot[a], {a, 0, 4, da}], 2]], 
        Frame -> True, 
        PlotLabel -> 
          StringForm["c[t] with  b = ``", b]];

Step 5.  There will be times when  DSolve[] will be unable to solve a differential equation analytically, i.e., come up with a closed, analytic solution in an algebraic form.  However, any differential equation can be solved numerically, although to do so requires explicit numerical values for the constant parameters involved in the problem.  Find the name of the Mathematica function that numerically solves differential equations and then get information about how to use this function.

Comment 5

1.  I came up with the name of the function by using Mathematica help in the form  ?*DSolve*.  Remember,  *   serves as a wildcard.
2.  It shouldn't come as a shock to you that the desired function is  NDSolve[].  Note that the argument for  NDSolve[] includes the range of the independent variable over which you want to find the numerical solution.

Step 6.  Look over the code below that numerically solves the differential equation in Problem 4 when  a = 0.7,  b = 0.75,   (Subscript[c, 0]) = 0.1  and  t  goes from  0.0  and  10.0.  Activate the code and try to figure out the meaning of the solution that's displayed.

Clear[a, b, c0, tMax, soln, c, t]
a  =  0.7; b = 0.75; c0 = 0.1; 
    tMax = 10.0;
soln = 
NDSolve[
{
c'[t] == 
            a*c[t] (1 - c[t]/b), 
c[0] == 
            c0
}, c[t], {t, 0, tMax}]

Comment 6

1.  As you might gather, an  InterpolatingFunction[] is one that gives approximate values to the actual function.  It does so by storing values of the function at a finite set of grid points and then using an interpolating algorithm to evaluate the function at intermediate points between the grid values.
2.  Once a function has been defined via InterpolatingFunction[],  Mathematica treats it as a regular function and you can plot, integrate, and even differentiate it.

Step 7.  Since a function produced by  InterpolatingFunction[]  is a legitimate function, you can define the approximate function for the concentration in Problem 4 just as you did for the exact function in Step 2.  Look over the code below and then activate it to define the approximate function  cApprox[].

Clear[cApprox, t]
cApprox[t_] = c[t] /. soln[[1]]

Step 8.  Graph cExact[t, a, b, c0]  and   cApprox[t]  on the same plot for the parameters stated in Step 6.    Have each curve in a different color so you can tell the curves apart.  This composite plot should give you some sense of how the actual and approximate solutions differ.

Comment 8

1.  When I created the plot, I could see no visual difference between the exact and the approximate solutions.

Step 9.  Sometime you can get a sense of a dynamical system, such as the one above, by creating a "phase plot" for the system, i.e., a plot of  c[t]  versus  dc[t]/dt.  Produce such a plot for our system with the parameters stated in Step 6, using the approximate function from Step 5 for the concentration.  You'll first need to define a function corresponding to the derivative of  c[t]  with respect to  t  and then use  ParametricPlot[] to create the phase plot.

Comment 9

1.  Note that Mathematica has taken the derivative of your approximate function without any complaints.
2.  Your curve should look parabolic.

Part 4 - Finite Difference Equations

Problem 5

Suppose it's 1960 and you have to numerically solve the equation  dc[t]/dt  =  g[c[t]]  over the time interval  {0, (Subscript[t, max])}  subject to the initial condition  c[0] = (Subscript[c, 0]),  where  g[c]  is a known function of  c.  For example,  for exponential growth,  g[c]  =  r c  and for the "logistic equation" in Part 3,   g[c]  =  a c (1 - c / b).   What you might do is to break the time interval  {0, (Subscript[t, max])}  into  N  small intervals  {(Subscript[t, n]), (Subscript[t, n + 1])},  each of the same size,  δt.  You would then find  c[(Subscript[t, n])]  by approximating   dc[t]/dt  with  (c[(Subscript[t, i + 1])] -c[(Subscript[t, i])]) / δt.  This approximation, coming from the definition of the derivative, leads to the iterative, difference equation

			    		     c[(Subscript[t, n + 1])]  =  (1 + a δt) c[(Subscript[t, n])](1 - a δt c[(Subscript[t, n])] / b)
			
from which you can find  c[(Subscript[t, n + 1])]  knowing  c[(Subscript[t, n])].  If you let 

					 d[(Subscript[t, n])] = a δt c[(Subscript[t, n])] / b    and    r = (1 + a δt) ,
			 
the difference equation simplifies to the form

					d[(Subscript[t, n + 1])]  =  r d[(Subscript[t, n])](1 - d[(Subscript[t, n])]) .

The symbol  t  is now superfluous and you end up with
 
 					d[n+1] = r d[n] (1- d[n]), 
 	
with  r  and  d[0]  known.  Create some code that allows you to solve this difference equation for different values of  r  and  d[0],  and for any number of total interval steps  (Subscript[n, Max]).

Topics and Skills

1.  One dimensional array variables
2.  Animation using  (apple)(y)

Step 1.  Define the function  g[d, r] = r d (1 - d)  and plot the function for  0.0  < d  < 1.0.  Use your plot to convince yourself that  r  must lie between  0.0  and  4.0  in order for  d  to lie between  0.0  and  1.0.

Step 2.  You're going to represent  d[n]   in Mathematica by a one-dimensional array  d[n].  In actual fact,  d[n]  is a function, but one that only takes on values when  n  is an integer.  The following code is one way in which to solve the iterative equation stated in Problem 5 for given values of  r  and  (Subscript[d, 0]) for the first  (Subscript[n, Max]) values.  Look over the code below and active the cell containing the code after you understand what the code is doing.

Clear[r, d, nMax]
r = 0.7; d[0]  =  0.1; nMax = 10;
Table[d[n] = r*d[n - 1]*(1 - d[n - 1]), {n, 1, nMax}]

Step 3.  You'd like to plot the values contained in the table and the value of  d[0]  as a function of  n.  To do so, use  ListPlot[] to plot a list whose elements are of the form  {n, d[n]}.  Use the  Table[]  command to generate the list, call your list  dList,  and plot the results.

Step 4.  Look over the code below that combines Step 1- 3 in the same cell and then activate the cell.

Clear[r, d, nMax, ptSize, dList]
r = 0.7; d[0]  =  0.1; nMax = 10; 
    ptSize = 0.01;
Table[
        d[n] = r*d[n - 1]*(1 - d[n - 1]), {n, 1, 
          nMax}];

dList = Table[{n, d[n]}, {n, 0, nMax}];
ListPlot[dList, Frame -> True, 
        PlotStyle -> {RGBColor[0, 0, 1], PointSize[ptSize]}, 
        PlotLabel -> "d[n]"];

Step 5.  If you want to produce a number of such plots for different values of  r,  (Subscript[d, 0]),  and  (Subscript[n, Max]),  it would be nice to have a single function that did it all at once.  In Tutorial 06, you learned about the Mathematica function   Module[] that allowed you to wrap a section of code all together to produce a single function call.  The code below does this.  Look over this code and active it.  Try out your function for the parameters given in Step 2.

Clear[makeOrbit, r, d0, nMax, ptSize]
makeOrbit[r_, d0_, nMax_, 
        ptSize_] := 
Module[
{d, n, 
          dList}, 

d[0] = 
          d0; 
Table[
          d[n] = r*d[n - 1]*(1 - d[n - 1]), {n, 1, 
            nMax}]; 
dList = 
          Table[{n, d[n]}, {n, 0, nMax}]; 
ListPlot[dList, 
          Frame -> True, 
          PlotStyle -> {RGBColor[1, 0, 0], PointSize[ptSize]}, 
          PlotLabel -> "c[n]", 
          PlotRange -> {0, 1}]
]

Step 6.  Use your function  makeOrbit[] and Mathematica's  Table[]  command to produce a series of plots corresponding to  (Subscript[c, 0]) = 0.1,  (Subscript[n, Max]) = 100,  and  r  going from  0.0  to  4.0  in steps of  0.1.

Step 7.  It's very difficult to page through the plots and see the general trend as  r  varies.  However, Mathematica can animate the process by running through all the frames consecutively to produce a movie.  To do so, double click on the cell bracket surrounding all the frames, after which only a single frame will show.  Then highlight the cell bracket and press  (apple)(y)  or go to the Mathematica menu and select  Cell / Animate Selected Graphics.  Now, for a truly exciting experience, go back to Step 6 and change the increment in  r  to  0.01  instead of  0.1.  Produce the frames (this will take a while since there are lots of plots to produce) and animate the results.

Comment 7

1.  When the move runs through once, it recycles back to the beginning and starts again.
2.  While the animation is running, there are some controls in the lower left hand corner that you can play with that affect the direction the movie plays and its speed.

Part 5 - Save Your Notebook (OPTIONAL)

Step 1.  Remove all output cells from this notebook, and to do one final save of your notebook before you quit Mathematica.