Use %%PRETTY-FNAME in more quadpack error messages
[maxima.git] / tests / wester_problems / test_calculus.mac
blob81e706acae26cc044f3afc4b457f1557b8330673
1 /* Original version of this file copyright 1999 by Michael Wester,
2  * and retrieved from http://www.math.unm.edu/~wester/demos/Calculus/problems.macsyma
3  * circa 2006-10-23.
4  *
5  * Released under the terms of the GNU General Public License, version 2,
6  * per message dated 2007-06-03 from Michael Wester to Robert Dodier
7  * (contained in the file wester-gpl-permission-message.txt).
8  *
9  * See: "A Critique of the Mathematical Abilities of CA Systems"
10  * by Michael Wester, pp 25--60 in
11  * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
12  * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
13  */
14 /* ----------[ M a c s y m a ]---------- */
15 /* ---------- Initialization ---------- */
16 showtime: all$
17 prederror: false$
18 /* ---------- Calculus ---------- */
19 /* Calculus on a non-smooth (but well defined) function => x/|x| or sign(x)
20    */
21 diff(abs(x), x);
22 /* Calculus on a piecewise defined function */
23 a(x):= if x < 0 then -x else x$
24 /* => if x < 0 then -1 else 1 */
25 diff(a(x), x);
26 a(x):= -x*unit_step(-x) + x*unit_step(x)$
27 ratsimp(diff(a(x), x));
28 remfunction(a)$
29 /* Derivative of a piecewise defined function at a point [Herbert Fischer].
30    f(x) = x^2 - 1 for x = 1 otherwise x^3.  f(1) = 0 and f'(1) = 3 */
31 f(x):= if x = 1 then x^2 - 1 else x^3$
32 f(1);
33 diff(f(x), x);
34 subst(x = 1, %);
35 remfunction(f)$
36 /* d^n/dx^n(x^n) => n! */
37 declare(n, integer)$
38 diff(x^n, x, n);
39 makefact(%);
40 remove(n, integer)$
41 /* Apply the chain rule---this is important for PDEs and many other
42    applications => y_xx (x_t)^2 + y_x x_tt */
43 derivabbrev: true$
44 /* Set up implicit functional dependencies */
45 depends(y, x, x, t);
46 diff(y, t, 2);
47 remove([y, x], dependency)$
48 derivabbrev: false$
49 /* => f(h(x)) dh/dx - f(g(x)) dg/dx */
50 'integrate(f(y), y, g(x), h(x));
51 diff(%, x);
52 /* Exact differential => d(V(P, T)) => dV/dP DP + dV/dT DT */
53 diff(V(P, T));
54 /* Implicit differentiation => dy/dx = [1 - y sin(x y)] / [1 + x sin(x y)] */
55 y = cos(x*y) + x;
56 solve(diff(%), del(y))[1] / del(x);
57 /* => 2 (x + y) g'(x^2 + y^2) */
58 eqn: 'diff(f(x, y), x) + 'diff(f(x, y), y);
59 ev(subst(f(x, y) = g(x^2 + y^2), eqn), diff);
60 load(ndiff)$
61 ev(subst(f(x, y) = g(x^2 + y^2), eqn), diff);
62 subst(f(x, y) = g(x^2 + y^2), eqn);
63 factor(ev(%, diff));
64 convert_to_de(%);
65 remvalue(eqn)$
66 /* Residue => - 9/4 */
67 residue((z^3 + 5)/((z^4 - 1)*(z + 1)), z, -1);
68 /* Differential forms */
69 init_cartan([x, y, z])$
70 /* (2 dx + dz) /\ (3 dx + dy + dz) /\ (dx + dy + 4 dz) => 8 dx /\ dy /\ dz */
71 (2*dx + dz) ~ (3*dx + dy + dz) ~ (dx + dy + 4*dz);
72 /* d(3 x^5 dy /\ dz + 5 x y^2 dz /\ dx + 8 z dx /\ dy)
73    => (15 x^4 + 10 x y + 8) dx /\ dy /\ dz */
74 factor(ext_diff(3*x^5 * dy ~ dz + 5*x*y^2 * dz ~ dx + 8*z * dx ~ dy));
75 /* => 1 - 3/8 2^(1/3) = 0.5275296 */
76 f(x):= x^4 - x + 1$
77 min_bracket(0.0, 1.0, f);
78 errcatch(apply('min_function, endcons(f, %[1])));
79 /* => [0, 1] */
80 /*[minimize(1/(x^2 + y^2 + 1)), maximize(1/(x^2 + y^2 + 1))];*/
81 /* Minimize on [-1, 1] x [-1, 1]:
82    => min(a - b - c + d, a - b + c - d, a + b - c - d, a + b + c + d) */
83 /*minimize(a + b*x + c*y + d*x*y, [x, -1, 1], [y, -1, 1]);*/
84 /* => [-1, 1] */
85 /*[minimize(x^2*y^3, [x, -1, 1], [y, -1, 1]),
86    maximize(x^2*y^3, [x, -1, 1], [y, -1, 1])];*/
87 /* Linear programming: minimize the objective function z subject to the
88    variables xi being non-negative along with an additional set of constraints.
89    See William R. Smythe, Jr. and Lynwood A. Johnson, _Introduction to Linear
90    Programming, with Applications_, Prentice Hall, Inc., 1966, p. 117:
91    minimize z = 4 x1 - x2 + 2 x3 - 2 x4 => {x1, x2, x3, x4}  = {2, 0, 2, 4}
92    with zmin = 4 */
93 lp_by_simplex(-(4*x1 - x2 + 2*x3 - 2*x4), [2*x1 + x2 + x3 + x4 <= 10,
94               x1 - 2*x2 - x3 + x4 >= 4, x1 + x2 + 3*x3 - x4 >= 4],
95               [x1, x2, x3, x4]);