Halfway symbolic halfway numeric evaluation on wxMAXIMA

Berns Buenaobra <berns.buenaobra at gmail.com> writes:
> Hello members:
>
> I have just recently adapted MAXIMA in my teaching for science education
> and physics majors in class - the symbolic math is key and the biggest
> motivator for students to have computer aided solution (to augment the
> previous class training on long hand calculus and ODE etc.).
>
> It seemed that when I started using real values on the my sheet towards the
> of end of my worksheet the earlier cells simply adapted and lost all my
> symbols and computed numerically?  Is the effect of say having assigned a
> value to a symbol say R:100, C:100e-6, L:35e-3 (from a problem in network
> analysis) global?

Stavros has suggested one sensible approach. Another approach is not to
bind R, C, L at all. Something like this:

  my_scary_formula: V = V0 * exp(-R*C*t);
  example1_vars: [R=100, C=100e-6, L=35e-3];
  example1_instance: subst(example1_vars, my_scary_formula);

Then you can fiddle with the example1_vars list and update the formulas
that depend on it. As with Stavros, I don't actually use wxMaxima, but I
guess that this method will work perfectly well there.

I find that I very rarely bind the "mathematical variables" like R,C,L
and instead end up bind names for expressions using them
(my_scary_formula etc.)

I wonder if anyone else on the list has neat tricks for the problem of
distinguishing general formulas from their specific instances?


Rupert


PS: I'm attaching an example of a calculation I did reasonably recently
    in case it gives you ideas. I've not cleaned it up at all, so if
    something looks weird... it's probably wrong! I'm actually
    calculating how big / stiff a spring I need for a personal DIY
    project, in case you're wondering. General equations are things like
    eq_moment, eq_x etc. Then it looks like I substitute in for things
    like num_theta, num_W etc. Note I have lists of numerical constants
    in the variable "measurements".

PPS: If anyone is wondering why I used but_nth and sub_from_nth rather
     than just using eliminate, so am I... I think it was because I
     wanted to control which formula got used or something. It's all a
     little hazy!

-------------- next part --------------
load(ezunits)$
load(quadratic)$

/*

  Assuming that the weight of the support structure is negligible, you
  end up with tension 

  Start with two parallel rods of length l hanging off a
  wall. Attached to them is a drill, with weight W, centre of gravity
  unspecified (since it doesn't matter. Keeping the whole thing up is
  a compression spring attached to the wall and lower rod a distance a
  from the hinge. Write F for the force from the spring. Write theta
  for the acute angle between the rod and the wall.

  Write L for the resolved force perpendicular to the rods that the
  lower hinge exerts on the drill. Resolving in this direction for the
  translation of the drill, notice that the upper rod only exerts a
  force parallel to the rods. As such, we have the simple equation:

    L = W * sin(theta)

  Now take moments of forces on the lower rod about the other
  hinge. You get

    l * L = a * F * sin(gamma)

  where gamma is the angle between the rod and the spring. The
  triangle created by the spring is isoceles and a quick check shows
  that sin(gamma) = cos(theta/2). So:
*/
eq_moment: l*W*sin(theta) = a*F*cos(theta/2);

/*
  And the cosine rule gives a calculation of x in terms of a,b and
  theta.
*/
eq_x: x = a*sqrt(2 - 2*cos(theta));

/*
  Hooke's law gives an expression for F in terms of x,s,k.
*/
eq_hooke: F = (s-x)*k;

initial_eqs: [eq_moment, eq_x, eq_hooke];

/*
  We want to try and find the correct a for a given spring, so eliminate
  psi, F, x
*/
but_nth (n, lst) :=
  block([i: 0, acc: []],
    for x in lst do if not((i: i+1) = n) then acc: cons(x, acc),
    reverse (acc))$

sub_from_nth (n, solve_for, eqns) :=
  block([sym: if atom(solve_for) then solve_for else gensym(), solns],
    solns: solve(subst(sym, solve_for, eqns[n]), sym),
    if not (length (solns) = 1) then
      error ("Complicated solution set: ", solns),
    subst (first (solns), subst(sym, solve_for, but_nth(n, eqns))))$

assume(a>0);
eqs2: sub_from_nth (2, x, initial_eqs);
eqs3: sub_from_nth (1, F, eqs2);

eq_for_a: first(eqs3);

/*
  Naively, We get two solutions because if you jam the thing really
  close to the pivot (and assume that everything stays elastic) then
  the force will be massive, but the mechanical advantage will be
  terrible, and they will cancel each other out. To ignore that
  solution, use the quadratic formula to solve the equation by hand.
*/
assume (W >= 0);
assume (k > 0);
assume (theta < %pi/2);
assume (cos(theta/2) >= 0);
assume (x > 0);
eq_a_analysis: analyse_quadratic (eq_for_a, a);

eq_a:
  (if not (assoc (solorder, eq_a_analysis) = true) then
     error ("Unknown solution order"),
   a = second(assoc (solutions, eq_a_analysis)));

/* ***************** */

/*
  Some actual numerics (I know theta, F2 etc.)

  I have a drill of length h connected below the end of the rod and
  want to have a default clearance of c above the floor. The pivot
  is a distance H above the floor.
*/
clearance_expr: H - h - cos(theta)*l = c;
eq_theta: first(solve(clearance_expr, theta));

measurements: [H = 19`cm, h = 12`cm, c = 3`cm, l=21`cm];
num_theta: rhs(expand(subst(measurements, eq_theta)));

/*
  Suppose that the drill has a mass of 750g.
*/
num_W: (0.75`kg) * (9.8`N/kg);

first_expansion (x) :=
  subst(ev(append ([theta = num_theta, W = num_W],
                   measurements),
           numer), x)$

/* ***************** */

ineq_discriminant: assoc (discriminant, eq_a_analysis) > 0;
num_ineq_discriminant: first_expansion (ineq_discriminant);

/* Inputting a metric spring */
Nmm_spring (k, rest_len, max_load, min_len) :=
  [spring, 'k = (k`N/mm``N/cm),
           's = (rest_len`mm``cm),
           'min_len = (min_len`mm``cm),
           'max_load = (max_load`N)]$

spring_calculate_a (spring) :=
  (if not (qty (subst (rest (spring), num_ineq_discriminant))) then
     error ("No solution for the spring ", spring),
   rhs (subst (rest (spring), first_expansion (eq_a))))$

spring_calculate_F (spring, a) :=
  block([Feq:
           first_expansion (
             first (solve (first (eliminate([eq_hooke, eq_x], ['x])),
                           F)))],
    subst (['a = a,
            's = assoc('s, rest (spring)),
            'k = assoc('k, rest (spring))],
           rhs (Feq)))$

spring_calculate_min_theta (spring, a) :=
  acos(subst(['x = assoc ('min_len, rest (spring)), 'a = a],
        rhs (first (solve (eq_x, cos(theta))))));

spring_downwards_movement (spring, a) :=
  first_expansion (l*(cos(spring_calculate_min_theta (spring, a)) -
                   cos(theta)))$

spring_calculate (spring) :=
  block([a: spring_calculate_a (spring), x, s, F, travel],
    x: subst ('a = a, first_expansion (rhs (eq_x))) `` cm,
    s: assoc ('s, rest (spring)) `` cm,
    min_len: assoc (min_len, rest (spring)) `` cm,
    if is (qty (x) > qty (s)) then
      error ("Spring not in compression (", x > s),
    if is (qty (x) < qty (min_len)) then
      error ("Spring too compressed (", x < min_len),
    F: spring_calculate_F (spring, a),
    travel: spring_downwards_movement (spring, a),
    ['a = a, 'x = x, 'F = F, 'travel = travel])$

/* RS!

121-220 ?8.92 for 10
(%i80) spring_calculate(Nmm_spring (4.04, 53.5, 100, 21.5))
(%o80) [a = 2.87631603480254 ` cm, x = 3.65987635551046 ` cm, 
                    F = 68.28099523737744 ` N, travel = 11.13331170542237 ` cm]

*/
1;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 315 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20130417/f9ed0af6/attachment-0001.pgp>;