(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "NeXT Mathematica Notebook Front End Version 2.2"; NeXTStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e8, 24, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, L1, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L1, a20, 18, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L1, a15, 14, "Times"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L1, a12, 12, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L1, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L1, 12, "Courier"; ; fontset = name, inactive, noPageBreakInGroup, nohscroll, preserveAspect, M7, italic, B65535, L1, 10, "Times"; ; fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L1, 12, "Times"; ; fontset = leftheader, 12; fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, italic, L1, 12, "Times"; ; fontset = leftfooter, 12; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Courier"; ; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; paletteColors = 128; automaticGrouping; currentKernel; ] :[font = title; inactive; noKeepOnOnePage; preserveAspect; startGroup] Polution Concentrations in Industry Lake :[font = section; inactive; preserveAspect; startGroup] BRIEF ABSTRACT :[font = subsection; inactive; preserveAspect; endGroup] A model for lake with pollutants is studied (1) first assuming instantaneous mixing and then (20 with diffusion in a two layered lake. Model results are compared for which is more realistic and why. :[font = section; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] GENERAL INFORMATION :[font = subsection; inactive; noKeepOnOnePage; preserveAspect; endGroup] FileName: POLCONC Full title: Polution Concentrations in Industry Lake Developer: Aaron Klebanoff, Department of Mathematics, Rose-Hulman Institute of Technology, Terre Haute IN 47803 USA. Contact: Brian J. Winkel, Department of Mathematics, Rose-Hulman Institute of Technology, Terre Haute IN 47803 USA. Phone: 812-877-8412. Email: winkel@rose-hulman.edu. FAX: 812-877-3198. Support: The production of this material is supported by the National Science Foundation under Division of Undergraduate Education grant DUE-9352849: Development Site for Complex, Technology-Based Problems in Calculus with Applications in Science and Engineering and the Arvin Foundation of Columbus IN. :[font = section; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] STATEMENT OF PROBLEM :[font = subsection; inactive; noKeepOnOnePage; preserveAspect] Industry Lake is bordered by a chemical plant which spews 3 tons (6000 pounds) of pollutants per minute into the lake. Numerous creeks run into the lake with a net inflow of 5000 gal/min, and a dam controls the outflow rate. The volume of the lake remains at approximately one million gallons. :[font = subsection; inactive; noKeepOnOnePage; preserveAspect] (1) If the dam allows mixed lake water to flow out at 5000 gal/min, find a function that describes the amount (in pounds) of pollutants in the lake at any time. Also, determine the long term behavior of the pollutants in the lake. Clearly state the assumptions that you are making in your model, with one of your assumptions being that the pollutants mix throughout the lake instantaneously. :[font = subsection; inactive; noKeepOnOnePage; preserveAspect] (2) Of course, the mixing is not instantaneous; in fact, studies have shown that the water in lakes is layered similar to the layering found in the atmosphere, and that small particulates actually stay close to the surface, while large bits of pollutants sink to the bottom. For a great simplification (although not as simplified as before), assume that the lake has essentially 2 layers (of equal volume) and that mixing is instantaneous only within each layer. Further assume that the transfer rate from one layer to the other is proportional to the rate is proportional to the pollution density level in the layer that the pollutants are LEAVING. (The constants of proportionality are called diffusivity constants.) You should still assume that mixing is instantaneous within each layer. Note: The assumptions outlined above are more general than assuming that pollutants diffusion from one cell to the other is proportional to the DIFFERENCE in between pollutant levels. In the case where the diffusivity constants are the same, you are making this assumption. (i) Write down the appropriate set of differential equations that models this situation. (You'll need two equations -- one for each layer.) (ii) Show that the total pollution concentration in the lake was over estimated by assuming instantaneous mixing throughout the lake. Why would this be an important result if you worked for the company? (iii) Assume that the diffusivity constants are the same. - What happens as the diffusivity constant goes to zero? - What happens as the diffusivity constant goes to infinity? Give a physical explanation for the results. :[font = subsection; inactive; noKeepOnOnePage; preserveAspect; endGroup] (3) Describe a general model that can account for noninstantaneous mixing throughout the lake. You are not expected to solve the equations, but you should be able to formulate them and carefully explain what all of the variables and parameters are. Hint: You might break the lake up into several compartments - a number in the upper level and a number in the lower level - all of which have diffusion with their neighboring compartment. :[font = section; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] KEYWORDS :[font = subsection; inactive; noKeepOnOnePage; preserveAspect; endGroup] Differential Equations; systems of linear, first order equations; Tank Problems; Pollutants; Diffusivity; Systems; Modeling. :[font = section; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] TEACHER NOTES :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] ISSUES RELATED TO THE PROBLEM :[font = subsubsection; inactive; noKeepOnOnePage; preserveAspect] This exercise encourages students to grapple with the meanings of individual terms in the differential equations. :[font = subsubsection; inactive; preserveAspect; endGroup] Questions (2) and (3) guide the students to more general models in which mixing is not assumed to be instantaneous (at least not over a large region). This helps students understand why the original seemingly unreasonable assumption was made --- namely, without the instantaneous mixing assumption, the problem is quite complicated. :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] Prerequisites :[font = subsubsection; inactive; noKeepOnOnePage; preserveAspect] For (1) and (2): Linear 1st order differential equations OR the derivative and a graphical DE Solver tool such as DSolve in Mathematica. ;[s] 3:0,0;126,1;137,2;138,-1; 3:1,10,8,Times,1,12,0,0,0;1,10,8,Times,3,12,0,0,0;1,10,8,Times,1,12,0,0,0; :[font = subsubsection; inactive; preserveAspect; endGroup] For the entire problem, the students should have been introduced to systems of ordinary differential equations. :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] Time allotment - time management :[font = subsubsection; inactive; preserveAspect; endGroup] This project should take 1 day in class to finish most of (1) and (2); and then a week out of class for (3) and (4). :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] Expectations :[font = subsubsection; inactive; preserveAspect] Students should be able to handle the Problem 1, and should also be able to figure out appropriate models for Problem 2 and 3 assuming that they have worked with other models of 1st order differential equations. :[font = subsubsection; inactive; preserveAspect; endGroup] Some students may arrive at non-linear models. You may wish to suggest out right that all models should be kept linear. In this case, there really aren't many reasonable choices for models. :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] Future payoffs :[font = subsubsection; inactive; preserveAspect] Reinforces the modeling concept for differential equations: The changes of a state variable is equal to the incoming rate minus the outgoing rate. :[font = subsubsection; inactive; preserveAspect; endGroup] Reinforces the notion that the equations are models -- not reality, and opens the door to experimentation with better models. :[font = subsection; inactive; noKeepOnOnePage; preserveAspect] Extensions :[font = subsection; inactive; noKeepOnOnePage; preserveAspect; endGroup] References and Sources :[font = section; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] POSSIBLE SOLUTION(S) :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] Problem 1. In and Out the Same --- Instantaneous Mixing :[font = subsubsection; inactive; preserveAspect] Key Assumptions: Pollutants dilute and mix uniformly throughout lake instantaneously! No evaporation (and thus, the lake volume is constant, V) Then, if P(t) = pounds of pollutants at time t (in minutes), V = volume of the lake (in gallons), Po = initial amount of pounds of pollutants at time t (in minutes); then the differential equation describing the amount of polution is: dP/dt = (incoming rate) - (outgoing rate) = 6000 - (P(t) [lb])/(V [gal])(5000 [gal/min] ) = 6000 - (5000/V) P(t) P(0) = Po. Note that since the right hand side of the DE is time dependent (i.e., the DE is autonomous), P_eq = 6000 1000000/5000 = 1.2 million pounds is an equilibrium solution. Furthermore, it's attracting, i.e., as t -> infinity, P(t) -> P_eq. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Here's how to solve the IVP (assuming the lake starts out with fresh water...) :[font = input; preserveAspect] Clear[p, t]; :[font = input; preserveAspect] V = 1000000; eqn = p'[t] == 6000 - 5000/V p[t]; ic = p[0] == 0; ans = Flatten[DSolve[{eqn, ic}, p[t], t]]; :[font = input; preserveAspect] sol = p[t] /. ans; :[font = input; preserveAspect; endGroup; endGroup] psol = Plot[sol, {t, 0, 100}] :[font = subsection; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] Problem 2. Two layers. :[font = subsubsection; inactive; preserveAspect] Assume that p1[t] = pounds of pollutant in the top layer; p2[t] = pounds of pollutant in the bottom layer; d12 = rate (in gallons/minute) that pollutants transfer from top to bottom; d21 = rate (in gallons/minute) that pollutants transfer from bottom to top; V = volume of lake (in gallons); :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Then, the differential equations that govern the pollutants in each layer are: :[font = input; preserveAspect; endGroup] V = 1000000; e1 = p1'[t] == 6000 - 5000/(V/2) p1[t] - d12/(V/2) p1[t] + d21/(V/2) p2[t]; e2 = p2'[t] == d12/(V/2) p1[t] - d21/(V/2) p2[t]; :[font = subsubsection; inactive; preserveAspect] The Fixed Point Analysis :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] shows that the unique fixed solution :[font = input; Cclosed; preserveAspect; startGroup] eq1 = 6000 - 5000/(V/2) p1[t] - d12/(V/2) p1[t] + d21/(V/2) p2[t]==0; eq2 = d12/(V/2) p1[t] - d21/(V/2) p2[t]==0; Flatten[Solve[{eq1, eq2}, {p1[t], p2[t]}]] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {p2[t] -> (600000*d12)/d21, p1[t] -> 600000} ;[o] 600000 d12 {p2[t] -> ----------, p1[t] -> 600000} d21 :[font = subsubsection; inactive; preserveAspect] is attracting. Note that the limiting weight of pollutants in the top layer (p1[t]) is independent of the diffusivities for this fixed solution. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] For those familiar with matrices, one can show that the solutions are attracting by use of the Jacobian matrix which is the matrix associated with the linear system: :[font = input; Cclosed; preserveAspect; startGroup] jac = {{- 5000/(V/2) - d12/(V/2), d21/(V/2)}, {d12/(V/2), -d21/(V/2)}} :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {{-1/100 - d12/500000, d21/500000}, {d12/500000, -d21/500000}} ;[o] 1 d12 d21 d12 -d21 {{-(---) - ------, ------}, {------, ------}} 100 500000 500000 500000 500000 :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] and the eigenvalues of this Jacobian matrix are :[font = input; Cclosed; preserveAspect; startGroup] Simplify[Eigenvalues[jac]] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {(-500000 - 100*d12 - 100*d21 - (-200000000*d21 + 10000*(5000 + d12 + d21)^2)^(1/2))/100000000\ , (-500000 - 100*d12 - 100*d21 + (-200000000*d21 + 10000*(5000 + d12 + d21)^2)^(1/2))/100000000} ;[o] {(-500000 - 100 d12 - 100 d21 - 2 Sqrt[-200000000 d21 + 10000 (5000 + d12 + d21) ]) / 100000000, (-500000 - 100 d12 - 100 d21 + 2 Sqrt[-200000000 d21 + 10000 (5000 + d12 + d21) ]) / 100000000} :[font = subsubsection; inactive; preserveAspect] which always have negative real part. So, the equilibrium solutions are stable (attracting.) :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Initial Conditions and Corresponding Solutions :[font = text; inactive; preserveAspect] If we assume that the lake is initially clean, we have the following initial conditions: :[font = input; preserveAspect] ic1 = p1[0] == 0; ic2 = p2[0] == 0; :[font = text; inactive; preserveAspect] Here's the general solutions: sol1 and sol2 (in terms of the unspecified parameters d12, d21, and V.) :[font = input; preserveAspect] ans = Flatten[DSolve[{e1, e2, ic1, ic2}, {p1[t], p2[t]}, t]]; :[font = input; preserveAspect; endGroup] sol1 = p1[t] /. ans[[1]]; sol2 = p2[t] /. ans[[2]]; :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Here are some sample graphs of solutions and comments on the results... :[font = text; inactive; preserveAspect] Here we show that when d12 = d21 is small, p1(t) and p2(t) are quite different and p1(t) > p2(t) (which makes sense since most particulate matter remains near the surface where it came into the lake water), while as d12 = d21 -> infinity, mixing becomes instantaneous throughout the lake and p1(t) and p2(t) coalesce. The animation further demonstrates that p2(t) increases from nearly zero (the initial weight of pollutants) as the diffusivity increases, while p1(t) decreases with increasing diffusivity since more of the pollutants are able to cross to the other layer. (In the animation, the parameter d12 = d21 increases exponentially. The thicker curve is p1(t).) :[font = input; preserveAspect; endGroup] Do[Plot[{sol1 /. {d12 -> N[Exp[i]], d21 -> N[Exp[i]]}, sol2 /. {d12 -> N[Exp[i]], d21 -> N[Exp[i]]}}, {t, 0, 100}, PlotStyle -> {AbsoluteThickness[2], AbsoluteThickness[1]}, AxesLabel -> {"minutes", "pollutant weight"}, PlotRange -> {{0, 100}, {0, 350000}}], {i, 5, 13}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Here, we consider the affects of changing the ratio of the diffusivities. This is done by letting alpha = d12/d21 and then fixing -- say d21 = 1 while varying d12 both above and below 1000. In the animation, the parameter alpha increases exponentially. The thicker curve is p1(t). Note that when alpha is small, most pollutants remain in the top layer, but when alpha is very large, most pollutants end up in the bottom layer. :[font = input; preserveAspect; endGroup; animationSpeed = 48] Do[Plot[{sol1 /. {d21 -> 1., d12 -> N[Exp[alpha]]}, sol2 /. {d21 -> 1., d12 -> N[Exp[alpha]]}}, {t, 0, 100}, PlotStyle -> {AbsoluteThickness[2], AbsoluteThickness[1]}, AxesLabel -> {"minutes", "pollutant weight"}, PlotRange -> {{0, 100}, {0, 350000}}], {alpha, 6, 12, 0.5}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Below we compare the graphs of the solutions when the lake is assumed to be in two layers vs. when the lake is assumed to mix entirely instantaneously. The bold graph corresponds to no layers (instantaneous mixing throughout), while the thinner graph corresponds to instantaneous mixing within each layer. This animation shows that as the diffusivities increase and thus mixing occurs more readily, that the simplified model becomes a reasonable approximation. The parameters d12 = d21 increase exponentially. This is an important result for the company dumping pollutants since the simplified model over estimates the amount of damage done to the environment. The solution of the simplified model (from #1) is the thicker curve and does not change. :[font = input; preserveAspect; endGroup; endGroup] Do[Plot[{sol1+sol2 /. {d21 -> N[Exp[i]], d12 -> N[Exp[i]]}, sol}, {t, 0, 100}, PlotStyle -> {AbsoluteThickness[1], AbsoluteThickness[2]}, AxesLabel -> {"minutes", "pollutant weight"}, PlotRange -> {{0, 100}, {0, 500000}}], {i, 5, 12, 0.5}] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Problem 3. Multiple Layers and Cells. :[font = subsubsection; inactive; preserveAspect] There are many levels of sophistication that one can impose here. To start simply, one might partition the lake by making lots of horizontal slabs (say, of equal volume) and assuming that there is instantaneous mixing occurs within each horizontal slab. If p_i(t) = pounds of pollutants in slab i, then the model would be p_1'(t) = 6000 - (n/V) { 5000 p_1(t) + d21 p_2(t) - d12 p_1(t) } p_2'(t) = (n/V) { d12 p_1(t) + d32 p_3(t) - (d21 + d23) p_2(t) } p_3'(t) = (n/V) { d23 p_2(t) + d43 p_4(t) - (d32 + d34) p_3(t) } p_4'(t) = (n/V) { d34 p_3(t) + d54 p_5(t) - (d43 + d45) p_4(t) } . . . p_n'(t) = ... It's a bit easier to see from the Jacobian matrix that this system is tridiagonal. At a higher level of complexity, one can extend to two or even three dimensions. In these cases, the modeler must decide how pollutants can diffuse in higher dimensions. A reasonable assumption is that there is only diffusion across borders of dimension one less than the dimension of the cell. For example, in two dimensions (the lake is partitioned into thin cylinders with small rectangular cross sections), pollutants move across the lines (1D) between neighboring cells, but not vertices at the diagonals. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Code for Grid :[font = input; preserveAspect; endGroup; endGroup; endGroup] Show[Graphics[{Line[{{0, 0}, {4, 0}}], Line[{{0, 0}, {0, 4}}], Line[{{0, 4}, {4, 4}}], Line[{{4, 4}, {4, 0}}], Line[{{1, 0}, {1, 4}}], Line[{{2, 0}, {2, 4}}], Line[{{3, 0}, {3, 4}}], Line[{{0, 1}, {4, 1}}], Line[{{0, 2}, {4, 2}}], Line[{{0, 3}, {4, 3}}], Text["Diffusion This Way", {1.3, 1.4}], Line[{{1.8, 1.5}, {2.2, 1.5}}], Line[{{2.1, 1.4}, {2.2, 1.5}}], Line[{{2.1, 1.6}, {2.2, 1.5}}], Text["But Not This Way", {2.3, 2.4}], Line[{{2.2, 2.2}, {1.8, 1.8}}], Line[{{1.8, 1.8}, {1.85, 1.9}}], Line[{{1.8, 1.8}, {1.9, 1.85}}]}, AspectRatio -> Automatic]] :[font = section; inactive; Cclosed; noKeepOnOnePage; preserveAspect; startGroup] ISSUES IN SOLUTION :[font = subsection; inactive; preserveAspect] In (3), students shouldn't have any trouble realizing that each horizontal sheet of water can in turn be partitioned. Students can then grapple with the problems of how diffusion should work in two or even three dimensions. Some will find it helpful if you suggest considering just one vertical division to have four compartments at first. :[font = subsection; inactive; noKeepOnOnePage; preserveAspect; endGroup; endGroup] It is far more important that the students are allowed time to grapple with modeling issues instead of the method of solution. While the methods involved with solving intermediate parts of the problem are interesting in their own right, if the students are forced to dwell on the details in this problem, they will never get to the major issues of the problem. Thus, use of a CAS such as Mathematica or Maple to solve the DEs quickly is advisable. ;[s] 3:0,0;390,1;401,2;449,-1; 3:1,12,9,Times,1,14,0,0,0;1,13,10,Times,3,14,0,0,0;1,12,9,Times,1,14,0,0,0; ^*)