restart; with(LinearAlgebra): #3D steady velocity field U3 in first-order Taylor approximation uuf := uu + uu_x*x + uu_y*y + uu_z*z; vvf := vv + vv_x*x + vv_y*y + vv_z*z; wwf := ww + ww_x*x + ww_y*y + ww_z*z; U3 := Vector([uuf,vvf,wwf]); #current particle velocity V3 := Vector([u,v,w]); #gravity Vector G3 := Vector([ug,vg,wg]); #3D Spatial Jacobian matrix J3: J3 := Matrix([ [diff(uuf,x),diff(uuf,y),diff(uuf,z)], [diff(vvf,x),diff(vvf,y),diff(vvf,z)], [diff(vvf,x),diff(wwf,y),diff(wwf,z)]]); #6D vector field U6 dx := u; dy := v; dz := w; du := (uuf-u)/r + ug; dv := (vvf-v)/r + vg; dw := (wwf-w)/r + wg; U6 := Vector([dx,dy,dz,du,dv,dw]); #6D Jacobian matrix J6 J6 := Matrix([ [diff(dx,x),diff(dx,y),diff(dx,z),diff(dx,u),diff(dx,v),diff(dx,w)], [diff(dy,x),diff(dy,y),diff(dy,z),diff(dy,u),diff(dy,v),diff(dy,w)], [diff(dz,x),diff(dz,y),diff(dz,z),diff(dz,u),diff(dz,v),diff(dz,w)], [diff(du,x),diff(du,y),diff(du,z),diff(du,u),diff(du,v),diff(du,w)], [diff(dv,x),diff(dv,y),diff(dv,z),diff(dv,u),diff(dv,v),diff(dv,w)], [diff(dw,x),diff(dw,y),diff(dw,z),diff(dw,u),diff(dw,v),diff(dw,w)]]); x := 0; y := 0; z := 0; #compare (9) in the paper: J6; #compare (11) in the paper: W := (U3-V3)/r + G3; #compare (12) in the paper: Multiply(J6,U6);