WOLFRAM 
   TECHNOLOGY 
         CONFERENCE 2004

Computational Science with Mathematica

12 Easy Pieces: Essential Operations in Computational Science

James L. Noyes
Wittenberg University

Introduction:

Computational Science (COSC) is a scientific investigation of physical processes on the computer by the use of computational models and methods.  Hypotheses are often formulated using mathematical models that can then be used to compute quantities of scientific interest and to gain new insights.  COSC is now on par with the development of scientific theory and the use of experimentation in order to understand the real world.  Mathematica is an ideal tool to use in this endeavor and also to teach the essential types of operations [7].  These types of operations may be placed into three categories: Evaluation, Simulation, and Optimization.  Within each category, fundamental operations will be described.  One or more examples will be presented for each operation.  This material is currently used to teach an undergraduate course in COSC.

Evaluation: This is the most basic category since the operations in this category are used in both simulation and optimization.

1. Analytic (Symbolic) vs. Numeric Computation.  The primary advantage of analytic ("closed form") solution methods is that one may examine this form (expressions, equations, and formulas) and investigate the behavior by changing parameters and variables and observing the consequences.  In general, one should try to obtain a simplified analytic solution whenever possible (using numerical solutions to aid the validation process) [5].  Unfortunately, for many problems, no elementary analytic solution exists, or it is impossible to obtain such a solution.  Hence a numeric solution must be found.  While numeric solutions are frequently possible, they are sometimes subject to truncation error, subject to numeric error, and often must be re-computed for changes in variables and parameters.  Examples:

A. Using Hooke's Law to Determine Work Formula for Stretching a Spring

Clear[h, k, x, a, b]
h[x_] := k x

Simplify[(∫_a^b) h[x] \[DifferentialD]x]

B. Integrating the Standard Normal Probability Density Function

Clear[ϕ, θ, 
      z, τ, μ, σ]
ϕ[
        z_] := (1/(Sqrt[2  π])) (E^(-((z^2)/2)))

Plot[ϕ[z], {z, -3, 3}];

(∫_(-∞)^∞) 
        z ϕ[z] \[DifferentialD]z

(∫_(-∞)^∞) (z^2) ϕ[
          z] \[DifferentialD]z

θ[
        z_] := (∫_(-∞)^z) ϕ[τ] \[DifferentialD]τ

θ[((x - μ)/σ)]

C. Integrating an Arbitrary Function (Simple Trapezoidal Rule)

Clear[f, x, Trap]
f[x_] := Sin[x (E^(0.25  (x^2)))]

(∫_0^π) f[x] \[DifferentialD]x

Trap[a_, b_, 
        n_] := ((b - a)/(2  n)) (f[a] + 
              f[b]) + ((b - a)/n) (∑_(j = 1)^(n - 1))f[
              a + (((b - a)/n)) j]

Trap[0, π, 5000]
N[(∫_0^π) f[x] \[DifferentialD]x]

2. Data Input/Output and Visualization.  Computational Science must deal with both small and large amounts of file data to be read or written.  Visualization can be critical to the interpretation of data, especially large amounts of it.  Often this data must be reformatted or filtered in some fashion before it can be used.  Many times a problem or function is so complex that it is required to have some type of graphical representation in order to understand the effect of variable and parameter changes.  Examples:

A. Input and Display of a Small Data File

Clear[inData]
inData = 
      ReadList["c:\\Witt\\WolframConf\\infile.txt", {Word,  Number, 
          Number, Number}]
inData // TableForm
inData[[2]]

B. Create and Output a Data File of Elevations (ft.) Over a Given Area (mi.)

Clear[lcOutData, myStream]
SeedRandom[1];
lcOutData = 
          8000.0 + 
            1000.0  Table[
                Sin[x y] + Random[Real, {-0.4, 0.4}], {x, 
                  0, ((3  π)/2), (π/15)}, {y, 
                  0, ((3  π)/2), (π/15)}];;
myStream =  
        OpenWrite["c:\\Witt\\WolframConf\\lakecity.txt"];
Write[myStream, lcOutData];
Close[myStream];

C. Input and Visualize a File of Elevations (ft.) Over a Given Area (mi.)

Clear[lcInData]
lcInData = 
        ReadList["c:\\Witt\\WolframConf\\lakecity.txt"];
lcInData[[1]] // TableForm
ListPlot3D[lcInData[[1]]];

D. Visualizing a Parameterized Bivariate Function

Clear[ρ]
ρ = 
      0.0;                                                                                (* 
      Try ρ = 0.0,  ρ = 0.1,  ρ = 
        0.5 *) 
Plot3D[(E^Sin)[ρ x y], {x, 0, 2  π}, {y, -π, π}];

3. Physical Units and Conversions.  Most realistic problems involve physical units (mass, length, time, current, temperature, etc.).  An incorrect specification of units in a model is just as bad as an incorrect numerical quantity.  Dimensional analysis can be used, whereby units are treated as algebraic quantities so that quantitative scientific relationships may be converted, simplified, and verified.  Example:

A. Compute the Speed of Light in Furlongs per Fortnight

(299792458 
           (m/s)) (((100  in)/(2.54  m))) ((yd/(36  in))) ((fur/(220  yd))) (((3600  s)/hr)) (((24  hr)/day)) (((14  day)/ftn))

4. Linear vs. Nonlinear Systems.  Linear systems of the form Ax=b are ubiquitous in modeling situations and their related solution methods.  Not only do they arise directly, but they frequently arise as sub-problems of a more complex problem (e.g., optimization and Markov chains).  Error formulas may be used to estimate or bound the relative error in the solution, ||Δx||/||x||, by computing norms and a condition number.  This is especially important if an ill-conditioned system is suspected.  Eigenvalues and eigenvectors that satisfy Au=λu, can also be important.  Linear systems can often be applied to approximating and solving nonlinear systems (e.g., using a multi-variable Taylor Series approximation).  Examples:

A. Solving an Ill-Conditioned Vandermonde Linear System

Clear[n, e, V, x, b, y, W, z]
n = 
      6;                                                                                        (* 
      Try n = 6 and n = 15 *) 
e = 
      1.0;
V = 
        Table[(j^(i - 1)), {i, 1, n}, {j, 1, 
            n}];
x = Table[100, {i, 1, n}];
If[n ≤ 6, V // MatrixForm]
b = V . x;
(y = LinearSolve[V, b]) // MatrixForm
W = Inverse[V];
(z = LinearSolve[(e V), b]) // 
      MatrixForm
Print["Abs.Err(Exact Arith.): ", 
      Norm[(y - x), ∞]]
Print["Rel.Err(Exact Arith.): ", 
      Norm[(y - x), ∞]/
        Norm[x, ∞]]
Print["Abs.Err(FP Arith.): ", 
      Norm[(z - x), ∞]]
Print["Rel.Err(FP Arith.): ", 
      Norm[(z - x), ∞]/
        Norm[x, ∞]]
Print["Matrix Cond.No.: ", (Norm[V, ∞] 
            Norm[W, ∞]) // N]

B. Linearizing a Nonlinear Bivariate Function About a Given Point

Clear[n, x, y, δ, l]
(Subscript[x, 0]) =. ; (Subscript[y, 0]) =. ; (Subscript[l, 0]) =. ; (Subscript[l, x]) =. ; 
    (Subscript[l, y]) =. ;
n[x_, y_] := (E^(-x y)) 
        Sin[x y]
(Subscript[x, 0]) = 1. ;
(Subscript[y, 0]) = 2. ;
δ = 
      3.0;                                                               (* 
      Try δ = 0.5,  δ = 1.0,  δ = 
        3.0 *) 
(Subscript[l, 0]) = 
      n[x , y] /. {x -> (Subscript[x, 0]), 
          y -> (Subscript[y, 0])};
(Subscript[l, x]) = (∂_x) n[x, y] /. {x -> (Subscript[x, 0]), 
            y -> (Subscript[y, 0])};
(Subscript[l, y]) = (∂_y) n[x, y] /. {x -> (Subscript[x, 0]), 
            y -> (Subscript[y, 0])};
l[x_, y_] = 
      Simplify[(Subscript[l, 0]) + (Subscript[l, x]) (x - (Subscript[x, 0])) + (Subscript[l, y]) (y - 
                (Subscript[y, 0]))]
Plot3D[
        Abs[n[x, y] - l[x, y]], {x, (Subscript[x, 0]) - δ, (Subscript[x, 0]) + δ}, {y, 
          (Subscript[y, 0]) - δ, (Subscript[y, 0]) + δ}];

5. Data and Function Approximation.  It is not uncommon to require that a data set be replaced by a single function (e.g., using the least-squares criterion or spline approximation) or that sometimes a complex function be approximated by a simpler function (e.g., a Taylor or Chebyshev approximation).  In both cases, there is an art to selecting the appropriate functional form, and visualization is usually necessary to help do this.  Once the functional form has been chosen, the determination and verification of its optimal parameters is often straightforward.  Examples:

A. Determine the Least-Squares Fit of a Linear Parameter Model

Clear[xydata, xyplot, g, x, gplot]
xydata = {{0.0, 7.3}, {0.5, 9.2}, {1.0, 11.0}, {1.5, 12.3}, {2.0, 
            13.1}, {2.5, 13.4}, {3.0, 13.6}, {3.5, 13.7}, {4.0, 14.2}, {4.5, 
            15.4}, {5.0, 17.3}};
xyplot = ListPlot[xydata];

g[x_] = Fit[xydata, {Sin[x], (E^(x/5))}, {x}]

gplot = Plot[g[x], {x, 0, 5}];
Show[xyplot, gplot];

(∑_(i = 1)^(Length[xydata]))(((g[xydata[[i, 1]]] - xydata[[i, 2]]))^2)

B. Determine the Least-Squares Fit of a Nonlinear Parameter Model [4]

Clear[pvdata, pvplot, params, v, p, a, b, vplot, 
      pred]
pvdata = {{341948. , 4.81}, {1092759. , 5.88}, {5491. , 
            3.31}, {49375. , 4.90}, {1340000. , 5.62}, {365. , 
            2.76}, {2500. , 2.27}, {78200. , 3.85}, {867023. , 
            5.21}, {14000. , 3.70}, {23700. , 3.27}, {70700. , 
            4.31}, {304500. , 4.42}, {138000. , 4.39}, {2602000. , 
            5.05}};
Print ["P vs. V Data: ",  
        TableForm[pvdata]];
pvplot = ListPlot[pvdata];

params = FindFit[pvdata, a (p^b), {a, b}, p]

v[p_] := a (p^b) /. params
v[p]

vplot = Plot[v[p], {p, 0, 2602000}];
Show[pvplot, vplot];

(∑_(i = 1)^(Length[pvdata]))(((v[pvdata[[i, 1]]] - pvdata[[i, 2]]))^2)

pred = 
        Table[{pvdata[[i, 2]], v[pvdata[[i, 1]]]}, {i, 1, 
            Length[pvdata]}];
Print["Observed vs. Predicted V Data:", 
        TableForm[pred]];

6. Computer Arithmetic.  Although integer arithmetic is exact in Mathematica, neither integer nor "real" (floating-point) computer arithmetic is exactly like its mathematical counterpart when using compiled languages such as C/C++ and Fortran.  FP arithmetic is especially vulnerable to roundoff error and it can be important to assess the effect of roundoff error in modeling situations (especially if single precision is to be used).  Example:

A. Evaluate a(b - c) in Simulated 4-Digit FP Arithmetic With 5 Roundoffs

Clear[fl]
<< NumericalMath`ComputerArithmetic`

SetArithmetic[4, 10, RoundingRule -> RoundToInfinity, 
        ExponentRange -> {-50, 50}, 
        MixedMode -> True];

fl[x_] := ComputerNumber[x]

fl[fl[1234.0]*fl[fl[567.84] - fl[567.75]]]
fl[567.84]
fl[567.75]

7. Uncertainty and Sensitivity in Computation.  It is not uncommon to have measured or computed quantities to be used in a formula, where these values are not known exactly, but may be approximately known or are known to be within a given domain.  In the first case, one may wish to know the first-order error estimates, so partial derivatives can be used in a first-order multivariate Taylor Series.  In the second case, one may want to know the range of the answer, so interval arithmetic can be used.  Examples:

A. First-Order Cylinder Volume Sensitivity

Clear[V, R, 
      H, ΔR, ΔH]
V = π (R^2) H;
((∂_R) V) ΔR  +  ((∂_H) 
            V) ΔH
Collect[%, π] /. {R -> 2, H -> 6}

B. Interval Arithmetic to Assess Equivalent Resistance Formulas [9]

(* (((Subscript[R, 1]) (Subscript[R, 2]))/((Subscript[R, 1]) + (Subscript[R, 2]))) with (Subscript[R, 1]) 
        in [(4/5), 1] and (Subscript[R, 2]) in [2, 
          (5/2)] *) 
((Interval[{4/5, 1}] Interval[{2, 5/2}])/(Interval[{4/5, 1}] + Interval[{2, 5/2}]))

    %[[1, 2]] - %[[1, 1]] // N

(* 
      (1/((1/(Subscript[R, 1])) + (1/(Subscript[R, 2])))) with (Subscript[R, 1]) 
        in [(4/5), 1] and (Subscript[R, 2]) in [2, 
          (5/2)] *) 
(1/((1/Interval)[{4/5, 1}] + (1/Interval)[{2, 5/2}]))

    %[[1, 2]] - %[[1, 1]] // N

Simulation: This is the most general category where models must often be implemented by writing application-specific programs.

8. Discrete Equation Models.  These describe the relationship between the state of some system or process between discrete intervals.  Difference equation models are among the simplest models to implement in code and to study.  They can be deterministic or stochastic.  They tend to be straightforward (when special coding techniques are not necessary) and intuitive.  Examples:

A. Simple Deterministic Biological Population Dynamics

Clear[old, maxperiod, growthrate, popseq, 
      new]
old = 10;
maxperiod = 50;
growthrate = 0.03;
popseq = {old};
Do[ (new = old + growthrate old;  
          popseq = Append[popseq, new];  
          old = new), 
{maxperiod}];
ListPlot[popseq, PlotStyle -> PointSize[0.02]];

9. Differential and Integral Equation Models.  Differential equations describe the relationship between the state of something and its instantaneous rate of change with respect to some variable.  These causal models are often involved, possibly requiring systems of ordinary or partial differential equations and boundary conditions to describe the physical problem.  Analytic solutions are desirable, but often the solutions must be numeric.  Integral equation models describe the relationship between the state of something and the accumulative changes in all of its previous states.  These models can also be quite involved, requiring the solution to single and multiple integrals.  Examples:

A. Investigate Newton's Law of Cooling

According to Newton's law of cooling, the time rate of change in temperature, du/dt, is proportional to the temperature difference between the current temperature u(t) and the given lower steady-state temperature (Subscript[u, s]).  For a decrease, the given constant of proportionality is -k where k≥0.  This first order ODE is given by the following relationship: du/dt = -k(u -  (Subscript[u, s])), with the given initial condition of u(t) = u(0) = (Subscript[u, 0]) at time t=0.

1. Solve this ODE for u(t) in simplest form when (Subscript[u, 0]) = 70 Fahrenheit and k = 0.462 and plot u(t) for 0 ≤t≤3 minutes.  2. Determine the temperature at t = 3 minutes.  3. Determine the time in minutes when the temperature is 15 degrees Fahrenheit

Clear[ans, u, t, ua]
ans = 
      DSolve[{u'[t] ==  -0.462 (u[t] - 10), 
          u[0] == 70}, u[t], t]
ua[t_] = 
      Simplify[ans[[1, 1]][[2]]]
Plot[ua[t], {t, 0, 3}];

Print["u[3] = ",  ua[3]];

Solve[ua[t] == 15, {t}]

B. Approximate the Behavior of an ODE/IVP System [1]

(dx/dt)  =  10 (y  -  x);               (Subscript[x, 0]) = 
      1
(dy/dt)  =  28  x  -  y  -  x z;      (Subscript[y, 0]) = 
      5
(dz/dt)  =  
      x y  -  (8/3) z;                (Subscript[z, 0]) = 
      10

Clear[x, y, z, xn, yn, zn, t]
(Subscript[x, 0]) =. ; (Subscript[y, 0]) =. ; (Subscript[z, 0]) =. ;

(Subscript[x, 0]) = 1 ; (Subscript[y, 0]) = 5 ; (Subscript[z, 0]) = 
        10 ;

ans = 
        NDSolve[{x'[t]  ==   10 (y[t]  -  x[t]), 
            x[0] ==  (Subscript[x, 0]), 
y'[t]  ==  
              28  x[t]  -  y[t]  -  x[t] z[t], 
            y[0] ==  (Subscript[y, 0]), 
z'[t]  ==  
              x[t] y[t]  -  (8/3) z[t], z[0] ==  (Subscript[z, 0])}, {x[t], 
            y[t], z[t]}, {t, 0, 5}];

xn[t_] = ans[[1, 1, 2]]
yn[t_] = ans[[1, 2, 2]]
zn[t_] = ans[[1, 3, 2]]
Plot[{xn[t], yn[t], zn[t]}, {t, 0, 1}];

Print["x(0.5)= ", xn[0.5]]
Print["y(0.5)= ", yn[0.5]]
Print["z(0.5)= ", zn[0.5]]

Plot[{xn[t], yn[t], zn[t]}, {t, 0, 5}];

10. Monte Carlo Models.  These models can be used when analytic and numeric solutions are both computationally impractical or when only an approximate solution is needed.  Here one simulates an experiment by employing pseudo-random numbers, which are used to generate samples.  It is crucial that a good PRN generator is chosen, that is compatible with the computer arithmetic being used.  Typically, several simulations are done for each estimate.  Variance-reducing techniques can be employed.  One can inject randomness into other types of models by replacing deterministic values by stochastic values, but in a Monte Carlo simulation, randomness is the primary focus.  Examples:

A. A Simple Stochastic Simulation on a Generic Sigmoid Activation [6]

Clear[f, a, b, c, d, x]
<< Statistics`ContinuousDistributions`

f[a_, b_, c_, d_, 
        x_] := (a/(1 + (E^(-b x + c))))+d         (* b > 0,  
          f :  R  ->  (d, 
              a + d) *) 
SeedRandom[1]

Plot[f[2, 0.5, 0, 1, x], {x, -10, 10}];

Plot[
        f[Random[NormalDistribution[2, 0.1]], 0.5, 0, 1, x], {x, -10, 
          10}];

Plot[
        f[2, 0.5, Random[NormalDistribution[0, 0.1]], 1, x], {x, -10, 
          10}];

B. Monte Carlo Volume Bounded by an Elliptic Paraboloid and xy Limits

Clear[z, x, y, vol, zbig, cuboid4th, ntry, xs, 
      ys]
z[x_, y_] := 4 - (x^2) - (1/4) (y^2)

Plot3D[z[x, y], {x, -2, 2}, {y, -4, 4}];

SeedRandom[1];

vol = 0. ;
zbig = 4;
cuboid4th = 2*4*zbig;
ntry = 100000;
Do[(xs = Random[Real, {0, 2}]; ys = Random[Real, {0, 4}]; 
          If[Random[Real, {0, zbig}] ≤ 
              z[xs, ys], vol++]), {ntry}];
Print["Vol. Est.: ", 4 ((vol/ntry)) cuboid4th]

4 (∫_0^2)(∫_0^(2 (Sqrt[ 4 - (x^2)])))(4 - 
                (x^2) - (1/4) 
                  (y^2)) \[DifferentialD]y \[DifferentialD]x
% // N

Optimization: This important category contains models that are designed to produce the best solution and typically require many controlled evaluations, to determine a minimum or maximum.

11. Continuous Optimization Models.  These models involve optimizing a function of one or more continuous independent variables, where the variables may be subject to one or more functional constraints.  Constrained optimization problems are also known as Mathematical Programming problems and come in several forms (e.g., Linear Programming, Quadratic Programming, and Nonlinear Programming).  The complexity increases with the number of variables, the amount of nonlinearity present in the functions, and the number of constraints.  Examples:

A. The Shekel Unconstrained Bivariate Maximization Problem [3]

Clear[f, x, y];
f = 
        (1/((((x - 4))^2) + (((y - 4))^2) + 0.1)) + 
          (1/((((x - 1))^2) + (((y - 1))^2) + 0.2)) + 
          (1/((((x - 8))^2) + (((y - 8))^2) + 0.2));

Plot3D[f, {x, 0, 10}, {y, 0, 10}];
ContourPlot[f, {x, 0, 10}, {y, 0, 10}];

Maximize[f, {x, y}]

FindMaximum[f, {x, 8}, {y, 8}]

FindMaximum[f, {x, 4}, {y, 4}]

B. A Simple and Small Linear Programming Problem

Clear[x, y]
Maximize[{3  x + 2  y, x + 12  y ≤ 20, x - y ≤ 2, 
        0 ≤ x ≤ 5, 
        0 ≤ y ≤ 10}, {x, y}]
% // N

C. A Small Nonlinear Programming Problem [2]

Clear[x, y]
Minimize[{Abs[x - 2] + Abs[y - 2], (y^2) ≤ x, 
        (x^2) + (y^2) == 1}, {x, y}]
% // N

12. Discrete Optimization Models.  These models are like their continuous counterparts, except that some or all of the variables are required to take on discrete values (e.g., Integer Linear Programming).  Combinatorial techniques must also be used in solving these models.  Example:

A. A Simple and Small Integer Linear Programming Problem

Clear[x, y]
NMaximize[{3  x + 2  y, x + 12  y ≤ 20, 
        x - y ≤ 2, 0 ≤ x ≤ 5, 
        0 ≤ y ≤ 10,  x  ∈  Integers,  
        y ∈  Integers}, {x, y}]

References:

1. Chicone, Carmen, Book Review, Introduction to Dynamical Systems, by M. Brin and G. Stuck, Cambridge University Press, 2002.  Appeared in SIAM Review, Vol. 46, No. 1, March 2004, pp. 157-162, SIAM Press.

2. Fiacco, Anthony, V. and Garth McCormick, Nonlinear Programming: Sequential Unconstrained Minimization Techniques, John Wiley & Sons, New York, NY, 1968, 210 pages.

3. Floudas, Christodoulos, A., Ed., Handbook of Test Problems in Local and Global Optimization, Kluwer Academic Publishers, Boston, MA, 1999, 441 pages.

4. Giordano, Frank R., et al., A First Course in Mathematical Modeling, 3rd Ed., Thompson Brooks/Cole, Pacific Grove, CA, 2003, 538 pages.

5. Noyes, James L., Artificial Intelligence with Common Lisp: Fundamentals of Symbolic and Numeric Processing, D.C. Heath & Co., Lexington, MA, 1992, 622 pages.

6. Noyes, James L., "Neural Network Training," Handbook of Neural Computation, Chapter B3, Oxford University Press, New York, NY, 1996, pp. B3:1-B3:31.

7. Noyes, James L., "A First Course in Computational Science (Why A Math Book Isn't Enough)," SIGCSE Bulletin, Chapter B3, Vol. 34, No. 1, March 2002, pp. 18-22, ACM Press.

8. Silver, Brian L., The Ascent of Science, Oxford University Press, New York, NY, 1998, 534 pages.

9. Skeel, Robert D. and Jerry B. Keiper, Elementary Numerical Computing with Mathematica, Stipes Publishing L.L.C., Champaign, IL, 2001, 434 pages.