/* file sphere.mac Edwin L Woollett, 2008, woollett@charter.net http://www.csulb.edu/~woollett See tutorial notes: Maxima by Example, Ch. 6: Differential Calculus for a discussion of this method. this file sets up and uses spherical polar coordinates: batch(sphere) sets up environment and calculates laplacian(f), divergence(Bvec), and the (r,theta,phi) components of grad(f) and curl(Bvec) */ " spherical polar coordinates (r,theta,phi ) = (r,t,p) "$ " replacement rules x,y,z to r,t,p "$ s3rule : [x = r*sin(t)*cos(p), y = r*sin(t)*sin(p), z = r*cos(t) ]$ s3sub(expr) := (subst(s3rule,expr),trigsimp(%%) )$ assume(r > 0, sin(t) > 0 )$ rxyz : sqrt(x^2 + y^2 + z^2)$ " partial derivatives of r, theta, and phi wrt x, y, and z "$ drdx : (diff(rxyz,x),s3sub(%%) ); drdy : (diff(rxyz,y),s3sub(%%) ); drdz : (diff(rxyz,z),s3sub(%%) ); dtdx : (diff( acos(z/rxyz),x ),s3sub(%%) ); dtdy : (diff( acos(z/rxyz),y ),s3sub(%%) ); dtdz : (diff( acos(z/rxyz),z ),s3sub(%%) ); dpdx : (diff( atan(y/x) ,x),s3sub(%%) ); dpdy : (diff( atan(y/x) ,y),s3sub(%%) ); " tell Maxima r=r(x,y,z), t = t(x,y,z), "$ " and p = p(x,y) "$ ( gradef(r,x,drdx),gradef(r,y,drdy),gradef(r,z,drdz) )$ ( gradef(t,x,dtdx),gradef(t,y,dtdy),gradef(t,z,dtdz) )$ ( gradef(p,x,dpdx),gradef(p,y,dpdy) )$ " tell Maxima to treat scalar function f as an "$ " explicit function of (r,t,p) "$ depends(f,[r,t,p]); " calculate the Laplacian of the scalar function f(r,t,p) "$ " using the cartesian definition "$ ( diff(f,x,2) + diff(f,y,2) + diff(f,z,2),trigsimp(%%), scanmap('multthru,%%) ); grind(%)$ /* ============================ */ " prepare to find derivatives which involve vectors "$ " cartesian unit vectors "$ ( xu : [1,0,0], yu : [0,1,0], zu : [0,0,1] )$ " spherical polar coordinate unit vectors "$ ru : sin(t)*cos(p)*xu + sin(t)*sin(p)*yu + cos(t)*zu; tu : cos(t)*cos(p)*xu + cos(t)*sin(p)*yu - sin(t)*zu; pu : -sin(p)*xu + cos(p)*yu; " cartesian def. of gradient of scalar f "$ fgradient : diff(f,x)*xu + diff(f,y)*yu + diff(f,z)*zu $ " r, theta, and phi components of grad(f) "$ fgradient_r : ( ru . fgradient, trigsimp(%%) ); fgradient_t : ( tu . fgradient, trigsimp(%%) ); fgradient_p : ( pu . fgradient, trigsimp(%%) ); " divergence of bvec "$ bvec : bx*xu + by*yu + bz*zu; " three equations which relate spherical polar components"$ " of bvec to the cartesian components "$ eq1 : br = ru.bvec; eq2 : bt = tu.bvec; eq3 : bp = pu.bvec; " invert these equations "$ sol : (linsolve([eq1,eq2,eq3],[bx,by,bz]), trigsimp(%%) ); [bx,by,bz] : map('rhs,sol)$ " tell Maxima to treat spherical polar components as "$ " explicit functions of (r,t,p) "$ depends([br,bt,bp],[r,t,p]); " calculate the divergence of bvec "$ bdivergence : ( diff(bx, x) + diff(by, y) + diff(bz, z),trigsimp(%%), scanmap('multthru,%%) ); " cartesian curl(vec) definition "$ bcurl : (diff(bz,y) - diff(by,z) )*xu + (diff(bx,z) - diff(bz,x) )*yu + (diff(by,x) - diff(bx,y) )*zu$ " find spherical polar components of curl(bvec) "$ bcurl_r : (ru.bcurl,trigsimp(%%),scanmap('multthru,%%) ); bcurl_t : (tu.bcurl,trigsimp(%%),scanmap('multthru,%%) ); bcurl_p : (pu.bcurl,trigsimp(%%),scanmap('multthru,%%) ); /* (%i1) batch(sphere); batching #pc:/work2/sphere.mac (%i2) spherical polar coordinates (r,theta,phi ) = (r,t,p) (%i3) replacement rules x,y,z to r,t,p (%i4) s3rule : [x = r sin(t) cos(p), y = r sin(t) sin(p), z = r cos(t)] (%i5) s3sub(expr) := (subst(s3rule, expr), trigsimp(%%)) (%i6) assume(r > 0, sin(t) > 0) 2 2 2 (%i7) rxyz : sqrt(z + y + x ) (%i8) partial derivatives of r, theta, and phi wrt x, y, and z (%i9) drdx : (diff(rxyz, x), s3sub(%%)) (%o9) cos(p) sin(t) (%i10) drdy : (diff(rxyz, y), s3sub(%%)) (%o10) sin(p) sin(t) (%i11) drdz : (diff(rxyz, z), s3sub(%%)) (%o11) cos(t) z (%i12) dtdx : (diff(acos(----), x), s3sub(%%)) rxyz cos(p) cos(t) (%o12) ------------- r z (%i13) dtdy : (diff(acos(----), y), s3sub(%%)) rxyz sin(p) cos(t) (%o13) ------------- r z (%i14) dtdz : (diff(acos(----), z), s3sub(%%)) rxyz sin(t) (%o14) - ------ r y (%i15) dpdx : (diff(atan(-), x), s3sub(%%)) x sin(p) (%o15) - -------- r sin(t) y (%i16) dpdy : (diff(atan(-), y), s3sub(%%)) x cos(p) (%o16) -------- r sin(t) (%i17) (gradef(r, x, drdx), gradef(r, y, drdy), gradef(r, z, drdz)) (%i18) (gradef(t, x, dtdx), gradef(t, y, dtdy), gradef(t, z, dtdz)) (%i19) (gradef(p, x, dpdx), gradef(p, y, dpdy)) (%i20) tell Maxima to treat scalar function f as an (%i21) explicit function of (r,t,p) (%i22) depends(f, [r, t, p]) (%o22) [f(r, t, p)] (%i23) calculate the Laplacian of the scalar function f(r,t,p) (%i24) (diff(f, z, 2) + diff(f, y, 2) + diff(f, x, 2), trigsimp(%%), scanmap('multthru, %%)) 2 2 d f d f df --- df --- -- cos(t) 2 2 -- 2 2 dt dp dr dt d f (%o24) --------- + ---------- + ---- + --- + --- 2 2 2 r 2 2 r sin(t) r sin (t) r dr (%i25) grind(%) 'diff(f,t,1)*cos(t)/(r^2*sin(t))+'diff(f,p,2)/(r^2*sin(t)^2)+2*'diff(f,r,1)/r +'diff(f,t,2)/r^2+'diff(f,r,2)$ (%i26) prepare to find derivatives which involve vectors (%i27) cartesian unit vectors (%i28) (xu : [1, 0, 0], yu : [0, 1, 0], zu : [0, 0, 1]) (%i29) spherical polar coordinate unit vectors (%i30) ru : cos(t) zu + sin(t) sin(p) yu + sin(t) cos(p) xu (%o30) [cos(p) sin(t), sin(p) sin(t), cos(t)] (%i31) tu : - sin(t) zu + cos(t) sin(p) yu + cos(t) cos(p) xu (%o31) [cos(p) cos(t), sin(p) cos(t), - sin(t)] (%i32) pu : cos(p) yu - sin(p) xu (%o32) [- sin(p), cos(p), 0] (%i33) cartesian def. of gradient of scalar f (%i34) fgradient : diff(f, z) zu + diff(f, y) yu + diff(f, x) xu (%i35) r, theta, and phi components of grad(f) (%i36) fgradient_r : (ru . fgradient, trigsimp(%%)) df (%o36) -- dr (%i37) fgradient_t : (tu . fgradient, trigsimp(%%)) df -- dt (%o37) -- r (%i38) fgradient_p : (pu . fgradient, trigsimp(%%)) df -- dp (%o38) -------- r sin(t) (%i39) divergence of bvec (%i40) bvec : bz zu + by yu + bx xu (%o40) [bx, by, bz] (%i41) three equations which relate spherical polar components (%i42) of bvec to the cartesian components (%i43) eq1 : br = ru . bvec (%o43) br = by sin(p) sin(t) + bx cos(p) sin(t) + bz cos(t) (%i44) eq2 : bt = tu . bvec (%o44) bt = - bz sin(t) + by sin(p) cos(t) + bx cos(p) cos(t) (%i45) eq3 : bp = pu . bvec (%o45) bp = by cos(p) - bx sin(p) (%i46) invert these equations (%i47) sol : (linsolve([eq1, eq2, eq3], [bx, by, bz]), trigsimp(%%)) (%o47) [bx = br cos(p) sin(t) + bt cos(p) cos(t) - bp sin(p), by = br sin(p) sin(t) + bt sin(p) cos(t) + bp cos(p), bz = br cos(t) - bt sin(t)] (%i48) [bx, by, bz] : map('rhs, sol) (%i49) tell Maxima to treat spherical polar components as (%i50) explicit functions of (r,t,p) (%i51) depends([br, bt, bp], [r, t, p]) (%o51) [br(r, t, p), bt(r, t, p), bp(r, t, p)] (%i52) calculate the divergence of bvec (%i53) bdivergence : (diff(bz, z) + diff(by, y) + diff(bx, x), trigsimp(%%), scanmap('multthru, %%)) dbp dbt --- --- bt cos(t) dp dt 2 br dbr (%o53) --------- + -------- + --- + ---- + --- r sin(t) r sin(t) r r dr (%i54) cartesian curl(vec) definition (%i55) bcurl : (diff(by, x) - diff(bx, y)) zu + (diff(bx, z) - diff(bz, x)) yu + (diff(bz, y) - diff(by, z)) xu (%i56) find spherical polar components of curl(bvec) (%i57) bcurl_r : (ru . bcurl, trigsimp(%%), scanmap('multthru, %%)) dbt dbp --- --- bp cos(t) dp dt (%o57) --------- - -------- + --- r sin(t) r sin(t) r (%i58) bcurl_t : (tu . bcurl, trigsimp(%%), scanmap('multthru, %%)) dbr --- dp bp dbp (%o58) -------- - -- - --- r sin(t) r dr (%i59) bcurl_p : (pu . bcurl, trigsimp(%%), scanmap('multthru, %%)) dbr --- bt dt dbt (%o59) -- - --- + --- r r dr (%o60) c:/work2/sphere.mac (%i61) */