##--------------------------------------- # In this Maple sheet, we demonstrate the linear dependence of the components in vector p (Eq. 21) on the particle velocity v, see Section 4.5. ##--------------------------------------- # Restart maple restart; # Include dependent libraries with(plots): with(LinearAlgebra): # Initialize the random number generator, with a free rolls. roll := rand(1..20): roll(); roll(); roll(); roll(); roll(); roll(); ##--------------------------------------- # Rather than demonstrating the linear dependence symbolically, we show it for a random vector field. # The symbolic expressions quickly become very complicated, since a linear 18x18 system must be solved along the way. ##--------------------------------------- # Create the u component and all its first and second-order derivatives randomly. uc := (roll()-10)/10: u_x := (roll()-10)/10: u_y := (roll()-10)/10: u_z := (roll()-10)/10: u_t := (roll()-10)/10: u_xx := (roll()-10)/10: u_yy := (roll()-10)/10: u_zz := (roll()-10)/10: u_tt := (roll()-10)/10: u_xy := (roll()-10)/10: u_xz := (roll()-10)/10: u_xt := (roll()-10)/10: u_yz := (roll()-10)/10: u_yt := (roll()-10)/10: u_zt := (roll()-10)/10: # Create the v component and all its first and second-order derivatives randomly. vc := (roll()-10)/10: v_x := (roll()-10)/10: v_y := (roll()-10)/10: v_z := (roll()-10)/10: v_t := (roll()-10)/10: v_xx := (roll()-10)/10: v_yy := (roll()-10)/10: v_zz := (roll()-10)/10: v_tt := (roll()-10)/10: v_xy := (roll()-10)/10: v_xz := (roll()-10)/10: v_xt := (roll()-10)/10: v_yz := (roll()-10)/10: v_yt := (roll()-10)/10: v_zt := (roll()-10)/10: # Create the w component and all its first and second-order derivatives randomly. wc := (roll()-10)/10: w_x := (roll()-10)/10: w_y := (roll()-10)/10: w_z := (roll()-10)/10: w_t := (roll()-10)/10: w_xx := (roll()-10)/10: w_yy := (roll()-10)/10: w_zz := (roll()-10)/10: w_tt := (roll()-10)/10: w_xy := (roll()-10)/10: w_xz := (roll()-10)/10: w_xt := (roll()-10)/10: w_yz := (roll()-10)/10: w_yt := (roll()-10)/10: w_zt := (roll()-10)/10: # We create the test vector field in Taylor expansion around the space-time point (0,0,0,0). t0 := 0; x0 := 0; y0 := 0; z0 := 0; # We form the Taylor expansion of our vector field (uf, vf, wf), using the derivatives that we just defined above. uf := uc + u_x*(x-x0) + u_y*(y-y0) + u_z*(z-z0) + u_t*(t-t0) + (1/2)*u_xx*(x-x0)*(x-x0) + (1/2)*u_yy*(y-y0)*(y-y0) + (1/2)*u_zz*(z-z0)*(z-z0) + (1/2)*u_tt*(t-t0)*(t-t0) + (1/1)*u_xy*(x-x0)*(y-y0) + (1/1)*u_xz*(x-x0)*(z-z0) + (1/1)*u_xt*(x-x0)*(t-t0) + (1/1)*u_yz*(y-y0)*(z-z0) + (1/1)*u_yt*(y-y0)*(t-t0) + (1/1)*u_zt*(z-z0)*(t-t0); vf := vc + v_x*(x-x0) + v_y*(y-y0) + v_z*(z-z0) + v_t*(t-t0) + (1/2)*v_xx*(x-x0)*(x-x0) + (1/2)*v_yy*(y-y0)*(y-y0) + (1/2)*v_zz*(z-z0)*(z-z0) + (1/2)*v_tt*(t-t0)*(t-t0) + (1/1)*v_xy*(x-x0)*(y-y0) + (1/1)*v_xz*(x-x0)*(z-z0) + (1/1)*v_xt*(x-x0)*(t-t0) + (1/1)*v_yz*(y-y0)*(z-z0) + (1/1)*v_yt*(y-y0)*(t-t0) + (1/1)*v_zt*(z-z0)*(t-t0); wf := wc + w_x*(x-x0) + w_y*(y-y0) + w_z*(z-z0) + w_t*(t-t0) + (1/2)*w_xx*(x-x0)*(x-x0) + (1/2)*w_yy*(y-y0)*(y-y0) + (1/2)*w_zz*(z-z0)*(z-z0) + (1/2)*w_tt*(t-t0)*(t-t0) + (1/1)*w_xy*(x-x0)*(y-y0) + (1/1)*w_xz*(x-x0)*(z-z0) + (1/1)*w_xt*(x-x0)*(t-t0) + (1/1)*w_yz*(y-y0)*(z-z0) + (1/1)*w_yt*(y-y0)*(t-t0) + (1/1)*w_zt*(z-z0)*(t-t0); # Gravity (ug,vg,wg) and particle response time r are also chosen randomly. ug := (roll()-10)/10; vg := (roll()-10)/10; wg := (roll()-10)/10; r := roll()/20; ##--------------------------------------- # Next, we compute the Jacobian and the temporal derivative of our vector field. ##--------------------------------------- # First, we put the three vector field components (uf,vf,wf) into a vector Uf. Uf := Vector([uf,vf,wf]); # The Jacobian is computed by arranging the spatial derivatives as column vectors in a matrix. Jf := Matrix([ [diff(Uf[1],x), diff(Uf[1],y), diff(Uf[1],z)], [diff(Uf[2],x), diff(Uf[2],y), diff(Uf[2],z)], [diff(Uf[3],x), diff(Uf[3],y), diff(Uf[3],z)]]); # Finally, the temporal derivative is computed. Uf_t := Vector([diff(Uf[1],t),diff(Uf[2],t),diff(Uf[3],t)]); ##--------------------------------------- # Given the vector field and its derivatives, we can not setup the matrix M, see Eq. (22). ##--------------------------------------- # First, we introduce a few abbreviations for the location X, the particle velocity V and the gravity vector G. X := Vector([x,y,z]); V := Vector([u,v,w]); G := Vector([ug,vg,wg]); # Next, we apply the operator _p of Eq. (1) to location X, flow velocity U, inertial particle velocity V and gravity G. # Define x_p XM := Matrix([ [ 0 ,-X[3], X[2]], [ X[3], 0 , -X[1]], [-X[2], X[1], 0 ]]); # Define u_p UM := Matrix([ [ 0 ,-Uf[3], Uf[2]], [ Uf[3], 0 , -Uf[1]], [-Uf[2], Uf[1], 0 ]]); # Define g_p GM := Matrix([ [ 0 ,-G[3], G[2]], [ G[3], 0 , -G[1]], [-G[2], G[1], 0 ]]); # Define v_p VM := Matrix([ [ 0 ,-V[3], V[2]], [ V[3], 0 , -V[1]], [-V[2], V[1], 0 ]]); # Given all symbols above, we can now define matrix M according to Eq. (22). # Using a small helper variable: H1 := Multiply(-Jf,XM) + UM; M := (1/r)*Matrix([ [H1[1,1],H1[1,2],H1[1,3], Jf[1,1],Jf[1,2],Jf[1,3], XM[1,1],XM[1,2],XM[1,3], -1,0,0, 0,0,0, 0,0,0 ], [H1[2,1],H1[2,2],H1[2,3], Jf[2,1],Jf[2,2],Jf[2,3], XM[2,1],XM[2,2],XM[2,3], 0,-1,0, 0,0,0, 0,0,0 ], [H1[3,1],H1[3,2],H1[3,3], Jf[3,1],Jf[3,2],Jf[3,3], XM[3,1],XM[3,2],XM[3,3], 0,0,-1, 0,0,0, 0,0,0 ]]) + Matrix([ [GM[1,1],GM[1,2],GM[1,3], 0,0,0, 2*VM[1,1],2*VM[1,2],2*VM[1,3], 0,0,0, XM[1,1],XM[1,2],XM[1,3], -1,0,0 ], [GM[2,1],GM[2,2],GM[2,3], 0,0,0, 2*VM[2,1],2*VM[2,2],2*VM[2,3], 0,0,0, XM[2,1],XM[2,2],XM[2,3], 0,-1,0 ], [GM[3,1],GM[3,2],GM[3,3], 0,0,0, 2*VM[3,1],2*VM[3,2],2*VM[3,3], 0,0,0, XM[3,1],XM[3,2],XM[3,3], 0,0,-1 ] ]); ##--------------------------------------- # Now that we have an analytic expression for matrix M, we compute the system matrix \hat{M} in Eq. (25) and the right hand side Eq. (26). ##--------------------------------------- # We sample the neighborhood U in the spatial domain n times n := 50; # The random 6d sample points all have the same particle velocity V, since we sample in space only. for i from 1 by 1 to n do xx0[i] := (roll()-10)/10: yy0[i] := (roll()-10)/10: zz0[i] := (roll()-10)/10: uu0[i] := u0: vv0[i] := v0: ww0[i] := w0: end do: # Evaluate M and the temporal derivative U_t of the vector at the sample points for i from 1 by 1 to n do M0[i] := eval(M,{x=xx0[i],y=yy0[i],z=zz0[i],t=t0, u=uu0[i], v=vv0[i],w=ww0[i] }): U0_t[i] := eval(Uf_t,{x=xx0[i],y=yy0[i],z=zz0[i],t=t0, u=uu0[i], v=vv0[i],w=ww0[i] }): end do: # Compute the system matrix M^T M and the right hand side M^T U_t by summing up all the samples (Eqs. 25 and 26), which gives an 18x18 system matrix and a 18x1 right hand side vector. MM := Matrix(18,18); YY := Vector(18); for i from 1 by 1 to n do MM := MM + Multiply(Transpose(M0[i]),M0[i]): YY := YY + (1/r)*Multiply(Transpose(M0[i]),U0_t[i]): end do: # Compute the rank of the system matrix, which gives a rank of 12. Rank(MM); # Following Section 4.4, we introduce a regularization matrix E. EE := Matrix(18,18); EE[13,13] := 1: EE[14,14] := 1: EE[15,15] := 1: EE[16,16] := 1: EE[17,17] := 1: EE[18,18] := 1: # Apply the regularization as in Eq. (30), which gives a full rank of 18. MM := MM + EE: Rank(MM); ##--------------------------------------- # Now comes the interesting part. How does the solution depend on particle velocity v? ##--------------------------------------- # In our right hand side of the linear system only the components 7-9 depend on v. YY[1]; YY[2]; YY[3]; YY[4]; YY[5]; YY[6]; YY[7]; YY[8]; YY[9]; YY[10]; YY[11]; YY[12]; YY[13]; YY[14]; YY[15]; YY[16]; YY[17]; YY[18]; # Now, we solve the linear system in Eq. (24), which produces the optimal reference frame parameters in vector p: P := Multiply(MatrixInverse(MM),YY); # In vector p, only components 10-12 depend on the particle velocity. simplify(P[1]); simplify(P[2]); simplify(P[3]); simplify(P[4]); simplify(P[5]); simplify(P[6]); simplify(P[7]); simplify(P[8]); simplify(P[9]); simplify(P[10]); simplify(P[11]); simplify(P[12]); simplify(P[13]); simplify(P[14]); simplify(P[15]); simplify(P[16]); simplify(P[17]); simplify(P[18]); # Looking at Eq. (21), we see that these components correspond to \ddot{d}. # In other words, only \ddot{d} depends on the particle velocity. # Further, we can see that \dddot{d} and \dddot{S} both become zero as expected from the regularization.