restart; with(LinearAlgebra): #3D unsteady velocity field U3 in first-order Taylor approximation uuf := uu + uu_x*x + uu_y*y + uu_z*z + uu_t*t; vvf := vv + vv_x*x + vv_y*y + vv_z*z + vv_t*t; wwf := ww + ww_x*x + ww_y*y + ww_z*z + ww_t*t; 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)]]); #7D vector field P7, defined in (19) dx := u; dy := v; dz := w; du := (uuf-u)/r + ug; dv := (vvf-v)/r + vg; dw := (wwf-w)/r + wg; dt := 1; P7 := Vector([dx,dy,dz,du,dv,dw,dt]); #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]); #7D Jacobian matrix J7 of P7, (20) J7 := Matrix([ [diff(dx,x),diff(dx,y),diff(dx,z),diff(dx,u),diff(dx,v),diff(dx,w),diff(dx,t)], [diff(dy,x),diff(dy,y),diff(dy,z),diff(dy,u),diff(dy,v),diff(dy,w),diff(dy,t)], [diff(dz,x),diff(dz,y),diff(dz,z),diff(dz,u),diff(dz,v),diff(dz,w),diff(dz,t)], [diff(du,x),diff(du,y),diff(du,z),diff(du,u),diff(du,v),diff(du,w),diff(du,t)], [diff(dv,x),diff(dv,y),diff(dv,z),diff(dv,u),diff(dv,v),diff(dv,w),diff(dv,t)], [diff(dw,x),diff(dw,y),diff(dw,z),diff(dw,u),diff(dw,v),diff(dw,w),diff(dw,t)], [diff(dt,x),diff(dt,y),diff(dt,z),diff(dt,u),diff(dt,v),diff(dt,w),diff(dt,t)] ]); #6D Jacobian matrix J6 of U6 (9) 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; t := 0; #compare (20) in the paper J7; #the field in (7) in the paper F3 := Vector([ufff,vfff,wfff]); F6 := Vector([ufff,vfff,wfff,0,0,0]); #compare (11) in the paper: W := (U3-V3)/r + G3; #compare (25) in the paper: Multiply( J6, U6-F6 ); #search Eigenvector E7 of J7 belonging to Eigenvalue 0: E7 := Vector([ufff,vfff,wfff,0,0,0,1]); sol1 := solve({ Multiply(J7,E7)[1]=0, Multiply(J7,E7)[2]=0, Multiply(J7,E7)[3]=0, Multiply(J7,E7)[4]=0, Multiply(J7,E7)[5]=0, Multiply(J7,E7)[6]=0, Multiply(J7,E7)[7]=0},{ufff,vfff,wfff}); assign(sol1); #This is exactly the vector in (7): F3 := Vector([ufff,vfff,wfff]); #compare (9) in the paper: J6; #compare (11) in the paper: W := (U3-V3)/r + G3; #compare (12) in the paper: Multiply(J6,U6);