(* Content-type: application/mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 6.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       145,          7]
NotebookDataLength[     16988,        443]
NotebookOptionsPosition[     15804,        401]
NotebookOutlinePosition[     16238,        420]
CellTagsIndexPosition[     16195,        417]
WindowFrame->Normal
ContainsDynamic->False*)

(* Beginning of Notebook Content *)
Notebook[{
Cell[BoxData["\[IndentingNewLine]"], "Input",
 CellChangeTimes->{3.441646421872959*^9}],

Cell[CellGroupData[{

Cell["\<\
Explicit and implicit methods for the catalytic converter model (eqs 5.19-23)\
\
\>", "Subtitle",
 CellChangeTimes->{{3.441466087791348*^9, 3.441466094325899*^9}, {
  3.445250472391621*^9, 3.445250504189272*^9}, {3.445266580359993*^9, 
  3.445266580670422*^9}}],

Cell[TextData[StyleBox["We' ll discuss the algorithms this code is based on \
on Tuesday. Let me know if you' re working with this before then and have \
questions or find bugs.",
 FontSize->14]], "Text",
 CellChangeTimes->{{3.445323907870487*^9, 3.445323961765668*^9}}],

Cell[CellGroupData[{

Cell["Definitions", "Subsection",
 CellChangeTimes->{{3.441970934922043*^9, 3.441970935706514*^9}}],

Cell[TextData[StyleBox["Beware of bugs! I haven' t vetted this very \
carefully.",
 FontSize->14]], "Text",
 CellChangeTimes->{{3.4452505181235447`*^9, 3.445250531182003*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"dd", "[", 
   RowBox[{"f_", ",", "s_"}], "]"}], ":=", 
  RowBox[{
   RowBox[{"Append", "[", 
    RowBox[{
     RowBox[{"Drop", "[", 
      RowBox[{"f", ",", "1"}], "]"}], ",", 
     RowBox[{"Last", "[", "f", "]"}]}], "]"}], "-", 
   RowBox[{"2", "f"}], "+", 
   RowBox[{"Prepend", "[", 
    RowBox[{
     RowBox[{"Drop", "[", 
      RowBox[{"f", ",", 
       RowBox[{"-", "1"}]}], "]"}], ",", "s"}], "]"}]}]}]], "Input",
 CellChangeTimes->{{3.44146335106499*^9, 3.441463513484474*^9}, {
  3.4451706190981092`*^9, 3.445170624518355*^9}, {3.445170706489233*^9, 
  3.445170730772644*^9}, {3.4452170851972437`*^9, 3.4452170873357077`*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"invddmat", "[", 
   RowBox[{"n_", ",", "h_"}], "]"}], ":=", 
  RowBox[{"Inverse", "[", 
   RowBox[{
    RowBox[{"IdentityMatrix", "[", "n", "]"}], "-", 
    RowBox[{"h", " ", 
     RowBox[{"Table", "[", 
      RowBox[{
       RowBox[{"Switch", "[", 
        RowBox[{
         RowBox[{"Abs", "[", 
          RowBox[{"i", "-", "j"}], "]"}], ",", "1", ",", "1", ",", "0", ",", 
         RowBox[{"-", "2"}], ",", "_", ",", "0"}], "]"}], ",", 
       RowBox[{"{", 
        RowBox[{"i", ",", "n"}], "}"}], ",", 
       RowBox[{"{", 
        RowBox[{"j", ",", "n"}], "}"}]}], "]"}]}]}], "]"}]}]], "Input",
 CellChangeTimes->{{3.445214338412882*^9, 3.445214357593511*^9}, {
  3.445214390822422*^9, 3.4452144847350283`*^9}, {3.44524981295238*^9, 
  3.445249823846754*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"scaleu", "[", 
   RowBox[{"t_", ",", "\[CapitalDelta]x_"}], "]"}], ":=", 
  RowBox[{"Drop", "[", 
   RowBox[{
    RowBox[{"FoldList", "[", 
     RowBox[{
      RowBox[{
       RowBox[{
        RowBox[{"(", 
         RowBox[{"1", "-", 
          RowBox[{"\[CapitalDelta]x", 
           FractionBox["#2", 
            RowBox[{"1", "+", "#2"}]]}]}], ")"}], "#1"}], "&"}], ",", "1", 
      ",", "t"}], "]"}], ",", "1"}], "]"}]}]], "Input",
 CellChangeTimes->{{3.445169334287265*^9, 3.445169422541273*^9}, {
   3.4451694873614063`*^9, 3.445169557617146*^9}, {3.445169606039524*^9, 
   3.44516960977407*^9}, {3.445169795047518*^9, 3.445169821264207*^9}, 
   3.445169863999428*^9}],

Cell[TextData[StyleBox["cattemp' s inputs are a function s of t (time), a \
\"bump function\" b of x (space) that should be concentrated near 0 (b \
smoothes the heat point source used in the text, to avoid horrific stability \
issues for coarse meshes), the raw exhaust constant u0, the length of the \
converter and number of spatial mesh points to use, the time interval and \
number of spatial mesh points, and the method specification \"exp\" \
(explicit) or \"imp\" (implicit). cattemp returns the approximate temperatue \
over the temporal/spatial mesh.",
 FontSize->14]], "Text",
 CellChangeTimes->{{3.445322917155437*^9, 3.445323209098837*^9}, {
  3.445323305782783*^9, 3.4453233156395187`*^9}, {3.4453238712554903`*^9, 
  3.44532388479429*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"cattemp", "[", 
   RowBox[{"s_", ",", "b_", ",", "u0_", ",", 
    RowBox[{"{", 
     RowBox[{"length_", ",", "xn_"}], "}"}], ",", 
    RowBox[{"{", 
     RowBox[{"tmax_", ",", "tn_"}], "}"}], ",", "method_"}], "]"}], ":=", 
  RowBox[{"Block", "[", 
   RowBox[{
    RowBox[{"{", 
     RowBox[{
      RowBox[{"\[CapitalDelta]x", "=", 
       RowBox[{"N", "[", 
        FractionBox["length", "xn"], "]"}]}], ",", 
      RowBox[{"\[CapitalDelta]t", "=", 
       RowBox[{"N", "[", 
        FractionBox["tmax", "tn"], "]"}]}], ",", "mat", ",", "update"}], 
     "}"}], ",", 
    RowBox[{
     RowBox[{"If", "[", 
      RowBox[{
       RowBox[{"method", "==", "\"\<imp\>\""}], ",", 
       RowBox[{
        RowBox[{"mat", "=", 
         RowBox[{"invddmat", "[", 
          RowBox[{
           RowBox[{"xn", "+", "1"}], ",", 
           FractionBox["\[CapitalDelta]t", 
            SuperscriptBox["\[CapitalDelta]x", "2"]]}], "]"}]}], ";", 
        RowBox[{"update", "=", 
         RowBox[{
          RowBox[{"Drop", "[", 
           RowBox[{
            RowBox[{"mat", ".", 
             RowBox[{"Prepend", "[", 
              RowBox[{
               RowBox[{"#1", "+", " ", 
                RowBox[{"\[CapitalDelta]t", " ", "u0", " ", 
                 RowBox[{"scaleu", "[", 
                  RowBox[{"#1", ",", "\[CapitalDelta]x"}], "]"}]}]}], ",", 
               "#2"}], "]"}]}], ",", "1"}], "]"}], "&"}]}]}], ",", 
       RowBox[{"update", "=", 
        RowBox[{
         RowBox[{"#1", "+", 
          RowBox[{"\[CapitalDelta]t", 
           RowBox[{"(", " ", 
            RowBox[{
             FractionBox[
              RowBox[{"dd", "[", 
               RowBox[{"#1", ",", "#2"}], "]"}], 
              SuperscriptBox["\[CapitalDelta]x", "2"]], "+", 
             RowBox[{"u0", " ", 
              RowBox[{"scaleu", "[", 
               RowBox[{"#1", ",", "\[CapitalDelta]x"}], "]"}]}]}], ")"}]}]}], 
         "&"}]}]}], "]"}], ";", "\[IndentingNewLine]", 
     RowBox[{"FoldList", "[", 
      RowBox[{"update", ",", 
       RowBox[{
        RowBox[{"(", 
         RowBox[{"s", "/.", 
          RowBox[{"t", "\[Rule]", "0"}]}], ")"}], 
        RowBox[{"Array", "[", 
         RowBox[{
          RowBox[{
           RowBox[{"(", 
            RowBox[{"b", "/.", 
             RowBox[{"x", "\[Rule]", 
              RowBox[{"\[CapitalDelta]x", "#"}]}]}], " ", ")"}], "&"}], ",", 
          "xn"}], "]"}]}], ",", 
       RowBox[{"Array", "[", 
        RowBox[{
         RowBox[{
          RowBox[{"(", 
           RowBox[{"s", "/.", 
            RowBox[{"t", "\[Rule]", 
             RowBox[{"\[CapitalDelta]t", "#"}]}]}], " ", ")"}], "&"}], ",", 
         "tn"}], "]"}]}], "]"}]}]}], "]"}]}]], "Input",
 CellChangeTimes->{{3.441463670749092*^9, 3.441464003865231*^9}, {
   3.441464181953949*^9, 3.441464207185364*^9}, {3.441465071991949*^9, 
   3.441465073837398*^9}, {3.441465697379727*^9, 3.441465769528563*^9}, {
   3.441465819936953*^9, 3.4414658343690557`*^9}, {3.4414662161819153`*^9, 
   3.441466272030188*^9}, {3.4419294681490107`*^9, 3.44192946928447*^9}, {
   3.441929509966228*^9, 3.441929510093254*^9}, {3.441929753781815*^9, 
   3.441929792585766*^9}, {3.441929823028234*^9, 3.441929830511952*^9}, {
   3.4419300472683697`*^9, 3.441930052716896*^9}, {3.4421949894461813`*^9, 
   3.442195034721233*^9}, {3.442195100999946*^9, 3.442195111046646*^9}, {
   3.442195145173451*^9, 3.442195156196533*^9}, {3.4421953894861813`*^9, 
   3.442195426742148*^9}, {3.442195537991158*^9, 3.442195546496901*^9}, {
   3.4421956325578136`*^9, 3.442195634565958*^9}, {3.4421957498965397`*^9, 
   3.4421958391874237`*^9}, {3.4421959113246*^9, 3.4421959135115833`*^9}, {
   3.442196005649219*^9, 3.4421960074846582`*^9}, {3.44516795834851*^9, 
   3.445167977054208*^9}, {3.445168330431795*^9, 3.445168335603718*^9}, {
   3.445168577077656*^9, 3.445168621169547*^9}, {3.4451688996452847`*^9, 
   3.4451689053175287`*^9}, {3.4451690318086033`*^9, 3.445169033699854*^9}, {
   3.4451699922998753`*^9, 3.4451699938958883`*^9}, {3.4451700461126213`*^9, 
   3.445170114716936*^9}, {3.445170199272728*^9, 3.445170204906458*^9}, {
   3.445170823260037*^9, 3.4451709462222033`*^9}, {3.445171014610458*^9, 
   3.4451711347282963`*^9}, {3.4451712663348703`*^9, 3.445171279469141*^9}, {
   3.445171335477211*^9, 3.445171373262429*^9}, {3.445172089672163*^9, 
   3.44517215166234*^9}, {3.445214568123427*^9, 3.445214570680687*^9}, {
   3.445214605077454*^9, 3.445214735352178*^9}, {3.445214834520021*^9, 
   3.44521483495984*^9}, {3.44521489497087*^9, 3.445214904610998*^9}, {
   3.4452150415001593`*^9, 3.445215056458786*^9}, {3.445215114676882*^9, 
   3.445215120942988*^9}, {3.445215246456839*^9, 3.4452152650997066`*^9}, {
   3.445215440280065*^9, 3.4452154410752487`*^9}, {3.445215783953497*^9, 
   3.445215922507839*^9}, {3.445215954271927*^9, 3.445216051171646*^9}, {
   3.445216081510882*^9, 3.44521614837773*^9}, 3.445216180547165*^9, {
   3.445217090688387*^9, 3.445217091322757*^9}}],

Cell[TextData[StyleBox["emred' s inputs are a temperature grid (e.g. an \
output of cattemp), and the associated mesh increments.  emred returns a \
scalar: the total relative emission from the converter over the grid time \
interval (multiply by the raw exhaust constant u0 to get the actual total \
emissions; divide by the length tmax of the time interval to get the ratio of \
pollutant released with the converter to that which would be released without \
it) and a plot of the instantaneous emissions over the time interval.",
 FontSize->14]], "Text",
 CellChangeTimes->{{3.44532323560203*^9, 3.4453232867842903`*^9}, {
  3.4453233280900593`*^9, 3.4453234476923227`*^9}, {3.4453236075703707`*^9, 
  3.445323673252591*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"emred", "[", 
   RowBox[{"temp_", ",", 
    RowBox[{"{", 
     RowBox[{"\[CapitalDelta]x_", ",", "\[CapitalDelta]t_"}], "}"}]}], "]"}], 
  ":=", 
  RowBox[{
   RowBox[{
    RowBox[{"(", 
     RowBox[{
      RowBox[{"Print", "[", 
       RowBox[{"\[CapitalDelta]t", " ", 
        RowBox[{"Plus", "@@", "#"}]}], "]"}], ";", 
      RowBox[{"ListPlot", "[", 
       RowBox[{"#", ",", 
        RowBox[{"Joined", "\[Rule]", "True"}], ",", 
        RowBox[{"PlotRange", "\[Rule]", "All"}]}], "]"}]}], ")"}], "&"}], "[",
    
   RowBox[{
    RowBox[{
     RowBox[{"Last", "[", 
      RowBox[{"scaleu", "[", 
       RowBox[{"#", ",", "\[CapitalDelta]x"}], "]"}], "]"}], "&"}], "/@", 
    "temp"}], "]"}]}]], "Input",
 CellChangeTimes->{{3.445216609544056*^9, 3.445216697474658*^9}, {
  3.445216820055265*^9, 3.445216918595252*^9}, {3.44521695560662*^9, 
  3.445216967433098*^9}, {3.445217006292432*^9, 3.445217018251143*^9}, {
  3.4452173577801313`*^9, 3.445217360019902*^9}}],

Cell[TextData[StyleBox["plotdata generates a 3 - D plot of a data grid.",
 FontSize->14]], "Text",
 CellChangeTimes->{{3.445323466899743*^9, 3.4453235369778013`*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"plotdata", "[", 
   RowBox[{"data_", ",", 
    RowBox[{"{", 
     RowBox[{"length_", ",", "xn_"}], "}"}], ",", 
    RowBox[{"{", 
     RowBox[{"tmax_", ",", "tn_"}], "}"}]}], "]"}], ":=", 
  RowBox[{"ListPlot3D", "[", 
   RowBox[{
    RowBox[{"Join", "@@", 
     RowBox[{"Table", "[", 
      RowBox[{
       RowBox[{"{", 
        RowBox[{
         RowBox[{"i", " ", 
          FractionBox["length", "xn"]}], ",", 
         RowBox[{"j", " ", 
          FractionBox["tmax", "tn"]}], ",", 
         RowBox[{"data", "\[LeftDoubleBracket]", 
          RowBox[{"j", ",", "i"}], "\[RightDoubleBracket]"}]}], "}"}], ",", 
       RowBox[{"{", 
        RowBox[{"i", ",", "xn"}], "}"}], ",", 
       RowBox[{"{", 
        RowBox[{"j", ",", "tn"}], "}"}]}], "]"}]}], ",", 
    RowBox[{"PlotRange", "\[Rule]", "All"}], ",", 
    RowBox[{"Mesh", "\[Rule]", "False"}], ",", 
    RowBox[{"Boxed", "\[Rule]", "False"}]}], "]"}]}]], "Input",
 CellChangeTimes->{{3.44146493522213*^9, 3.4414651792440977`*^9}, 
   3.4414652539022427`*^9, {3.441465374856699*^9, 3.441465375353839*^9}, {
   3.44146542424645*^9, 3.44146542551751*^9}, {3.4414655064732637`*^9, 
   3.441465533372279*^9}, {3.441465855435968*^9, 3.4414659898892193`*^9}, {
   3.4419387778695307`*^9, 3.4419387804526157`*^9}, {3.445173451107696*^9, 
   3.445173495855803*^9}, {3.4452504169009647`*^9, 3.4452504250024757`*^9}}]
}, Open  ]],

Cell[CellGroupData[{

Cell["Sample run", "Subsection",
 CellChangeTimes->{{3.441970934922043*^9, 3.441970935706514*^9}, {
  3.44525025871791*^9, 3.445250260395986*^9}}],

Cell[TextData[StyleBox["N.B. I didn' t experiment to find good \
\[CapitalDelta]x and \[CapitalDelta]t; please play with that to find \
reasonable values. Also, if you have time, please experiment with other \
choices of s (t) beyond the two very simple choices given in the homework \
problem.",
 FontSize->14]], "Text",
 CellChangeTimes->{{3.445323685811496*^9, 3.445323814123921*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"sample", "=", 
   RowBox[{"cattemp", "[", 
    RowBox[{
     RowBox[{"If", "[", 
      RowBox[{
       RowBox[{"t", "<", ".5"}], ",", ".5", ",", "1.5"}], "]"}], ",", 
     SuperscriptBox["E", 
      RowBox[{
       RowBox[{"-", "10"}], 
       SuperscriptBox["x", "2"]}]], ",", "1", ",", 
     RowBox[{"{", 
      RowBox[{"10", ",", "100"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"1", ",", "200"}], "}"}], ",", "\"\<imp\>\""}], "]"}]}], 
  ";"}]], "Input",
 CellChangeTimes->{{3.445250022008005*^9, 3.445250030761363*^9}}],

Cell[BoxData[
 RowBox[{"Animate", "[", 
  RowBox[{
   RowBox[{"ListPlot", "[", 
    RowBox[{
     RowBox[{"sample", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}],
      ",", 
     RowBox[{"Joined", "\[Rule]", "True"}], ",", 
     RowBox[{"PlotRange", "\[Rule]", 
      RowBox[{"{", 
       RowBox[{
        RowBox[{"{", 
         RowBox[{"0", ",", "100"}], "}"}], ",", 
        RowBox[{"{", 
         RowBox[{"0", ",", "1"}], "}"}]}], "}"}]}]}], "]"}], ",", 
   RowBox[{"{", 
    RowBox[{"k", ",", "1", ",", "200", ",", "20"}], "}"}]}], "]"}]], "Input",
 CellChangeTimes->{{3.445250037305208*^9, 3.445250096332851*^9}}],

Cell[BoxData[
 RowBox[{"plotdata", "[", 
  RowBox[{"sample", ",", 
   RowBox[{"{", 
    RowBox[{"10", ",", "100"}], "}"}], ",", 
   RowBox[{"{", 
    RowBox[{"1", ",", "200"}], "}"}]}], "]"}]], "Input",
 CellChangeTimes->{{3.445250276738266*^9, 3.445250288545065*^9}}],

Cell[BoxData[
 RowBox[{"emred", "[", 
  RowBox[{"sample", ",", 
   RowBox[{"{", 
    RowBox[{
     FractionBox["1", "10"], ",", 
     FractionBox["1", "200"]}], "}"}]}], "]"}]], "Input",
 CellChangeTimes->{{3.445171203450555*^9, 3.4451712144680147`*^9}, {
   3.445171324391914*^9, 3.4451713408625937`*^9}, {3.445171383791641*^9, 
   3.445171397907793*^9}, {3.445171433839017*^9, 3.44517143844263*^9}, {
   3.445171635632938*^9, 3.445171658939595*^9}, {3.445171736659377*^9, 
   3.445171737963193*^9}, 3.4451719907644978`*^9, {3.445172174137066*^9, 
   3.445172339489159*^9}, 3.445173232714061*^9, {3.445214047088072*^9, 
   3.445214055686425*^9}, {3.44521483009029*^9, 3.4452148318237534`*^9}, {
   3.445215284526841*^9, 3.445215290682805*^9}, {3.445216153008727*^9, 
   3.445216165388529*^9}, {3.445217238621779*^9, 3.445217242826477*^9}, {
   3.445217469398917*^9, 3.445217512448413*^9}, {3.4452501587892723`*^9, 
   3.445250159453452*^9}}]
}, Open  ]]
}, Open  ]]
},
WindowSize->{992, 670},
WindowMargins->{{Automatic, -7}, {Automatic, 0}},
PrintingCopies->1,
PrintingPageRange->{1, Automatic},
ShowSelection->True,
Magnification->1.,
FrontEndVersion->"6.0 for Mac OS X PowerPC (32-bit) (April 20, 2007)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[568, 21, 87, 1, 43, "Input"],
Cell[CellGroupData[{
Cell[680, 26, 271, 6, 46, "Subtitle"],
Cell[954, 34, 270, 4, 28, "Text"],
Cell[CellGroupData[{
Cell[1249, 42, 99, 1, 34, "Subsection"],
Cell[1351, 45, 175, 3, 28, "Text"],
Cell[1529, 50, 670, 18, 27, "Input"],
Cell[2202, 70, 798, 21, 27, "Input"],
Cell[3003, 93, 709, 19, 46, "Input"],
Cell[3715, 114, 754, 11, 79, "Text"],
Cell[4472, 127, 5007, 108, 162, "Input"],
Cell[9482, 237, 728, 10, 79, "Text"],
Cell[10213, 249, 1002, 28, 27, "Input"],
Cell[11218, 279, 166, 2, 28, "Text"],
Cell[11387, 283, 1401, 33, 63, "Input"]
}, Open  ]],
Cell[CellGroupData[{
Cell[12825, 321, 146, 2, 34, "Subsection"],
Cell[12974, 325, 387, 6, 45, "Text"],
Cell[13364, 333, 565, 17, 40, "Input"],
Cell[13932, 352, 628, 17, 27, "Input"],
Cell[14563, 371, 268, 7, 27, "Input"],
Cell[14834, 380, 942, 17, 47, "Input"]
}, Open  ]]
}, Open  ]]
}
]
*)

(* End of internal cache information *)