Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 357 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Here is a complete implementation as a module. Please let me know if you have any questions. If you want to save this to a library, it may need a minor change: MyTypes might need to be global. If it doesn't work when saved to a library, let me know.

restart:

MyModule:= module()
uses TT= TypeTools;
global _T1, _T2L, _T2V, _T3L, _T3V, _MyType;
local
     MyTypes:= {_T1, _T2L, _T2V, _T3L, _T3V},
     AllMyTypes:= MyTypes union {_MyType},

     ModuleLoad:= proc()
     local
          g, #iterator over module globals
          e
     ;
          #op([2,6], ...) of a module is its globals.
          for g in op([2,6], thismodule) do
               e:= eval(g);
               if g <> e and e in AllMyTypes then
                    error "The name %1 must be globally available.", g
               end if
          end do;
          TT:-AddType(_T1, algebraic);
          TT:-AddType(_T2V, 'Vector(2, algebraic)');
          TT:-AddType(_T2L, [algebraic $ 2]);
          TT:-AddType(_T3V, 'Vector(3, algebraic)');
          TT:-AddType(_T3L, [algebraic $ 3]);
          TT:-AddType(_MyType, MyTypes)
     end proc,

     ModuleUnload:= proc()
     local T;
          for T in AllMyTypes do TT:-RemoveType(T) end do
     end proc,

     MyDispatch:= overload([
          proc(A::_T1, B, C)
          option overload;
          local r:= "A, B, C are T1."; #unnecessary; just an example.
               #statements to process this type
          end proc,

          proc(A::_T2L, B, C)
          option overload;
          local r:= "A, B, C are T2L.";
               #
          end proc,

          proc(A::_T2V, B, C)
          option overload;
          local r:= "A, B, C are T2V.";
               #
          end proc,

          proc(A::_T3L, B, C)
          option overload;
          local r:= "A, B, C are T3L.";
               #         
          end proc,

          proc(A, B, C)
          local r:= "A, B, C are T3V.";
               #
          end proc
     ]),

     ModuleApply:= proc(
          A::And(
               _MyType,
               satisfies(A-> andmap(T-> A::T implies B::T and C::T, MyTypes))
          ),
          B::_MyType, C::_MyType
     )
          MyDispatch(args)
     end proc
;
     ModuleLoad()    
end module:

#Example usage:                                          
MyModule(x,y,z);

 

The elementwise operator ~ can be used like zip, but for arbitrary arity. In this way it is more powerful than any single spelled out command in Maple such as map, zip, map2, etc. Observe:

Vals:= [x= [1,2,3], y= [4,5,6], z= [7,8,9]]:
f~(eval([x,y,z], Vals)[]);

     [f(1, 4, 7), f(2, 5, 8), f(3, 6, 9)]

I think that you'll be able on your own to apply this to your situation. However, if you need further assistance with that, just ask.

I agree with all that Tom said. In addition, if you're trying to use color to precisely encode numeric information in a plot, then having transparency and directional lighting (directional lighting is a default) makes that numeric information very difficult to read. So, in that case, I'd use this:

plots:-display(
     plottools:-cuboid([0,0,0],[1,1,1]),
     transparency= 0,
     colorscheme= ["xyzcoloring", (x,y,z)-> x^2 + y^2 - z^2],
     lightmodel= NONE
);

If your use of color is primarily aesthetic, then ignore this Answer.

Use add, not sum, because add allows a step value for the index. Also, your exponent is the inverse of what it should be. It should be (i-1)/2. So, the corrected command is

add(2*q*cos(2*i*x)/(i*Pi)*(-1)^((i-1)/2), i= 1..35, 2);

@sigl1982 It's quite refreshing to see such neatly written code posted here.

I suspect that there is some intentional randomization somewhere in the mod code that your program uses, perhaps in `mod/Factors`. Do this test with randomize: Make the first statement of your procedure randomize(1). Run the procedure several times with the same input. (Q1) Do you get identical results? Then change the first statement to randomize(2). (Q2) Do you now get different results?

  • If Q1 and Q2 are both "yes", that's very very strong evidence that intentional randomization is being used somewhere.
  • If Q1 is "no", then there's no point in asking Q2, and you have strong evidence that intentional randomization is not being used. Back to the drawing board.
  • If Q1 is "yes" and Q2 is "no", then try changing the 2 to other numbers and re-asking Q2 until the answer is "yes". If you can't make the answer "yes" by changing the argument to randomize, then the test is inconclusive.

I recommend that you turn off the ASSERTion checking while doing this simply to reduce the volume of infomation that you need to look through.

Randomized algorithms are commonly used in finite-field polynomial arithmetic. No doubt you can find several such algorithms in the book Modern Computer Algebra by von zur Gaithen and Gerhard. (Gerhard works at Maplesoft, by the way.)

 

As usual, it'd help if you post your code.

In general, the order of sets of immutable objects (such as algebraic expressions, numbers, sets, lists) is fixed, and it doesn't change between runs, even runs on different computers. However, if your sets contain mutable objects like Vectors, Matrices, or tables, then the ordering may be based on address---I'll need to research that---and of course addresses change between runs. Mutable objects also possibly include modules (which includes Objects and Records); I'll need to research that also.

However, even if your set order is based on something ephemeral like addresses, there's no way that the order is being changed by set operations such as union, minus, and indexing. The order within any given session (i.e., between restarts) is fixed.

In command-line Maple there's a command-line flag that causes Maple to use different set orders. The purpose of this is to test whether one's code depends on set order.

I think that you'll get a better efficiency gain by switching to irem on the inside:

mods(mul(irem(coeff(g[i], x, 0), p), i= S), p);

You just need to wait for the result. On my computer it took about 5 minutes.

Just a small change is needed: taylor(y(t__n+h), h, 4). This is equivalent to making h the main variable (the differentiation variable) and expanding around h=0.

The keys to doing this are using the Units package, the function-definition operator ->, and the commands solve and evalf. A benefit of using Units is that if the units come out correct, that's a fairly good check that you've set up the problem correctly.

The additional key to doing it neatly is to use eval to substitute numeric values so that that step can be put off until the end of each computation.

Glossary:

p = air density (kg/m^3)

v = velocity (m/s)

A = wing area (m^2)

C__d = drag coefficient (Units: 1)

F__d = drag force (N)

F_l = lift force (N)

C_ = lift coefficient (Units: 1)

g = acceleration due to gravity (N/kg)

m = plane's mass (kg)

 

restart:

F__d:= v-> 1/2*C__d*p*A*v^2;

proc (v) options operator, arrow; (1/2)*C__d*p*A*v^2 end proc

F__l:= v-> 1/2*C__l*p*A*v^2;

proc (v) options operator, arrow; (1/2)*C__l*p*A*v^2 end proc

with(Units:-Standard):

 

a) When an aircraft flies horizontally, the lift force F__l must balance the force of gravity m*g, where m is the aircraft’s mass and g = 9.81m/s^2 the gravitational acceleration. Compute the velocity, v__0, of a sport airplane at an altitude of h__0 = 1 000 m, where air’s density is p= 1.10 kg/m^3. The plane has mass m = 750 kg, wing area A = 18 m^2 and lift coefficient C__l = 0.35 .

Vals:= {g = 9.81*Unit(N/kg), p = 1.1*Unit(kg/m^3), m = 750*Unit(kg), A = 18*Unit(m^2), C__l = 0.35*Unit(1)};

{A = 18*Units:-Unit(('m')^2), C__l = .35, g = 9.81*Units:-Unit(('N')/('kg')), m = 750*Units:-Unit('kg'), p = 1.1*Units:-Unit(('kg')/('m')^3)}

solve({F__l(v__0) = m*g}, {v__0});

{v__0 = 2^(1/2)*(A*C__l*p*m*g)^(1/2)/(A*C__l*p)}, {v__0 = -2^(1/2)*(A*C__l*p*m*g)^(1/2)/(A*C__l*p)}

#Select positive branch; ignore negative branch.
evalf(eval(%[1], Vals));

{v__0 = 46.0801109306027*Units:-Unit(('m')/('s'))}

Vals:= Vals union %:

 

b) The propeller of an aircraft transmits power from the engine to the surrounding air. This power is used entirely to overcome the drag force. Compute the engine’s power P0 = F__d*v of the sport plane when it flies at a velocity v = v0 if the drag coefficient (based on the wing area) is C__d = 0.030.

Vals:= Vals union {C__d = 0.03*Unit(1)}:

P__0 = F__d(v__0)*v__0;

P__0 = (1/2)*C__d*p*A*v__0^3

evalf(eval(%, Vals));

P__0 = 29060.0928147354*Units:-Unit('W')

Vals:= Vals union {%}:

evalf[2](Vals);

{A = 18.*Units:-Unit(('m')^2), C__d = 0.3e-1, C__l = .35, P__0 = 0.29e5*Units:-Unit('W'), g = 9.8*Units:-Unit(('N')/('kg')), m = 0.75e3*Units:-Unit('kg'), p = 1.1*Units:-Unit(('kg')/('m')^3), v__0 = 46.*Units:-Unit(('m')/('s'))}

Something that bothers me about the above: Isn't the cross-sectional area in the lift direction different from the cross-sectional area in the drag direction?

 

 

Download airplane.mw

I've made an extensive investigation of this twice over the years. I haven't reinvestigated with Maple 2016, so this may have changed. What I found was:

  • The primary order of layering of the objects when a PLOT structure is displayed is determined by the object type (CURVES, POLYGONS, POINTS, etc.).
  • The order cannot be changed by changing the order that the objects appear in the PLOT structure.
  • The objects most in the background are those of two topological dimensions such as area fillls.
  • Next are objects of one topological dimension such as CURVES and the boundaries of POLYGONS.
  • Next are objects of zero topological dimension: POINTS.
  • Next is plotted text.
  • The object most in the foreground are the axes.
  • Within each class of objects, the order of layering can be changed by changing the order in the PLOT structure, which can in turn be controlled by the order that objects appear in a display command.

 

The plot that you want can be obtained through numeric differentiation with command fdiff:

fy:= y0-> fdiff((r-> rhs(r[-1]))@pds:-value(f(x,y), x=0), [1], [y0]):
plot(fy, 0..10, labels= [y, ``]);


You need to define an "epsilon cloud" around each integer so that if x is sufficiently close to an integer, then f(x) will be 1.

f:= x-> piecewise(abs((x-round(x))/x) < Float(6, 1-Digits), 1, 0);

The value Float(6, 1-Digits) = 6*10^(1-Digits) is often used as the target accuracy for floating-point computation.

Try any of

numtheory:-cfrac([1,1,2]);
seq(numtheory:-cfrac([1$k, 2]), k= 1..3);
Value(NumberTheory:-ContinuedFraction([1,1,2]));
seq(Value(NumberTheory:-ContinuedFraction([1$k, 2])), k= 1..3);

When using Maple's 2D Input, you can't have a space between a function name and the following left parenthesis. These extra spaces get translated as multiplication operators, which is why you see * after sphereplot and fieldplot above.

First 207 208 209 210 211 212 213 Last Page 209 of 395