(************** Content-type: application/mathematica ************** CreatedBy='Mathematica 5.0' Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 53379, 1421]*) (*NotebookOutlinePosition[ 83593, 2524]*) (* CellTagsIndexPosition[ 83549, 2520]*) (*WindowFrame->Normal*) Notebook[{ Cell["The Seesaw/Pendulum Process", "Title", CellMargins->{{0, -87}, {Inherited, Inherited}}, Evaluatable->False, CellLabelMargins->{{0, Inherited}, {Inherited, Inherited}}, TextAlignment->Left, AspectRatioFixed->True], Cell[TextData[StyleBox["Modeling, Controller Design, Code Generation,\nand \ Simulation", FontSize->25]], "Subtitle", TextAlignment->Left, TextJustification->0], Cell["Copyright \[Copyright] MathCore AB 1999-2004", "Copyright"], Cell[BoxData[ \(\(\(CSP\)\(=\)\(False\)\(\ \)\( (*\ Whether\ Control\ System\ Professional\ is\ installed\ *) \)\(\ \ \)\)\)], "Input"], Cell[CellGroupData[{ Cell["Introduction", "Section"], Cell["\<\ This notebook illustrates how to build a dynamic model of a \ physical system and then generate efficient code for simulation.\ \>", "Text"], Cell["\<\ The physical system is a seesaw where two carts C1 and C2 can be \ moved along two tracks by controlling the forces acting on them. Cart C1 \ carries a weight and cart C2 has an inverted pendulum attached by a friction \ free joint.\ \>", "Text"], Cell[TextData[{ "The model of the open loop system has the forces ", Cell[BoxData[ \(TraditionalForm\`F\_1\)]], " and ", Cell[BoxData[ \(TraditionalForm\`F\_2\)]], " acting on each cart as inputs. We first model the system by Lagrangian \ methodology and then linearize the obtained nonlinear model. The linearized \ model is used for control system design where the Control System Professional \ (CSP) application package is used. The closed loop system are then simulated \ both within ", StyleBox["Mathematica", FontSlant->"Italic"], " and by external code generated using ", StyleBox["MathCode C++", FontSlant->"Italic"], "." }], "Text"], Cell[GraphicsData["Bitmap", "\<\ CF5dJ6E]HGAYHf4PAg9QL6QYHggo1P22o`H0lOl0 0>Oo1P2>o`H0jol00>7o1P2Jo`H0iOl00=_o1P2Vo`H0gol00=Go1P2bo`H0 fOl0000fo`030?oo03?o00<0ool0nOl00=;o:`0mocL0nol00=;o00<0 ool09ol00`3oo`1Fo`030?oo0?oo5Ol00=;o00<0ool09ol00`3oo`1Fo`03 0?oo0?oo5Ol00=;o00<0ool09ol00`3oo`1Go`030?oo0?oo5?l00=;o00<0 ool09ol00`3oo`1Go`030?oo0?oo5?l00=;o00<0ool09ol00`3oo`1Go`03 0?oo0?oo5?l00=;o00<0ool09ol00`3oo`1Go`030?oo0?oo5?l00=;o00<0 ool09ol00`3oo`1Go`030?oo0?oo5?l00=;o00<0ool09ol00`3oo`1Ho`03 0?oo0?oo4ol00=;o00<0ool09ol00`3oo`1Ho`030?oo0?oo4ol00=;o00<0 ool09ol00`3oo`1Ho`030?oo0?oo4ol004Oo100o`00oomMo`030?oo 0?oo3_l00?ooGOl00`3oo`3oo`ko003ooego00<0ool0ool>o`00oomMo`03 0?oo0?oo3_l00?ooG_l00`3oo`3oo`go003ooeko00<0ool0ool=o`00oomN o`030?oo0?oo3Ol00?ooG_l00`3oo`3oo`go003ooeko00<0ool0ool=o`00 oomOo`030?oo0?oo3?l00?ooGol00`3oo`3oo`co003ooeoo00<0ool0ool< o`00oomOo`030?oo0?oo3?l00?ooGol00`3oo`3oo`co003oof3o00<0ool0 ool;o`00oomPo`030?oo0?oo2ol00?ooH?l00`3oo`3oo`_o003oof3o00<0 ool0ool;o`00oomPo`030?oo0?oo2ol00?ooHOl00`3oo`3oo`[o003oof7o 00<0ool0ool:o`00oomQo`030?oo0?oo2_l00?ooHOl00`3oo`3oo`[o003o of7o00<0ool0ool:o`00oomQo`030?oo0?oo2_l00?ooH_l00`3oo`3oo`Wo 003oof;o00<0ool0ool9o`00oomRo`030?oo0?oo2Ol00?ooH_l00`3oo`3o o`Wo003oof;o00<0ool0ool9o`00oomSo`030?oo0?oo2?l00?ooHol00`3o o`3oo`So003oof?o00<0ool0ool8o`00oomSo`030?oo0?oo2?l00?ooHol0 0`3oo`3oo`So003oofCo00<0ool0ool7o`00oomTo`030?oo0?oo1ol00?oo I?l00`3oo`3oo`Oo003oofCo00<0ool0ool7o`00oomTo`030?oo0?oo1ol0 0?ooI?l00`3oo`3oo`Oo003oofGo00<0ool0ool6o`00oomUo`030?oo0?oo 1_l00?ooIOl00`3oo`3oo`Ko003oofGo00<0ool0ool6o`00oomUo`030?oo 0?oo1_l00?ooI_l00`3oo`3oo`Go003oofKo00<0ool0ool5o`00oomVo`03 0?oo0?oo1Ol00?ooI_l00`3oo`3oo`Go003oofKo00<0ool0ool5o`00oomW o`030?oo0?oo1?l00?ooIol00`3oo`3oo`Co003oofOo00<0ool0ool4o`00 oomWo`030?oo0?oo1?l00?ooIol00`3oo`3oo`Co003oofSo00<0ool0ool3 o`00oomXo`030?oo0?oo0ol00?ooJ?l00`3oo`3oo`?o003oofSo00<0ool0 ool3o`00oomXo`030?oo0?oo0ol00?ooJ?l00`3oo`3oo`?o003oofWo00<0 ool0ool2o`00oomYo`030?oo0?oo0_l00?ooJOl00`3oo`3oo`;o003oofWo 00<0ool0ool2o`00oomYo`030?oo0?oo0_l00?ooJ_l00`3oo`3oo`7o003o of[o00<0ool0ool1o`00oomZo`030?oo0?oo0Ol00?ooJ_l00`3oo`3oo`7o 003oof[o00<0ool0ool1o`00oom[o`030?oo0?oo003oof_o00<0ool0ool0 0?ooJol00`3oo`3oo`00oom[o`030?oo0?oo003oof_o00<0ool0ool00?oo K?l00`3oo`3no`00oom/o`030?oo0?ko003oofco00<0ool0o_l00?ooK?l0 0`3oo`3no`00oom/o`030?oo0?ko003oofco00<0ool0o_l00?ooKOl00`3o o`3mo`00oom]o`030?oo0?go003oofgo00<0ool0oOl00?ooKOl00`3oo`3m o`00oom]o`030?oo0?go003oofko00<0ool0o?l00?ooK_l00`3oo`3lo`00 oom^o`030?oo0?co003oofko00<0ool0o?l00?ooK_l00`3oo`3lo`00oom_ o`030?oo0?_o003oofoo00<0ool0nol00?ooKol00`3oo`3ko`00oom_o`03 0?oo0?_o003oofoo00<0ool0nol00?ooKol00`3oo`3ko`00oom`o`030?oo 0?[o003oog3o00<0ool0n_l00?ooL?l00`3oo`3jo`00oom`o`030?oo0?[o 003oog3o00<0ool0n_l00?ooLOl00`3oo`3io`00oomao`030?oo0?Wo003o og7o00<0ool0nOl00?ooLOl00`3oo`3io`00oomao`030?oo0?Wo003oog;o 00<0ool0n?l00?ooL_l00`3oo`3ho`00oombo`030?oo0?So003oog;o00<0 ool0n?l00?ooL_l00`3oo`3ho`00oomco`030?oo0?Oo003oog?o00<0ool0 mol00?ooLol00`3oo`3go`00oomco`030?oo0?Oo003oog?o00<0ool0mol0 0?ooLol00`3oo`3go`00oomdo`030?oo0?Ko003oogCo00<0ool0m_l00?oo M?l00`3oo`3fo`00oomdo`030?oo0?Ko003oogCo00<0ool0m_l00?ooMOl0 0`3oo`3eo`00oomeo`030?oo0?Go003oogGo00<0ool0mOl00?ooMOl00`3o o`3eo`00oomeo`030?oo0?Go003oogKo00<0ool0m?l00?ooM_l00`3oo`3d o`00oomfo`030?oo0?Co003oogKo00<0ool0m?l00?ooM_l00`3oo`3do`00 oomgo`030?oo0??o003oogOo00<0ool0lol00?ooMol00`3oo`3co`00oomg o`030?oo0??o003oogOo00<0ool0lol00001\ \>"], "Text", ImageSize->{620, 291}, ImageMargins->{{0, 0}, {0, 1}}, ImageRegion->{{0, 1}, {0, 1}}] }, Open ]], Cell[CellGroupData[{ Cell["Variables and Constants", "Section", Evaluatable->False, AspectRatioFixed->True], Cell[TextData[{ "Variables:\n\t", Cell[BoxData[ \(TraditionalForm\`x\_1\)]], "\ttranslation of cart 1 from center of track\n\t", Cell[BoxData[ \(TraditionalForm\`x\_2\)]], "\ttranslation of cart 2 from center of track\n\t\[Theta]\tangle of seesaw \ with vertical\n\t\[Alpha]\tangle of pendulum with normal to track\n\t", Cell[BoxData[ \(TraditionalForm\`F\_1\)]], "\tforce applied to cart 1\n\t", Cell[BoxData[ \(TraditionalForm\`F\_2\)]], "\tforce applied to cart 2" }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell[TextData[{ "Constants:\n\tJ\tinertia of seesaw with track att height h\n\t", Cell[BoxData[ \(TraditionalForm\`m\_s\)]], "\tmass of seesaw with track\n\tc\tcenter of gravity\n\th\theight of track \ from pivot point\n\t", Cell[BoxData[ \(TraditionalForm\`m\_1\)]], "\tmass of cart 1 (weight cart, on the back track)\n\t", Cell[BoxData[ \(TraditionalForm\`m\_2\)]], "\tmass of cart 2\t (pendulum cart, on the front track)\n\t", Cell[BoxData[ \(TraditionalForm\`m\_p\)]], "\tmass of pendulum\n\t", Cell[BoxData[ \(TraditionalForm\`l\_p\)]], "\tcenter of mass of pendulum rod (half of full length)\n\tg\tgravitational \ acceleration" }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(physicalvalues := {\ \ J \[Rule] 1.6, m\_s \[Rule] 6.6, c \[Rule] 0.06, h \[Rule] 0.115, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tm\_1 \[Rule] 0.48 + 0.38, m\_2 \[Rule] 0.48, m\_p \[Rule] 0.2, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\ \ l\_p \[Rule] 0.61\/2 + 0.03, g \[Rule] 9.81}\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Mathematical Modeling", "Section", Evaluatable->False, AspectRatioFixed->True], Cell["\<\ The Lagrangian methodology for deriving dynamic equations of a \ mechanical system consists of introducing so called generalized coordinates \ that describes (angular) positions and (angular) velocities of subsystems of \ the true system. Using these coordinates expressions for the kinetic and \ potential energy for each subsystem can be formed. Once this is done the \ dynamic equations describing the whole system can be derived automatically. \ \ \>", "Text"], Cell[CellGroupData[{ Cell["Computation of the Lagrangian", "Subsection", Evaluatable->False, AspectRatioFixed->True], Cell["Coordinates of the center of mass of the seesaw", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(x\_s := c\ Sin[\[Theta][t]]\), "\n", \(y\_s := c\ Cos[\[Theta][t]]\)}], "Input"], Cell["The potential and kinetic energies of the seesaw", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(V\_s := m\_s\ g\ y\_s\), "\n", \(T\_s := Simplify[1\/2\ J\ \((\[PartialD]\_t \[Theta][t])\)\^2]\)}], "Input"], Cell["Coordinates of the center of track", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(x\_c := h\ Sin[\[Theta][t]]\), "\n", \(y\_c := h\ Cos[\[Theta][t]]\)}], "Input"], Cell["Coordinates of cart 1", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(x\_m1 := x\_c + x\_1[t]\ Cos[\[Theta][t]]\), "\n", \(y\_m1 := y\_c - x\_1[t]\ Sin[\[Theta][t]]\)}], "Input"], Cell["The potential and kinetic energies of cart 1", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(V\_m1 := m\_1\ g\ y\_m1\), "\n", \(T\_m1 := Simplify[1\/2\ m\_1\ \((\((\[PartialD]\_t x\_m1)\)\^2 + \ \((\[PartialD]\_t y\_m1)\)\^2)\)]\)}], "Input"], Cell["Coordinates of cart 2", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(x\_m2 := x\_c + x\_2[t]\ Cos[\[Theta][t]]\), "\n", \(y\_m2 := y\_c - x\_2[t]\ Sin[\[Theta][t]]\)}], "Input"], Cell["The potential and kinetic energies of cart 2", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(V\_m2 := m\_2\ g\ y\_m2\), "\n", \(T\_m2 := Simplify[1\/2\ m\_2\ \((\((\[PartialD]\_t x\_m2)\)\^2 + \ \((\[PartialD]\_t y\_m2)\)\^2)\)]\)}], "Input"], Cell["Coordinates of the center of mass of the pendulum", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(x\_p := x\_m2 + l\_p\ Sin[\[Alpha][t] + \[Theta][t]]\), "\n", \(y\_p := y\_m2 + l\_p\ Cos[\[Alpha][t] + \[Theta][t]]\)}], "Input"], Cell["The potential and kinetic energies of the pendulum", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(V\_p := m\_p\ g\ y\_p\), "\n", \(T\_p := Simplify[1\/2\ m\_p\ \((\((\[PartialD]\_t x\_p)\)\^2 + \ \((\[PartialD]\_t y\_p)\)\^2)\)]\)}], "Input"], Cell["The total potential energy", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(V\_tot := V\_s + V\_m1 + V\_m2 + V\_p\)], "Input"], Cell["The total kinetic energy", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(T\_tot := T\_s + T\_m1 + T\_m2 + T\_p\)], "Input"], Cell["The Lagrangian of the system ", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(L := T\_tot - V\_tot\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Computation of the Equations of Motion", "Subsection", Evaluatable->False, AspectRatioFixed->True], Cell[TextData[{ "The equations of motion are given by\n\t ", Cell[BoxData[ FormBox[ StyleBox[\(\(d\/\(d\ t\)\) \[PartialD]L\/\[PartialD]\(q\& . \)\_i - \ \[PartialD]L\/\[PartialD]q\_i \[Equal] Q\_i\), FontSize->18], TraditionalForm]]], "\n", "(This computation may take a while ...)" }], "Text"], Cell[BoxData[{ \(eqn\_1 = Simplify[\[PartialD]\_t\((\[PartialD]\_\(\[PartialD]\_t x\_1[t]\)L)\) - \ \[PartialD]\_\(x\_1[t]\)L == F\_1]; \), "\n", \(eqn\_2 = Simplify[\[PartialD]\_t\((\[PartialD]\_\(\[PartialD]\_t x\_2[t]\)L)\) - \ \[PartialD]\_\(x\_2[t]\)L == F\_2]; \), "\n", \(eqn\_3 = Simplify[\[PartialD]\_t\((\[PartialD]\_\(\[PartialD]\_t \[Theta][t]\)L)\ \) - \[PartialD]\_\(\[Theta][t]\)L == 0]; \), "\n", \(eqn\_4 = Simplify[\[PartialD]\_t\((\[PartialD]\_\(\[PartialD]\_t \[Alpha][t]\)L)\ \) - \[PartialD]\_\(\[Alpha][t]\)L == 0]; \)}], "Input"], Cell[TextData[{ "These equations are coupled ordinary differential equations (ODEs) of \ second order in the generalized coordinates ", Cell[BoxData[ \(TraditionalForm\`x\_1, \ \[Theta], \ x\_2, \ \[Alpha]\)]], ". To rewrite these in standard state space form (a system of first order \ ODEs) we introduce the following states:\n\t", Cell[BoxData[ FormBox[ RowBox[{"x", "=", RowBox[{"(", GridBox[{ {\(x\_1\), "\[Theta]", \(x\_2\), "\[Alpha]", \(\(x\& . \)\_1\), \(\[Theta]\& . \), \(\(x\& . \ \)\_2\), \(\[Alpha]\& . \)} }], ")"}]}], TraditionalForm]], FontSize->14], ".\nThe inputs to the system are\n\t", Cell[BoxData[ FormBox[ RowBox[{"u", "=", RowBox[{"(", GridBox[{ {\(F\_1\), \(F\_2\)} }], ")"}]}], TraditionalForm]]], ".\nWe observe that the second order time derivatives of the positions and \ angles appears linearly which makes it easy to solve for these in terms the \ states." }], "Text"], Cell["Solve for second order derivatives", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(sol = Solve[\ {\ eqn\_1, eqn\_2, eqn\_3, eqn\_4\ }, \n\t\t\t\t\t\t\t\t\t\t\t\t{\ \[PartialD]\_{t, 2}x\_1[ t], \[PartialD]\_{t, 2}x\_2[ t], \n\t\t\t\t\t\t\t\t\t\t\t\t\t\ \[PartialD]\_{t, 2}\[Theta][ t], \[PartialD]\_{t, 2}\[Alpha][t]\ }\ ]; \)], "Input"], Cell["\<\ A look at the solutions (these will become the right hand sides of \ the differential equations used for simulating the system)\ \>", "Text"], Cell[BoxData[ \(sol\)], "Input"], Cell["Turn off the spell checking", "Text"], Cell[BoxData[ \(Off[General::"\"]\)], "Input"], Cell["Conversion rules to get rid of explicit time dependence", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(ssvariables := {\ \ x\_1[t] \[Rule] x\_1, \[Theta][ t] \[Rule] \[Theta], \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tx\_2[ t] \[Rule] x\_2, \[Alpha][ t] \[Rule] \[Alpha], \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\ \[PartialD]\_t x\_1[t] \[Rule] dx\_1, \[PartialD]\_t \[Theta][t] \[Rule] d\[Theta], \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\[PartialD]\_t x\_2[ t] \[Rule] dx\_2, \[PartialD]\_t \[Alpha][t] \[Rule] d\[Alpha]\ }; \)], "Input"], Cell[TextData[{ "We compute the right hand side of the nonlinear state space equations:\t", Cell[BoxData[ FormBox[ StyleBox[\(\(x\& . \)[t] = f\ [x[t], u[t]]\), FontSize->16], TraditionalForm]]] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ RowBox[{ RowBox[{"f", "=", RowBox[{ RowBox[{ RowBox[{ RowBox[{"(", GridBox[{ {\(\[PartialD]\_t\ x\_1[t]\)}, {\(\[PartialD]\_t \[Theta][t]\)}, {\(\[PartialD]\_t x\_2[t]\)}, {\(\[PartialD]\_t \[Alpha][t]\)}, {\(\[PartialD]\_{t, 2}x\_1[t]\)}, {\(\[PartialD]\_{t, 2}\[Theta][t]\)}, {\(\[PartialD]\_{t, 2}x\_2[t]\)}, {\(\[PartialD]\_{t, 2}\[Alpha][t]\)} }], ")"}], "/.", \(sol\[LeftDoubleBracket]1\[RightDoubleBracket]\)}], "/.", "ssvariables"}], "//", "Flatten"}]}], ";"}]], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["\<\ The Linearized Model (ONLY if Control Systems Professional IS \ installed)\ \>", "Subsection", Evaluatable->False, AspectRatioFixed->True], Cell["\<\ Load the Control System Professional (CSP) package (if you don't \ have this package skip this subsection and evaluate the cells in \"The \ Linearized Model (Non Control Systems Professional Users)\" instead)\ \>", \ "Text"], Cell[BoxData[ \(If[CSP, Needs["\"]]; \)], "Input"], Cell[TextData[{ "Linearize the equations around the origin\n\t", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{"x", "=", RowBox[{"(", GridBox[{ {"0", "0", "0", "0", "0", "0", "0", "0"} }], ")"}]}], ",", " ", RowBox[{"u", "=", RowBox[{ RowBox[{"(", GridBox[{ {"0", "0"} }], ")"}], " ", "\[Implies]", " ", \(\(x\& . \)[t] \[Equal] A\ x[t] + \ B\ u[t]\)}]}]}], TraditionalForm]], FontSize->14], StyleBox[".", FontSize->14], "\nusing the Linearize function in CSP" }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(If[CSP, ss = Linearize[ f, {x\_1, \[Theta], x\_2, \[Alpha]}, \n\t\t{{x\_1, 0}, {\[Theta], 0}, {x\_2, 0}, {\[Alpha], 0}, {dx\_1, 0}, {d\[Theta], 0}, {dx\_2, 0}, {d\[Alpha], 0}}, \n\t\t{{F\_1, 0}, {F\_2, 0}}] // Simplify; \nMap[MatrixForm, ss]]\)], "Input"], Cell["\<\ We extract the matrices from the linearization and substitute with \ the physical values of the parameters\ \>", "Text"], Cell[BoxData[ \(If[CSP, {\[ScriptCapitalA], \[ScriptCapitalB], \[ScriptCapitalC], \ \[ScriptCapitalD]} = \(ss /. StateSpace \[Rule] List\) /. physicalvalues]; \)], "Input"], Cell[BoxData[ \(\[ScriptCapitalA] // MatrixForm\)], "Input"], Cell[BoxData[ \(\[ScriptCapitalB] // MatrixForm\)], "Input"], Cell[BoxData[ \(\[ScriptCapitalC] // MatrixForm\)], "Input"], Cell[BoxData[ \(\[ScriptCapitalD] // MatrixForm\)], "Input"], Cell["Check controllability", "Text"], Cell[BoxData[ \(If[CSP, Controllable[ss /. physicalvalues]]\)], "Input"], Cell["\<\ Compute a full state feedback controller using the LQG-method \ (LQRegulatorGains is a CSP command).\ \>", "Text"], Cell[BoxData[ \(If[CSP, \[ScriptCapitalL] = LQRegulatorGains[\n\t\t\t\t\t\tss /. physicalvalues, \n\t\t\t\t\t\tDiagonalMatrix[{1, 5. , 1, 5. , 0, 0, 0, 0}], \n\t\t\t\t\t\t0.1\ IdentityMatrix[ 2]\ ]; \n\[ScriptCapitalL] // MatrixForm]\)], "Input"], Cell[TextData[{ "The control law to be used is\n\t", Cell[BoxData[ \(TraditionalForm\`u[t] = \(-\[ScriptCapitalL] . x[t]\)\)]], " \nwhere ", Cell[BoxData[ \(TraditionalForm\`x[t]\)]], " are measurements of the states. This gives the following closed loop \ system to simulate\n\t", Cell[BoxData[ \(TraditionalForm\`\(x\& . \)[t] = f\ [x[t], \(-\[ScriptCapitalL] . x[t]\)]\)]], "." }], "Text"] }, Open ]], Cell[CellGroupData[{ Cell["\<\ The Linearized Model (ONLY if Control Systems Professional is NOT \ installed)\ \>", "Subsection"], Cell["\<\ This subsection only needs to be evaluated if the application \ package Control System Professional is not available.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"\[ScriptCapitalA]", "=", RowBox[{ RowBox[{"(", GridBox[{ {"0", "0", "0", "0", "1", "0", "0", "0"}, {"0", "0", "0", "0", "0", "1", "0", "0"}, {"0", "0", "0", "0", "0", "0", "1", "0"}, {"0", "0", "0", "0", "0", "0", "0", "1"}, {\(-\(\(g\ h\ m\_1\)\/J\)\), \(g - \(c\ g\ h\ m\_s\)\/J\), \ \(-\(\(g\ h\ \((m\_2 + m\_p)\)\)\/J\)\), "0", "0", "0", "0", "0"}, {\(\(g\ m\_1\)\/J\), \(\(c\ g\ m\_s\)\/J\), \(\(g\ \((m\_2 + m\_p)\)\)\/J\), "0", "0", "0", "0", "0"}, {\(-\(\(g\ h\ m\_1\)\/J\)\), \(g - \(c\ g\ h\ m\_s\)\/J\), \ \(-\(\(g\ h\ \((m\_2 + m\_p)\)\)\/J\)\), \(-\(\(g\ m\_p\)\/m\_2\)\), "0", "0", "0", "0"}, {\(-\(\(g\ m\_1\)\/J\)\), \(-\(\(c\ g\ m\_s\)\/J\)\), \ \(-\(\(g\ \((m\_2 + m\_p)\)\)\/J\)\), \(\(g\ \((m\_2 + m\_p)\)\)\/\(l\_p\ m\_2\)\), "0", "0", "0", "0"} }], ")"}], "/.", "physicalvalues"}]}], ";"}]], "Input"], Cell[BoxData[ RowBox[{ RowBox[{"\[ScriptCapitalB]", "=", RowBox[{ RowBox[{"(", GridBox[{ {"0", "0"}, {"0", "0"}, {"0", "0"}, {"0", "0"}, {\(h\^2\/J + 1\/m\_1\), \(h\^2\/J\)}, {\(-\(h\/J\)\), \(-\(h\/J\)\)}, {\(h\^2\/J\), \(h\^2\/J + 1\/m\_2\)}, {\(h\/J\), \(h\/J - 1\/\(l\_p\ m\_2\)\)} }], ")"}], "/.", "physicalvalues"}]}], ";"}]], "Input"], Cell[BoxData[ RowBox[{ RowBox[{"\[ScriptCapitalC]", "=", TagBox[ RowBox[{"(", GridBox[{ {"1", "0", "0", "0", "0", "0", "0", "0"}, {"0", "1", "0", "0", "0", "0", "0", "0"}, {"0", "0", "1", "0", "0", "0", "0", "0"}, {"0", "0", "0", "1", "0", "0", "0", "0"} }], ")"}], (MatrixForm[ #]&)]}], ";"}]], "Input"], Cell[BoxData[ RowBox[{ RowBox[{"\[ScriptCapitalD]", "=", TagBox[ RowBox[{"(", GridBox[{ {"0", "0"}, {"0", "0"}, {"0", "0"}, {"0", "0"} }], ")"}], (MatrixForm[ #]&)]}], ";"}]], "Input"], Cell[BoxData[ RowBox[{ RowBox[{"\[ScriptCapitalL]", "=", RowBox[{"(", GridBox[{ {"26.7929176340570451`", "52.940017572506921`", "21.3676324363225722`", "11.8260757522378035`", "7.23981708831001657`", "18.3074986680684386`", "6.54027930973907345`", "2.28919536473904816`"}, {\(-7.65182809411970588`\), \(-26.7908119385335119`\), \ \(-10.0775285345229015`\), \(-27.725924950006533`\), \ \(-1.62889424735209065`\), \(-8.77816698970899089`\), \(-5.38840857578496201`\ \), \(-4.54939001709865742`\)} }], ")"}]}], ";"}]], "Input"], Cell["The following standard package is needed later on.", "Text"], Cell[BoxData[ \(Needs["\"]\)], "Input"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Simulation", "Section", Evaluatable->False, AspectRatioFixed->True], Cell[CellGroupData[{ Cell[TextData[{ "Load ", StyleBox["MathCode C++", FontSlant->"Italic"], " and set Current Directory" }], "Subsection"], Cell[TextData[{ "Load the ", StyleBox["MathCode C++", FontSlant->"Italic"], " package" }], "Text"], Cell["Needs[\"MathCode`\"];", "Input"], Cell[TextData[{ "Since a number of temporary files will be generated by ", StyleBox["MathCode C++", FontSlant->"Italic"], " we create a temporary directory and set this to the current working \ directory." }], "Text"], Cell["SetDirectory[$MCRoot<>\"/Demos/Pendulum\"];", "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Preliminaries", "Subsection"], Cell["We store the states as a vector", "Text"], Cell[BoxData[ \(x\_ss := {x\_1, \[Theta], x\_2, \[Alpha], dx\_1, d\[Theta], dx\_2, d\[Alpha]}\)], "Input"], Cell["\<\ and close the feedback loop which gives the following right hand \ side of the state space equations\ \>", "Text"], Cell[BoxData[ \(f1 = f /. Thread[{F\_1, F\_2} \[Rule] \(-\[ScriptCapitalL] . x\_ss\)]; \)], "Input"], Cell["Turn off the spell checking", "Text"], Cell[BoxData[ \(Off[General::"\"]\)], "Input"], Cell["\<\ Change names of variables to only ASCII characters (needed for code \ generation)\ \>", "Text"], Cell[BoxData[ \(Cvariables\ = \n\t\t{x\_1 \[Rule] \ x1, \[Theta]\ \[Rule] \ theta, \n\t\t\tx\_2\ \[Rule] \ x2, \[Alpha]\ \[Rule] \ alpha, \n\t\t\tdx\_1\ \[Rule] \ dx1, d\[Theta]\ \[Rule] \ dtheta, \n\t\t\tdx\_2\ \[Rule] \ dx2, \ d\[Alpha] \[Rule] \ dalpha}; \)], "Input"], Cell["\<\ The right hand side of the closed loop state space equations \ (intended for code generation)\ \>", "Text"], Cell[BoxData[ \(rhs = \(f1 /. physicalvalues\)\ /. \ Cvariables; \)], "Input"], Cell[CellGroupData[{ Cell["The right hand side of the state space equations", "Subsubsection"], Cell["\<\ We define the function ClosedLoopPendulumEquationsRHS which has the \ huge rhs expression from above as its body. Note that the expression also can \ read from the file PendulumExpr.m which is created above (in case of memory \ problems). This file makes it quicker to try out the compilation in different \ sessions without recomputing rhs each time.\ \>", "Text"], Cell[BoxData[ \(ClosedLoopPendulumEquationsRHS[\n\tReal\ x1_, Real\ theta_, Real\ x2_, Real\ alpha_, \n\tReal\ dx1_, Real\ dtheta_, Real\ dx2_, Real\ dalpha_] \[Rule] Real[8] = rhs; \)], "Input"], Cell["\<\ Note that we have added type declarations to be able to generate \ code for this function later on.\ \>", "Text"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[TextData[{ "Simulation of the Nonlinear Model within ", StyleBox["Mathematica", FontSlant->"Italic"] }], "Subsection"], Cell[TextData[{ "We define ", Cell[BoxData[ \(TraditionalForm\`f\_cl[t]\)]], " to be a short notation for the ClosedLoopPendulumEquationsRHS" }], "Text"], Cell[BoxData[ \(f\_cl[t_] := ClosedLoopPendulumEquationsRHS[x1[t], theta[t], x2[t], alpha[t], dx1[t], dtheta[t], dx2[t], dalpha[t]]\)], "Input"], Cell["The state vector", "Text"], Cell[BoxData[ \(x[t_] := {x1[t], theta[t], x2[t], alpha[t], dx1[t], dtheta[t], dx2[t], dalpha[t]}; \)], "Input"], Cell["The differential equations for the closed loop system", "Text"], Cell[BoxData[ \(deq := Thread[\(x'\)[t] == f\_cl[t]]\)], "Input"], Cell["\<\ Initial conditions for the simulation. Cart 1 is positioned at -1, \ the seesaw is horizontal, cart 2 is positioned at +1 and the pendulum has a \ small deviation of 5.7\[Degree] from the vertical axis. The system is at rest \ at these positions (all time derivatives are zero).\ \>", "Text"], Cell[BoxData[ \(initeq := Thread[x[0] == {\(-1\), 0, 1, 0.1, 0, 0, 0, 0}]\)], "Input"], Cell["We simulate the system using NDSolve", "Text"], Cell[TextData[StyleBox["Standard usage", FontWeight->"Bold"]], "Text"], Cell[BoxData[{ \(\(dsolStandard = NDSolve[\ deq\ \[Union] \ initeq, {x1, theta, x2, alpha, dx1, dtheta, dx2, dalpha}, {t, 0, 10}\ ]; \) // Timing\), "\n", \(NDSolveStandardTime = First[%]; \)}], "Input"], Cell[TextData[StyleBox["Uncompiled equations", FontWeight->"Bold"]], "Text"], Cell[BoxData[{ \(\(dsolUncompiled = NDSolve[\ deq\ \[Union] \ initeq, {x1, theta, x2, alpha, dx1, dtheta, dx2, dalpha}, {t, 0, 10}, Compiled \[Rule] False\ ]; \) // Timing\), "\n", \(NDSolveUncompiledTime = First[%]; \)}], "Input"], Cell[TextData[{ StyleBox["Force ", FontWeight->"Bold"], StyleBox["Mathematica", FontWeight->"Bold", FontSlant->"Italic"], StyleBox[" to use 4-5 Runge-Kutta with step size 0.01", FontWeight->"Bold"] }], "Text"], Cell[BoxData[{ \(\(dsolRK = NDSolve[\ deq\ \[Union] \ initeq, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t{x1, theta, x2, alpha, dx1, dtheta, dx2, dalpha}, {t, 0. , 10}, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tMaxStepSize \[Rule] 0.01, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tMethod \[Rule] RungeKutta, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tMaxSteps \[Rule] 10000\ ]; \) // Timing\), "\n", \(NDSolveRKTime = First[%]; \)}], "Input"], Cell[TextData[{ StyleBox["Force ", FontWeight->"Bold"], StyleBox["Mathematica", FontWeight->"Bold", FontSlant->"Italic"], StyleBox[" to use 4-5 Runge-Kutta with step size 0.01 (uncompiled \ equations)", FontWeight->"Bold"] }], "Text"], Cell[BoxData[{ \(\(dsolRKUncompiled = NDSolve[\ deq\ \[Union] \ initeq, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t{x1, theta, x2, alpha, dx1, dtheta, dx2, dalpha}, {t, 0. , 10}, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tMaxStepSize \[Rule] 0.01, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tMethod \[Rule] RungeKutta, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tMaxSteps \[Rule] 10000, \n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tCompiled \[Rule] False\ ]; \) // Timing\), "\n", \(NDSolveRKUncompiledTime = First[%]; \)}], "Input"], Cell[TextData[StyleBox["The linear system", FontWeight->"Bold"]], "Text"], Cell[BoxData[ \(dlineq := Thread[\(x'\)[ t] == \((\[ScriptCapitalA] - \[ScriptCapitalB] . \ \[ScriptCapitalL])\) . x[t]] // Chop\)], "Input"], Cell["We simulate the linear system using NDSolve", "Text"], Cell[BoxData[ \(dsollin = NDSolve[\ dlineq\ \[Union] \ initeq, {x1, theta, x2, alpha, dx1, dtheta, dx2, dalpha}, {t, 0, 8}\ ]; \)], "Input"], Cell[TextData[{ "Plot the result for the first 4 state variables (the ", StyleBox["Mathematica", FontSlant->"Italic"], " standard usage of NDSolve) together with the result for the linear \ system. The other simulation methods above can be shown to give the same \ result as the standard usage of NDSolve." }], "Text"], Cell[BoxData[ \(Map[\n\t Plot[#, {t, 0, 8}, PlotRange \[Rule] All, DisplayFunction \[Rule] Identity, PlotStyle \[Rule] {{RGBColor[0, 0, 0]}, {RGBColor[1, 0, 0]}}] &, {Part[Flatten[x[t] /. dsolStandard], Range[4]], \n\t\tPart[Flatten[x[t] /. dsollin], Range[4]]} // Transpose\ ]; \)], "Input"], Cell[BoxData[ \(\(Show[GraphicsArray[%], DisplayFunction \[Rule] $DisplayFunction];\)\)], "Input"], Cell["\<\ We notice that the linear response (red) differs significantly from \ the nonlinear one.\ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["External Simulation of the Nonlinear Model", "Subsection"], Cell[CellGroupData[{ Cell["Code for External Simulation of the Nonlinear Model", "Subsubsection"], Cell["\<\ We first make an implementation of a state equation solver using \ the Runge-Kutta method\ \>", "Text"], Cell[BoxData[ \(RK[Integer\ n_, Real\ h_, Real\ t0_, Real[_]\ startv_, Integer\ dimen_, Integer\ dimPlusOne_] \[Rule] Real[n, dimPlusOne] := \n\t\t\tModule[{\n\t\t\t\t\tReal[dimen]\ {x, k1, k2, k3, k4}, \n\t\t\t\t\tInteger\ i, \n\t\t\t\t\tReal\ t, \n\t\t\t\t\t\ Real[n, dimPlusOne]\ res\ \n\t\t\t\t}, \n\t\t\t\tt = t0; \n\t\t\t\tx = startv; \n\t\t\t\tres\[LeftDoubleBracket]1, 1\[RightDoubleBracket] = t; \n\t\t\t\tres\[LeftDoubleBracket]1, 2 | dimPlusOne\[RightDoubleBracket] = x; \n\t\t\t\tFor[i = 1, \ i \[LessEqual] n - 1, \ \(i++\), \n\t\t\t\t\tk1 = h\ ClosedLoopPendulumEquationsRHS[ x\[LeftDoubleBracket]1\[RightDoubleBracket], x\[LeftDoubleBracket]2\[RightDoubleBracket], x\[LeftDoubleBracket]3\[RightDoubleBracket], x\[LeftDoubleBracket]4\[RightDoubleBracket], x\[LeftDoubleBracket]5\[RightDoubleBracket], x\[LeftDoubleBracket]6\[RightDoubleBracket], x\[LeftDoubleBracket]7\[RightDoubleBracket], x\[LeftDoubleBracket]8\[RightDoubleBracket]]; \n\t\t\t\t\tk2 \ = h\ ClosedLoopPendulumEquationsRHS[ x\[LeftDoubleBracket]1\[RightDoubleBracket] + k1\[LeftDoubleBracket]1\[RightDoubleBracket]/2, x\[LeftDoubleBracket]2\[RightDoubleBracket] + k1\[LeftDoubleBracket]2\[RightDoubleBracket]/2, x\[LeftDoubleBracket]3\[RightDoubleBracket] + k1\[LeftDoubleBracket]3\[RightDoubleBracket]/2, x\[LeftDoubleBracket]4\[RightDoubleBracket] + k1\[LeftDoubleBracket]4\[RightDoubleBracket]/2, x\[LeftDoubleBracket]5\[RightDoubleBracket] + k1\[LeftDoubleBracket]5\[RightDoubleBracket]/2, x\[LeftDoubleBracket]6\[RightDoubleBracket] + k1\[LeftDoubleBracket]6\[RightDoubleBracket]/2, x\[LeftDoubleBracket]7\[RightDoubleBracket] + k1\[LeftDoubleBracket]7\[RightDoubleBracket]/2, x\[LeftDoubleBracket]8\[RightDoubleBracket] + k1\[LeftDoubleBracket]8\[RightDoubleBracket]/ 2]; \n\t\t\t\t\tk3 = h\ ClosedLoopPendulumEquationsRHS[ x\[LeftDoubleBracket]1\[RightDoubleBracket] + k2\[LeftDoubleBracket]1\[RightDoubleBracket]/2, x\[LeftDoubleBracket]2\[RightDoubleBracket] + k2\[LeftDoubleBracket]2\[RightDoubleBracket]/2, x\[LeftDoubleBracket]3\[RightDoubleBracket] + k2\[LeftDoubleBracket]3\[RightDoubleBracket]/2, x\[LeftDoubleBracket]4\[RightDoubleBracket] + k2\[LeftDoubleBracket]4\[RightDoubleBracket]/2, x\[LeftDoubleBracket]5\[RightDoubleBracket] + k2\[LeftDoubleBracket]5\[RightDoubleBracket]/2, x\[LeftDoubleBracket]6\[RightDoubleBracket] + k2\[LeftDoubleBracket]6\[RightDoubleBracket]/2, x\[LeftDoubleBracket]7\[RightDoubleBracket] + k2\[LeftDoubleBracket]7\[RightDoubleBracket]/2, x\[LeftDoubleBracket]8\[RightDoubleBracket] + k2\[LeftDoubleBracket]8\[RightDoubleBracket]/ 2]; \n\t\t\t\t\tk4 = h\ ClosedLoopPendulumEquationsRHS[\n\t\t\t\t\t\tx\ \[LeftDoubleBracket]1\[RightDoubleBracket] + k3\[LeftDoubleBracket]1\[RightDoubleBracket], x\[LeftDoubleBracket]2\[RightDoubleBracket] + k3\[LeftDoubleBracket]2\[RightDoubleBracket], x\[LeftDoubleBracket]3\[RightDoubleBracket] + k3\[LeftDoubleBracket]3\[RightDoubleBracket], x\[LeftDoubleBracket]4\[RightDoubleBracket] + k3\[LeftDoubleBracket]4\[RightDoubleBracket], x\[LeftDoubleBracket]5\[RightDoubleBracket] + k3\[LeftDoubleBracket]5\[RightDoubleBracket], x\[LeftDoubleBracket]6\[RightDoubleBracket] + k3\[LeftDoubleBracket]6\[RightDoubleBracket], x\[LeftDoubleBracket]7\[RightDoubleBracket] + k3\[LeftDoubleBracket]7\[RightDoubleBracket], x\[LeftDoubleBracket]8\[RightDoubleBracket] + k3\[LeftDoubleBracket]8\[RightDoubleBracket]]; \n\t\t\t\t\t\ x = x + 1/6\ \((k1 + 2\ k2 + 2 k3 + k4)\); \n\t\t\t\t\tt = t + h; \n\t\t\t\t\tres\[LeftDoubleBracket]i + 1, 1\[RightDoubleBracket] = t; \n\t\t\t\t\tres\[LeftDoubleBracket]i + 1, 2 | dimPlusOne\[RightDoubleBracket] = x; \n\t\t\t\t]; \n\t\t\t\tres\n\t\t\t]\)], "Input"], Cell["We use RK in the following ODE solver", "Text"], Cell[BoxData[ \(ndSolve[initvalue_, maxstepsize_, {x_, xmin_, xmax_}] := \n\t Module[{n, dimen}, \n\t\tn = \ IntegerPart[\((xmax - xmin)\)/maxstepsize] + 1; \n\t\tdimen = Length[initvalue]; \n\t\tres = RK[n, maxstepsize, xmin, initvalue, dimen, dimen + 1]; \n\t\tArray[\((Interpolation[ AppendRows[ Transpose[{res\[LeftDoubleBracket]_, 1\[RightDoubleBracket]}], \n\t\t\t\t\t\t\t\t\t\ Transpose[{res\[LeftDoubleBracket]_, # + 1\[RightDoubleBracket]}]]])\) &, \ {8}]\n\t]\)], \ "Input"], Cell["\<\ To get a correct estimate of the time spent by the external code \ for simulating the system we define the folowing function which repeats the \ external simulation n times.\ \>", "Text"], Cell[BoxData[ \(RKRepeated[Integer\ n_, Real\ h_, Real\ t0_, Real[_]\ startv_, Integer\ dimen_, Integer\ dimPlusOne_, Integer\ loops_] \[Rule] Real[n, dimPlusOne] := \n\t\ Module[{Real[n, dimPlusOne]\ res}, \ \n\t\tDo[\t res\ = \ RK[n, h, t0, startv, dimen, dimPlusOne], {loops}]; \n\t\tres\n\t]\)], "Input"], Cell[BoxData[ \(ndSolveRepeated[\n\tinitvalue_, maxstepsize_, {x_, xmin_, xmax_}, loops_] := \n\t Module[{n, dimen}, \n\t\tn = \ IntegerPart[\((xmax - xmin)\)/maxstepsize] + 1; \n\t\tdimen = Length[initvalue]; \n\t\tres = RKRepeated[n, maxstepsize, xmin, initvalue, dimen, dimen + 1, \ loops]; \n\t\tArray[\((Interpolation[ AppendRows[ Transpose[{res\[LeftDoubleBracket]_, 1\[RightDoubleBracket]}], \n\t\t\t\t\t\t\t\t\t\ Transpose[{res\[LeftDoubleBracket]_, # + 1\[RightDoubleBracket]}]]])\) &, \ {8}]\n\t]\)], \ "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Compilation", "Subsubsection"], Cell["\<\ We are now ready to compile the package. The rhs of the state \ equations are optimized using common subexpression elimination\ \>", "Text"], Cell[BoxData[ \(CompilePackage[ EvaluateFunctions \[Rule] {ClosedLoopPendulumEquationsRHS}]\)], "Input"], Cell["The corresponding code", "Text"], Cell[BoxData[ \(\(!! Global.cc\)\)], "Input"], Cell[TextData[{ "All necessary files for building binaries are now generated. The \ MakeBinary command can build either a standalone application using the code \ specified by the option GenerateMainFileAndFunction or a ", StyleBox["MathLink", FontSlant->"Italic"], " version which can be easily installed into ", StyleBox["Mathematica", FontSlant->"Italic"], ". We build the ", StyleBox["MathLink", FontSlant->"Italic"], " version" }], "Text"], Cell[BoxData[ \(MakeBinary[]; \)], "Input"], Cell["\<\ To use the generated code for the simulation we only have to \ install the code using the Install command.\ \>", "Text"] }, Open ]], Cell[CellGroupData[{ Cell["External Simulation", "Subsubsection"], Cell["\<\ This function is a fix to be able to measure the time for external \ processes.\ \>", "Text"], Cell[BoxData[ \(SetAttributes[AbsTime, HoldFirst]; AbsTime[x_] := Module[{start, res}, \n\t\tstart = AbsoluteTime[]; res = x; {\((AbsoluteTime[] - start)\)\ Second, res}]; \)], "Input"], Cell[TextData[{ StyleBox["Uncompiled code (computations within ", FontWeight->"Bold"], StyleBox["Mathematica", FontWeight->"Bold", FontSlant->"Italic"], StyleBox[")", FontWeight->"Bold"] }], "Text"], Cell[BoxData[{ \(\(dsolUncompiled = ndSolve[{\(-1\), 0, 1, 0.1, 0, 0, 0, 0}, 0.01, {t, 0. , 10. }]; \) // Timing\), "\n", \(ndSolveUncompiledTime = First[%]; \)}], "Input"], Cell["We plot the result to verify the Runge-Kutta solver", "Text"], Cell[BoxData[ \(Map[ Plot[#[t], {t, 0, 8}, PlotRange \[Rule] All, DisplayFunction \[Rule] Identity] &, Part[dsolUncompiled, Range[4]]]; \)], "Input"], Cell[BoxData[ \(Show[GraphicsArray[%], DisplayFunction \[Rule] $DisplayFunction]; \)], "Input"], Cell[TextData[StyleBox["Compiled code (external computations)", FontWeight->"Bold"]], "Text"], Cell["We install the executable code", "Text"], Cell[BoxData[ \(InstallCode[]; \)], "Input"], Cell["\<\ The external simulation is performed 600 times to get a accurate \ measure of its timing\ \>", "Text"], Cell[BoxData[{ \(\(dsolCompiled = ndSolveRepeated[{\(-1\), 0, 1, 0.1, 0, 0, 0, 0}, 0.01, {t, 0. , 10. }, 600]; \) // AbsTime\), "\n", \(ndSolveCompiledTime = First[%]/600; \)}], "Input"], Cell["We plot the result to verify the compiled Runge-Kutta solver", "Text"], Cell[BoxData[ \(Map[ Plot[#[t], {t, 0, 8}, PlotRange \[Rule] All, DisplayFunction \[Rule] Identity] &, Part[dsolCompiled, Range[4]]]; \)], "Input"], Cell[BoxData[ \(Show[GraphicsArray[%], DisplayFunction \[Rule] $DisplayFunction]; \)], "Input"], Cell["Looks good!", "Text"] }, Open ]] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Comparing Performance", "Section"], Cell["\<\ We present the comparision of performance in the following \ table\ \>", "Text"], Cell[BoxData[ RowBox[{"TableForm", "[", RowBox[{ RowBox[{"(", GridBox[{ { "NDSolveStandardTime", \(NDSolveStandardTime/ ndSolveCompiledTime\)}, { "NDSolveUncompiledTime", \(NDSolveUncompiledTime/ ndSolveCompiledTime\)}, {"NDSolveRKTime", \(NDSolveRKTime/ndSolveCompiledTime\)}, { "NDSolveRKUncompiledTime", \(NDSolveRKUncompiledTime/ ndSolveCompiledTime\)}, { "ndSolveUncompiledTime", \(ndSolveUncompiledTime/ ndSolveCompiledTime\)}, { "ndSolveCompiledTime", \(ndSolveCompiledTime/ ndSolveCompiledTime\)} }], ")"}], ",", \(TableHeadings \[Rule] {{"\", "\ \", "\", "\", \ "\", "\"}, {"\