From f504a3078167289229690df1c3d2c3cb0417bfa5 Mon Sep 17 00:00:00 2001
From: Luke Longworth <34358809+LukeLongworth@users.noreply.github.com>
Date: Thu, 21 Mar 2024 15:26:31 +1300
Subject: [PATCH] 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 | 93 +++++++++++++++++++++++--
1 file changed, 86 insertions(+), 7 deletions(-)
diff --git a/stack/maxima/contrib/vectorcalculus.mac b/stack/maxima/contrib/vectorcalculus.mac
index 9570f3b6e..be8402151 100644
--- a/stack/maxima/contrib/vectorcalculus.mac
+++ b/stack/maxima/contrib/vectorcalculus.mac
@@ -1,6 +1,6 @@
/* 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);
--
GitLab