Maxima/Introduction By Example

Getting Used to Maxima By Way of an Example of Use The examples that follows are presented for the purpose of providing some beginning perspective, and hopefully some motivation to make the effort to get familiar with how Maxima works. The task at hand is relevant in the context of aerodynamics and aviation. The specific goal is to find the dynamic pressure at a true airspeed (VTAS) of 200 ft/s at sea level. Dynamic pressure is calculated using the expression:Where: = dynamic pressure lbf/ft² (pounds of force per square foot) = air density (at sea level is 0.0023769) slugs/ft³ = true airspeed in feet per second (ft/s). British Engineering units are used in this context to avoid multiple conversions due to the units used. It is better to convert in the end if necessary in this case. Note: The appearance of the following Maxima examples may vary depending on Unicode support in the version of Maxima being used and the Maxima user interface -- UI. The Maxima examples used in this book are from the wxMaxima GUI unless noted otherwise, and a version of Maxima with Unicode support. First, enter the Maxima expression for dynamic pressure q:(%i1) q : 1/2*\ρ*V^2; (q) (V^2*ρ)/2This Maxima expression uses the colon character to assign the expression for dynamic pressure to the identifier q. The identifier q is now simply a name for an expression. The equal sign ( = ) is not used for this operation in Maxima as is the case with some programming languages. This Maxima line of a input expression ends with a semicolon ( ; ) character. Each line of input must end with a semicolon or the dollar sign character ( $ ), the use of which will be described later. Next, enter an assignment expression for the numerical value of the air density at sea level:(%i2) ρ : 0.0023769; (ρ) 0.0023769Density ρ is in the units of slugs per cubic foot - slugs/ft³. With this next input expression, Maxima is asked to evaluate q. This is accomplished by two single quotation marks placed before q as shown next:(%i3) ''q; (%o3) 0.0033855446120918224*VTAS^2 The output indicates that q depends on VTAS -- true airspeed. The goal is to evaluate q, if VTAS = 200 ft/s is true. One way to do that is with the following Maxima expression:(%i4) ''q, VTAS = 200; (%o4) 135.42178448367287The result is about 135 pounds of force per square foot. There is another way to accomplish this that may be somewhat more convenient for determining sea level dynamic pressure given true airspeed. First, assign the Maxima floating point literal value of 0.0033855446120918224 to an identifier named c as follows:(%i5) c : 0.0033855446120918224; (c) 0.0033855446120918224Next, enter(%i6) q(VTAS) := c*VTAS^2; (%o6) q(VTAS):=c*VTAS^2(%i7) q(200); (%o7) 135.4217844836729This result is limited to sea level density. That does not make the result of no practical use. This relationship, or something similar is, used to calibrate aircraft airspeed indicators to display indicated airspeed (IAS) in the cockpit. However, it might be the case that it is necessary to determine the dynamic pressure at other altitudes besides sea level. In the standard atmosphere model, air density is a function of altitude (h), and is defined piece-wise based on the atmospheric layer. The density is derived from the Ideal Gas Law: where in British Engineering units: ‒ is density in slugs/ft³ ‒ is pressure in lbf/ft² ‒ is temperature in °R (Rankine) ‒ is the specific gas constant (for dry air) approximately 1716.56 (ft·lbf)/(slug·°R) (%i8) P = \ρ*R*T; (%o8) P=RTρThe backslash character ( \ ) in front of the character for in the input expression allows an identifier to use Unicode. This makes for some nicer looking output. This is a Maxima mathematical equation expression. This equation can be solved for as follows:(%i9) solve(%,ρ); (%o58) [ρ=P/(R*T)]The Maxima built-in identifier % names the last output expression.The output of this expression is a Maxima list object with the only instance being the equation expression solved for . Depending on the expression, the solve() operation may result in more than one solution expression. Maxima creates a list object for one or more solution expressions. This solution equation expression can be accessed as follows:(%i10) %[1]; (%o10) ρ=P/(R*T)The value of this expression is obtained by referencing the first (and only) solution of the Maxima list object of solutions ‒ [ρ=P*R*T]. The reference is obtained by use of the [1] index syntax as shown. There are other ways Maxima can be used to accomplish the same result: (%i11) igl : P = ρ*R*T; (igl) P=R*T*ρIn this case the expression for the ideal gas law was assigned an identifier igl. The advantage to this sort of assignment is that the equation can be referred to using the identifier as follows for example:(%i12) solve(igl,ρ); (%o12) [ρ=P/(R*T)]The identifier igl can be used at any point to refer to the equation expression: P = \ρ*R*T). An assignment operator can also be used to assign the identifier solution to an expression. For example:(%i13) solve(igl,ρ); (%o13) [ρ=P/(R*T)] (%i14) solution : solve(igl,ρ); (solution) [ρ=P/(R*T)] (%i15) solve(solution,P)[1]; (%o15) P=R*T*ρThe last input and output expression might be somewhat of a "head-scratcher" at first glance, but it is worthwhile to understand how Maxima works in this case. The identifier solution names an expression that uses the solve() operator on the identifier igl as an argument. This expression is evaluated first by Maxima. Then the solve() operator is applied to the result solving for P. The result is a list object with one solution expression. This expression is then referred to using the [1] index, and the output is the equation expression P=R*T*ρ . The expression for ρ (density) as a function of pressure (p) and temperature (°F) can be checked against known values using the following expressions: (%i16) \°R(°F) := °F + 459.67; (%o16) °R(°F):=°F+459.67 (%i17) P(p) := p*144.0; (%o17) P(p):=p*144.0 (%i18) R:1716.56; (R) 1716.56 (%i19) ρ(p,°F) := P(p)/(''R*\°R(°F)); (%o19) ρ(p,°F):=P(p)/(1716.56*°R(°F))The above Maxima expressions are examples of Maxima function expressions. For standard sea level conditions of pressure and temperature as follows, p is 14.696 lbf/ft² °F is 59.0° the calculated value for density is:(%i20) ρ(14.696,59.0); (%20) 0.0023769033250291326Compare this resulting value of density with the the value of the table below for 0 altitude and temperature of 59.0 °F: +U.S. Standard Atmosphere vs. Altitude - BE UnitsAltitude above Sea Level - h - (ft)Temperature - T - (oF)Acceleration of Gravity - g - (ft/s2)Pressure - p - (lb/in2)Density - ρ - (10-4 slugs/ft3)059.032.17414.69623.77500041.1732.15912.22820.481000023.3632.14310.10817.56150005.5532.1288.29714.9620000 -12.2632.1126.75912.6725000 -30.0532.0975.46110.6630000 -47.8332.0824.3738.9135000 -65.6132.0663.4687.38 It would be useful to have the means to round result values to a desired number of decimal places. Maxima does not have a built-in operator for rounding decimal fraction values. Also, there is no package that is included as part of the Maxima distribution that provides for this capability. However, Robert Dodier has created a Maxima package that includes a function for this purpose ‒ named excel_round.macDodier Robert maxima-packages/robert-dodier/excel_round/ https://github.com/maxima-project-on-github/maxima-packages/tree/master/robert-dodier/excel_round GitHub. However, to be able access the excel_round() function of this package, it must first be downloaded or copied into a location accessible to Maxima using the load() operator. The built-in Maxima load() operator is used to access specialized functions that are not built-in to Maxima. In order to use the load() operator, the location of the downloaded excel_round.mac file has to be part of the argument to the load() operator unless the file happens to be in the Maxima current working directory. The built-in operatingsystem package can provide some functions to help with determining the current working directory. This package is loaded as follows: (%i21) load(operatingsystem); (%o21) "/usr/local/share/maxima/5.49.0/share/contrib/operatingsystem/operatingsystem.mac" (%i22) getcurrentdirectory(); (%o22) "/home/username/"In this example username is the name of the home directory of the user. so if the downloaded excel_round.mac file is not located in the current working directory, then the directory can be changed to the location where the file is located: (%i22) chdir("/home/username/local/maxima/packages"); (%o22) "#P"/home/username/local/maxima/packages""Now the excel_round.mac file can be loaded assuming it is within the directory /home/username/local/maxima/packages as follows:(%i23) load("excel_round.mac"); (%o23) "excel_round.mac"Alternatively, if the location of the file is known it can be loaded using the full path name as follows for example:(%i24) load("/home/username/local/maxima/packages/excel_round.mac"); (%o24) "/home/username/local/maxima/packages/excel_round.mac"This method obviates the need for loading the operatingsystem package. Note: The wxMaxima GUI makes loading packages simply a matter of a menu selection. The density value above can be rounded to six decimal places as follows:(%i25) excel_round(ρ(14.696,59.0),6); (%o25) 0.002377This rounded value agrees with the sea level value of slugs/ft³ taken from the table of Standard Atmosphere values. For Standard Atmosphere values of pressure ( 6.758 lbf/in² ) and temperature ( -12.26 °F ) at 20,000 ft the density value is:(%i26) excel_round(ρ(6.759,-12.26),6); (%o26) 0.001267This rounded value agrees with the sea level value of slugs/ft³ taken from the table of Standard Atmosphere values. Pressure However, what is needed for the task at hand is density as a function of altitude. Both the pressure P and the temperature T can be expressed as a function of altitude h: Where in British Engineering units, the constants are: is Pressure in pounds per square foot lbf/ft² is (Sea Level Pressure): 2116.224 lbf/ft² or 14.696 lbf/in² is the specific gas constant (for dry air) approximately 1716.56 (ft·lbf)/(slug·°R) is (Sea Level Temperature): 518.67 °R (15 °C or 59 °F) is (Lapse Rate): 0.003563 °R/ft is 32.174 ft/s² is altitude in feet above sea level The exponent evaluates to approximately 5.26054:(%i27) 32.174/(0.003563*1716.56); (%o27) 5.260539027391964 (%i28) excel_round(%,5); (%o28) 5.26054 evaluates to approximately (%i29) 0.003563/518.67; (%o29) 6.869493126650858*10^-6 (%i30) excel_round(6.869493126650858*10^-6,11); (%o30) 6.86949*10^-6 Temperature is Temperature in °R is altitude in feet. Practical Formulas for Pressure and Temperature Thus, the practical formula for pressure and temperatures in terms of British Engineering units are: lbf/ft² °R Pressure (%i31) P(h) := 2116.224*(1-6.8649*10^-6*h)^5.26054; (%o31) P(h):=2116.224*(1-6.8649*10^(-6)*h)^5.26054 (%i32) p(P) := P/144.0; (%o32) p(P):=P/144.0 (%i33) P(20000); (%o33) 973.0940875140712 (%i34) p(%); (%o34) 6.757597829958828 (%i35) excel_round(%,3); (%o35) 6.758 Temperature (%i36) T(h) := 518.67 - 0.003563*h; (%o36) T(h):=518.67-0.003563*h (%i37) T(20000); (%o37) 447.40999999999997 (%i38) excel_round(%,2); (%o38) 447.41 (%i39) declare("°", alphabetic); (%o39) done (%i39) °F(°R) := °R -459.67; (%o39) °F(°R):=°R-459.67 (%i40) °F(447.41); (%o40) -12.259999999999991 (%i41) excel_round(%,2); (%o41) -12.26 Density The expressions for pressure and temperature as functions of altitude can be used with the expression for the ideal gas law for an expression of density as a function of altitude as follows:(%i42) ρ(h) := P(h)/(R*T(h)); (%o42) ρ(h):=P(h)/(R*T(h)) (%i42) excel_round(ρ(20000)*10^4,2); (%o42) 12.67This result agrees with the density value of the standard atmosphere at 20,000 ft. Using the expression for the air density as a function of altitude, an expression for dynamic pressure is as follows:(%i43) q(h,VTAS) := ((ρ(h)*VTAS^2))/2; (%o43) q(h,VTAS):=(ρ(h)*VTAS^2)/2Using the expression above, values of dynamic pressure can be calculated as a function of altitude and true airspeed as follows for 20,000 feet and 15,000 feet at 200 ft/s:(%i44) q(20000,200); (%o44) 25.34079027107967 (%i45) q(15000,200); (%o45) 29.922636037412556This value can be checked using the value for density from the Standard Atmosphere table at 15,000 ft and the expression for dynamic pressure as follows:(%i46) excel_round(%,2); (%o46) 29.92 (%i47) (200.0^2)*0.001496/2.0; (%o47) 29.919999999999998 (%i48) excel_round(%,2); (%o48) 29.92In the context of aviation and aircraft performance, airspeeds are commonly referenced to knots (nautical miles per hour) instead of ft/s. To determine dynamic pressure as a function of altitude and true airspeed in knots, the following expressions are used:(%i49) VTAS(knots) := 1.687811*knots; (%o49) VTAS(knots):=1.687811*knots (%i50) excel_round(VTAS(118.49676),2); (%o50) 200.0 (%i51) q(h,knots) := ((ρ(h)*VTAS(knots)^2))/2; (%o51) q(h,knots):=(ρ(h)*VTAS(knots)^2)/2 (%i52) excel_round(q(15000,118.49676),2); (%o52) 29.92So approximately 118 knots true airspeed is equivalent to 200 ft/s, and the above expression for dynamic pressure is in terms of a function of altitude in feet and true airspeed in knots. Plotting Maxima has several ways to do plotting. The image that follows is a simple plot of density as a function of altitude produced using the plotting operator wxplot2d():(%i53) T(h) := 518.688 - 0.0035616*h; (%o53) T(h):=518.688-0.0035616*h (%i54) P(h) := 2116.224*(1-6.8649*10^-6*h)^5.26054; (%o54) P(h):=2116.224*(1-6.8649*10^(-6)*h)^5.26054 (%i55) R : 1716.56; (R) 1716.56 (%i56) ρ(h) := P(h)/(R*T(h)); (%o56) ρ(h):=P(h)/(R*T(h)) (%i57) wxplot2d([ρ(h)], [h,-100,35000])$ thumb|876x876px|Air Density (slugs/ft³) As a Function of Altitude (ft)|none A Practical Application - True Airspeed (TAS) Having the display of TAS during actual flight conditions is very useful information. With the necessary sensors and electronics, TAS can be displayed continuously to pilots during changing flight conditions. One way to accomplish this is to have aircraft sensors that measure values of dynamic pressure, air pressure, and air temperature. These measured values can be used by electronic devices designed to calculate TAS for display and use by other equipment. The essentials of how this could be implemented can be demonstrated using Maxima by way of the expressions that follow:(%i1) dpe : q = 1/2*ρ*V^2; (dpe) q=(V^2*ρ)/2 (%i2) solutions:solve(dpe,V); (solutions) [V=-(sqrt(2)*sqrt(q/ρ)),V=sqrt(2)*sqrt(q/ρ)]The above expressions result in a list of two solution expressions. The second of these two is the one that is useful in this context.(%i3) solutions[2]; (%o3) V=sqrt(2)*sqrt(q/ρ) (%i4) igl : ρ = P/(R*T); (igl) ρ=P/(R*T) (%i5) eq : subst(rhs(igl),ρ,solutions[2]); (%o5) V=sqrt(2)*sqrt((R*T*q)/P)In the last of the above expressions the subst() operator is used to substitute an expression for density ρ in the expression for speed V.(%i6) R :1716.56; (R) 1716.56 (%i7) V(q,P,T) := ''(rhs(eq)); (%o7) V(q,P,T):=sqrt(2)*sqrt((R*T*q)/P)The value for R determines the units: dynamic pressure q is in lbf/ft² pressure P is in lbf/ft² T is in °R V is in ft/s The second of the two expressions above define a Maxima function for V in terms of values of dynamic pressure q, air pressure P, and temperature T. This expression for speed can be checked using known values. At 15,000 ft if: q is 29.92 lbf/ft² P is 1194.79 lbf/ft² T is 465.24 °R then:(%i8) load("excel_round.mac"); (%o8) "/home/davidmiller/Downloads/excel_round.mac" (%i9) excel_round(float(V(29.92,1194.79,465.24)),1); (%o9) 200.0This resulting value of 200 ft/s for the TAS agrees with previous calculations under the given conditions. However, this is TAS in ft/s. Displaying TAS in terms of knots is much more relevant. The conversion factor for this purpose is 1.687811. So TAS in terms of knots can be expressed as follows:(%i10) TAS(q,P,T) := V(q,P,T)/1.687811; (%o10) TAS(q,P,T):=V(q,P,T)/1.687811 (%i11) excel_round(float(TAS(29.92,1194.79,465.24)),1); (%o11) 118.5The result is 118.5 knots for TAS which agrees with previous calculated values. The expression for TAS could have been more general by using an identifier (cf, for example) for the conversion factor as follows:(%i12) TAS(q,P,T) := V(q,P,T)/cf; (%o12) TAS(q,P,T):=V(q,P,T)/cf (%i13) cf : 1.687811; (cf) 1.687811 (%i14) excel_round(float(TAS(29.92,1194.79,465.24)),1); (%o14) 118.5The conversion factor for meters per second (m/s) is: 3.28084(%i15) cf:3.28084; (cf) 3.28084 (%i16) excel_round(float(TAS(29.92,1194.79,465.24)),1); (%o16) 61.0This is TAS in meters per second. Thus, by changing the value of only the conversion factor, TAS could be determined for various units of measure. Summary This section provided examples of the use of Maxima expressions in an applied context. It was demonstrated that Maxima can be used effectively for numerical computations. However, unlike typical general-purpose programming languages, it was also demonstrated that Maxima is capable of providing much more capability than merely "number crunching" by way of symbolic operations. The symbolic expression capabilities of Maxima extend way beyond what was demonstrated in this section. The purpose of this section is not to provide a comprehensive survey of the symbolic expression capabilities of Maxima. Rather, it is intended to be an introduction to Maxima for the sake of familiarity with the sort of expressions that are typical for using Maxima in an applied context, with the syntax that is used for Maxima expressions, and generally the important role expressions have in interacting with the Maxima interface. Three Maxima operators that should be recognized as used frequently as evidenced by the expressions of this section are: = the equality operator (equation expressions) := the Maxima function definition operator : the assignment operator No matter the task at hand for which Maxima is being used, these three operators are likely to be essential, so having a clear understanding of the nature of each is important. Equations and Functions Here is what a wxMaxima "Tip of the Day" remarks about the first two:"Equations have several advantages over functions. For example they can be manipulated with factor(), expand() and similar functions. They can easily be introduced one into another."It is essential to recognize that Maxima function expressions defined using the := operator (and likewise for the define() operator, not discussed here) is not synonymous with the mathematical concept of a function, and this is also the case with typical general-purpose programming languages, and not only Maxima. The Maxima function (in the programming sense) operator can be used to define a function (in the mathematical sense), but this is not generally the case for all Maxima function expressions. Maxima function expressions can be entered that cannot be used to define a mathematical function in the technical sense. To define a mathematical function using Maxima, to be precise, would require something similar to the first two of the following expressions:(%i1) y = f(x); (%o1) y=f(x) (%i2) f(x) := x^3; (%o2) f(x):=x^3 (%i3) y = f(x), x=3; (%o3) y=27 (%i4) y = f(3); (%o4) y=27 (%i5) is(equalp(y,f(3))); (%o5) unknown (%i6) is(y = f(3)); (%o6) false (%i7) is(y = f(3)), y = 27; (%o7) true (%i8) f(3); (%o8) 27The first of the above expressions, if evaluated, will have a value of true or false, depending on the values of x and y. The second of these two expressions will have a numerical value that depends only on the numerical value of x. This difference is essential to understand. For a function, in the mathematical sense, it is necessary for values of x and y as applied to the expression y = f(x), to evaluate to true, that is the values of the ordered pairs of numbers (x,y) evaluate the expression y = f(x) to true if they are instances of the function, and false otherwise. The function f(x), in the programming sense, can be defined in all sorts of ways as a Maxima expressions, including for the purpose of the definition of a mathematical function, but not necessarily so. The following expressions from those above, provides some insight into the way Maxima as a CAS differs in an essential way from typical general-purpose programming languages:(%i5) is(equalp(y,f(3))); (%o5) unknown (%i6) is(y = f(3)); (%o6) false The value of the identifier y in these expressions is itself. because of this, y as an identifier represents the essence of a variable quantity, and not as merely as a name for another expression. The first of these two expressions has the value of unknown, because it simply is not possible to say if y and f(3) are equal. The second has a value of false simply because the value of y is itself and unless the value of f(3) is y, then this expression is obviously false.(%i1) g(x) := y; (%o1) g(x):=y (%i2) y = g(y); (%o2) y=y This is trivially obvious, but also clearly logical, and if Maxima is anything, it is logical. This is as it should be. Technically the := operator is a type of assignment. The other frequently used assignment operator is :, and a brief discussion of this operator follows. Assignment The : assignment operator is most typical used to name another value or expression. The word "identifier" in this case is understood to be a name that is in a sense identical to a value or expression. some examples follow:(%i08) x : 2; (x) 2 (%i09) x; (%o09) 2 (%i10) y : 3.14159; (y) 3.14159 (%i11) y; (%o11) 3.14159 (%i12) r : rationalize(0.1); (r) 3602879701896397/36028797018963968 (%i13) r; (%o13) 3602879701896397/36028797018963968 (%i14) s : "a string"; (s) "a string" (%i15) s; (%o15) "a string" (%i16) \9ofHearts : true; (\9ofHearts) true (%i17) \9ofHearts; (%o17) trueAssignment is not limited however to simple expressions such as these. Entire expressions from simple to more complex can be assigned to identifiers. For example:(%i18) z : m + n +3; (z) n+m+3 (%i19) z; (%o19) n+m+3 (%i20) ex1 : a*x^2 +b*x + c; (ex1) c+2*b+4*a (%i21) ex1; (%o21) c+2*b+4*a (%i22) ex2 : 'integrate( %e^(2*t+a)*sin(t)*t,t); (ex2) integrate(%e^(2*t+a)*t*sin(t),t) (%i23) ex2; (%o23) integrate(%e^(2*t+a)*t*sin(t),t) (%i24) solution : ex2, integrate; (solution) ((10*%e^(2*t+a)*t-3*%e^(2*t+a))*sin(t)+(4*%e^(2*t+a)-5*%e^(2*t+a)*t)*cos(t))/25 (%i25) solution; (%o25) ((10*%e^(2*t+a)*t-3*%e^(2*t+a))*sin(t)+(4*%e^(2*t+a)-5*%e^(2*t+a)*t)*cos(t))/25This sort of assignment is not possible using typical general-purpose programming languages, and will result in undefined variable error messages. For example in the case of Julia:julia> z = m + n + 3 ERROR: UndefVarError: `m` not defined in `Main` Suggestion: check for spelling errors or missing imports. Stacktrace: [1] top-level scope @ REPL[1]:1Julia is a very capable and useful programming language, especially for numerical computations, but symbolic expressions are not what Julia was intended to be used for. Maxima on the other hand, is quite at home with symbolic expressions. Some of the output expressions in these examples are rendered in a much more readable and conventional manner by the wxMaxima GUI using typical mathematical notation. The output for ex2 for example is rendered as: File:Image_of_wxMaxima_output.png The output for the solution identifier is rendered as: 500x500px The wxMaxima GUI can be configured to suppress the % character displayed before the special constant %e in the above output expressions if desired. The reader is referred to Sections 9.6 (Assignment operators) and 22 (Equations) of the Maxima Help document for more information related to these three operators and more.