Merge pull request #1145 from LukeLongworth/patch-10

Update vectorcalculus.mac

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);