A general theoretical and computational procedure for dealing with an exponential-logarithmic kinematic model for transformation stretch tensor in a multiphase phase field approach to stress- and temperature- induced martensitic transformations with N martensitic variants is developed for transformations between all possible crystal lattices. This kinematic model, where the natural logarithm of transformation stretch tensor is a linear combination of natural logarithm of the Bain tensors, yields isochoric variant-variant transformations for the entire transformation path. Such a condition is plausible and cannot be satisfied by the widely used kinematic model where the transformation stretch tensor is linear in Bain tensors. Earlier studies can handle commutative Bain tensors only. In the present treatment, the exact expressions for the first and second derivatives of the transformation stretch tensor with respect to the order parameters are obtained. Using these relations, the transformation work for austenite <=> martensite and variant <=> variant transformations is analyzed and the thermodynamic instability criteria for all homogeneous phases are explicitly expressed. The finite element procedure with an emphasis on the derivation of the tangent matrix for the phase field equations, which involves second derivatives of the transformation deformation gradients with respect to the order parameters, is developed. Change in anisotropic elastic properties during austenite-martensitic variants and variant-variant transformations is taken into account. The numerical results exhibiting twinned microstructures for cubic to orthorhombic and cubic to monoclinic-I transformations are presented.