Update vectorcalculus.mac

Added functions grad and dir_deriv to compute a gradient vector and directional derivative respectively. 

Changed all functions such that the variable list is an optional input, and in cases where it is omitted the variable list is taken from the given function. 

Added many unit tests.

stack/maxima/contrib/vectorcalculus.mac

 /*  Author Luke Longworth
     University of Canterbury
-    Copyright (C) 2023 Luke Longworth
+    Copyright (C) 2024 Luke Longworth
 
     This program is free software: you can redistribute it or modify
     it under the terms of the GNU General Public License version two.
@@ -16,16 +16,43 @@
 /****************************************************************/
 /*  Vector calculus functions for STACK                         */
 /*                                                              */
-/*  V1.0 June 2023                                              */
+/*  V2.0 March 2024                                             */
 /*                                                              */
 /****************************************************************/
 
+/* A flag used throughout the file. */
+/* If return_vect is true, then vector answers are returned as an nx1 matrix. */
+/* If return_vect is false, then vector answers are returned as a list. */
+return_vect: true;
+
+/****************************************************************/
+/* Calculate the gradient vector of a multivariate function     */
+/****************************************************************/
+grad(f, [vars]):= block([grad_vec],
+    vars: flatten(vars),
+    if emptyp(vars) then vars: listofvars(f),
+    /* TODO: confirm grad should always simplify? */
+    grad_vec: map(lambda([ex], ev(diff(f,vars[ex]), simp)), ev(makelist(ii,ii,1,length(vars)), simp)),
+    if return_vect then return(transpose(matrix(grad_vec))) else return(grad_vec)
+);
+
+s_test_case((return_vect:true, grad(x*y*z,[x,y,z])),matrix([y*z],[x*z],[x*y]));
+s_test_case((return_vect:true, grad(x*y*z,x,y,z)),matrix([y*z],[x*z],[x*y]));
+s_test_case((return_vect:true, grad(x*y*z)),matrix([y*z],[x*z],[x*y]));
+s_test_case((return_vect:false, grad(x*y*z,[x,y,z])),[y*z,x*z,x*y]);
+s_test_case((return_vect:false, grad(x*y*z,x,y,z)),[y*z,x*z,x*y]);
+s_test_case((return_vect:false, grad(x*y*z)),[y*z,x*z,x*y]);
+s_test_case((return_vect:false, grad(x^2 + x)),[2*x+1]);
+s_test_case((return_vect:true, grad(a+2*b+3*c+4*d+5*p)),matrix([1],[2],[3],[4],[5]));
+s_test_case((return_vect:true, grad(a+2*b+3*c+4*d+5*p,[p,d,c,b,a])),matrix([5],[4],[3],[2],[1]));
+
 /****************************************************************/
 /* Calculate the divergence of a vector-valued function         */
 /****************************************************************/
-div(u, vars):= block([div_vec],
-    /* TODO: error trapping: if not(listp(vars)) or emptyp(vars) then error("div: the second argument must be a list of variables."), */
+div(u, [vars]):= block([div_vec],
     if matrixp(u) then funcs: list_matrix_entries(u) else funcs: flatten(u),
+    vars: flatten(vars),
+    if emptyp(vars) then vars: listofvars(u),
     /* TODO: confirm div should always simplify? */
     div_vec: map(lambda([ex], ev(diff(funcs[ex],vars[ex]), simp)), ev(makelist(ii,ii,1,length(vars)), simp)),
     return(apply("+", div_vec))
@@ -34,15 +61,67 @@ div(u, vars):= block([div_vec],
 s_test_case(div([x^2*cos(y),y^3],[x,y]), 2*x*cos(y)+3*y^2);
 s_test_case(div(transpose(matrix([x^2*cos(y),y^3])),[x,y]), 2*x*cos(y)+3*y^2);
 s_test_case(div(matrix([x^2*cos(y),y^3]),[x,y]), 2*x*cos(y)+3*y^2);
+s_test_case(div([x^2*cos(y),y^3],[y,x]), -x^2*sin(y));
+s_test_case(div([y^3,x^2*cos(y)],[y,x]), 2*x*cos(y)+3*y^2);
+s_test_case(div([x^2*cos(y),y^3]), 2*x*cos(y)+3*y^2);
+s_test_case(div(transpose(matrix([x^2*cos(y),y^3]))), 2*x*cos(y)+3*y^2);
+s_test_case(div(matrix([x^2*cos(y),y^3])), 2*x*cos(y)+3*y^2);
+s_test_case(div([x^2*cos(y),y^3],x,y), 2*x*cos(y)+3*y^2);
+s_test_case(div(transpose(matrix([x^2*cos(y),y^3])),x,y), 2*x*cos(y)+3*y^2);
+s_test_case(div(matrix([x^2*cos(y),y^3]),x,y), 2*x*cos(y)+3*y^2);
 
 /****************************************************************/
 /* Calculate the curl of a vector-valued function               */
 /****************************************************************/
-curl(u,vars):= block([cux, cuy, cuz],
-    /* TODO: error trapping: if not(listp(vars)) or emptyp(vars) then error("curl: the second argument must be a list of 3 variables."), */
+curl(u, [vars]):= block([cux, cuy, cuz],
     if matrixp(u) then [ux,uy,uz]: list_matrix_entries(u) else [ux,uy,uz]: flatten(u),
+    vars: flatten(vars),
+    if emptyp(vars) then vars: listofvars(u),
     cux: diff(uz,vars[2]) - diff(uy,vars[3]),
     cuy: diff(ux,vars[3]) - diff(uz,vars[1]),
     cuz: diff(uy,vars[1]) - diff(ux,vars[2]),
-    return(transpose(matrix([cux,cuy,cuz])))
+    if return_vect then return(transpose(matrix([cux,cuy,cuz]))) else return([cux,cuy,cuz])
 );
+
+s_test_case((return_vect: true, curl([x*y*z,x*y*z,x*y*z],[x,y,z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl([x*y*z,x*y*z,x*y*z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: false, curl([x*y*z,x*y*z,x*y*z],[x,y,z])),[x*z-x*y,x*y-y*z,y*z-x*z]);
+s_test_case((return_vect: false, curl([x*y*z,x*y*z,x*y*z])),[x*z-x*y,x*y-y*z,y*z-x*z]);
+s_test_case((return_vect: true, curl([x*y*z,x*y*z,x*y*z],[y,z,x])),matrix([x*y-y*z],[y*z-x*z],[x*z-x*y]));
+s_test_case((return_vect: true, curl(matrix([x*y*z,x*y*z,x*y*z]),[x,y,z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z,x*y*z,x*y*z]),x,y,z)),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z,x*y*z,x*y*z]))),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z],[x*y*z],[x*y*z]),[x,y,z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z],[x*y*z],[x*y*z]),x,y,z)),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z],[x*y*z],[x*y*z]))),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+
+/*******************************************************************/
+/* Calculate the directional derivative of a multivariate function */
+/*******************************************************************/
+dir_deriv(f, u, [vars]):= block([unit_u, der],
+    if matrixp(u) then u: list_matrix_entries(u),
+    vars: flatten(vars),
+    if emptyp(vars) then vars: listofvars(f),
+    unit_u: u/sqrt(u . u),
+    der: ev(flatten(args(grad(f, vars))) . unit_u,simp),
+    return(der)
+);
+
+s_test_case((return_vect: false, dir_deriv(x*y*z,[1,2,2],[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2],[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,[1,2,2])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2],[y,z,x])),(2*y*z)/3+(x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2],x,y,z)),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,matrix([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,matrix([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]),[y,z,x])),(2*y*z)/3+(x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]),x,y,z)),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,transpose([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,transpose([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]),[y,z,x])),(2*y*z)/3+(x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]),x,y,z)),(y*z)/3+(2*x*z)/3+(2*x*y)/3);