Linear Quaternion Differential Equations: Basic Theory

and Fundamental Results

arXiv:1510.02224v5 [math.CA] 7 Sep 2017

†
Kit Ian Kou1,2∗ Yong-Hui Xia1,3
1. School of Mathematical Sciences, Huaqiao University, 362021, Quanzhou, Fujian, China.
xiadoc@163.com; xiadoc@hqu.edu.cn (Y.H.Xia)
2.Department of Mathematics, Faculty of Science and Technology, University of Macau, Macau
kikou@umac.mo (K. I. Kou)
3. Department of Mathematics, Zhejiang Normal University, Jinhua, 321004, China

∗
Kit Ian Kou acknowledges financial support from the National Natural Science Foundation of China
under Grant (No. 11401606), University of Macau (No. MYRG2015-00058-FST and No. MYRG099(Y1-L2)-
FST13-KKI) and the Macao Science and Technology Development Fund (No. FDCT/094/2011/A and No.
FDCT/099/2012/A3).
†
Corresponding author. Email: xiadoc@163.com. Yonghui Xia was supported by the National Natural
Science Foundation of China under Grant (No. 11671176 and No. 11271333), Natural Science Foundation of
Zhejiang Province under Grant (No. Y15A010022), Marie Curie Individual Fellowship within the European
Community Framework Programme(MSCA-IF-2014-EF, ILDS - DLV-655209), the Scientific Research Funds
of Huaqiao University and China Postdoctoral Science Foundation (No. 2014M562320).

1
 Abstract
Quaternion-valued differential equations (QDEs) is a new kind of
differential equations which have many applications in physics and life
sciences. The largest difference between QDEs and ODEs is the alge-
braic structure. Due to the non-commutativity of the quaternion alge-
bra, the algebraic structure of the solutions to the QDEs is completely
different from ODEs. It is actually a right free module, not a linear
vector space.
This paper establishes a systematic frame work for the theory
of linear QDEs, which can be applied to quantum mechanics, fluid
mechanics, Frenet frame in differential geometry, kinematic modelling,
attitude dynamics, Kalman filter design and spatial rigid body dynam-
ics, etc. We prove that the algebraic structure of the solutions to the
QDEs is actually a right free module, not a linear vector space. On the
non-commutativity of the quaternion algebra, many concepts and prop-
erties for the ordinary differential equations (ODEs) can not be used.
They should be redefined accordingly. A definition of Wronskian is in-
troduced under the framework of quaternions which is different from
standard one in the ordinary differential equations. Liouville formula
for QDEs is given. Also, it is necessary to treat the eigenvalue prob-
lems with left- and right-sides, accordingly. Upon these, we studied the
solutions to the linear QDEs. Furthermore, we present two algorithms
to evaluate the fundamental matrix. Some concrete examples are given
to show the feasibility of the obtained algorithms. Finally, a conclusion
and discussion ends the paper.

Keywords: quaternion; quantum; solution; noncommutativity; eigenvalue;

differential equation; fundamental matrix
2000 Mathematics Subject Classification: 34K23; 34D30; 37C60;
37C55;39A12

1 Introduction and Motivation

The ordinary differential equations (ODEs) have been well studied. Theory of ODEs is very
systematic and rather complete (see monographs, e.g. [1, 11, 12, 14, 23, 31, 44]). Quaternion-
valued differential equations (QDEs) is a new kind of differential equations which have many
applications in physics and life sciences. Recently, there are few papers trying to study QDEs
[7, 22, 29, 41, 42]. Up till now, the theory of the QDEs is not systemic. Many questions
remains open. For example, even the structure of the solutions to linear QDEs is unknown.
What kind of vector space is the solutions to linear QDEs? What is the difference between
QDEs and ODEs? First of all, why should we study the QDEs? In the following, we will
introduce the background of the quaternion-valued differential equations?

2
1.1 Background for QDEs

Quaternions are 4-vectors whose multiplication rules are governed by a simple non-commutative
division algebra. We denote the quaternion q = (q0 , q1 , q2 , q3 )T ∈ R4 by

q = q0 + q1 i + q2 j + q3 k,

where q0 , q1 , q2 , q3 are real numbers and i, j, k satisfy the multiplication table formed by

i2 = j2 = k2 = ijk = −1, ij = −ji = k, ki = −ik = j, jk = −kj = i.

The concept was originally invented by Hamilton in 1843 that extends the complex numbers
to four-dimensional space.

Quaternions have shown advantages over real-valued vectors in physics and engineering
applications for their powerful modeling of rotation and orientation. Orientation can be defined
as a set of parameters that relates the angular position of a frame to another reference frame.
There are numerous methods for describing this relationship. Some are easier to visualize than
the others. Each has some kind of limitations. Among them, Euler angles and quaternions
are commonly used. We give an example from [34, 35] for illustration. Attitude and Heading
Sensors (AHS) from CH robotics can provide orientation information using both Euler angles
and quaternions. Compared to quaternions, Euler angles are simple and intuitive (see Figure
1,2,3 and 4). The complete rotation matrix for moving from the inertial frame to the body
frame is given by the multiplication of three matrix taking the form

RIB (φ, θ, ψ) = RVB2 (φ)RVV12 (θ)RIV1 (ψ),

where RVB2 (φ), RVV12 (θ), RIV1 (ψ) are the rotation matrix moving from vehicle-2 frame to the body
frame, vehicle-1 frame to vehicle-2 frame, inertial frame to vehicle-1 frame, respectively.

On the other hand, Euler angles are limited by a phenomenon called “Gimbal Lock”.
The cause of such phenomenon is that when the pitch angle is 90 degrees(see Figure 5). An
orientation sensor that uses Euler angles will always fail to produce reliable estimates when
the pitch angle approaches 90 degrees. This is a serious shortcoming of Euler angles and can
only be solved by switching to a different representation method. Quaternions provide an
alternative measurement technique that does not suffer from gimbal lock. Therefore, all CH
Robotics attitude sensors use quaternions so that the output is always valid even when Euler
Angles are not.

The attitude quaternion estimated by CH Robotics orientation sensors encodes rotation

from the “inertial frame” to the sensor “body frame.” The inertial frame is an Earth-fixed
coordinate frame defined so that the x-axis points north, the y-axis points east, and the z-axis
points down as shown in Figure 1. The sensor body-frame is a coordinate frame that remains
aligned with the sensor at all times. Unlike Euler angle estimation, only the body frame and
the inertial frame are needed when quaternions are used for estimation.

3
 Figure 1: Inertial Frame

Figure 2: Vehicle-1 Frame.

Let the vector q = (q0 , q1 , q2 , q3 )T be defined as the unit-vector quaternion encoding rotation
from the inertial frame to the body frame of the sensor, where T is the vector transpose
operator. The elements q1 , q2 , and q3 are the “vector part” of the quaternion, and can be
thought of as a vector about which rotation should be performed. The element q0 is the
“scalar part” that specifies the amount of rotation that should be performed about the vector
part. Specifically, if θ is the angle of rotation and the vector e = (ex , ey , ez ) is a unit vector
representing the axis of rotation, then the quaternion elements are defined as
   
q0 cos(θ/2)
 q1   ex sin(θ/2) 
q=  q2  =  ey sin(θ/2)  .
  

q3 ez sin(θ/2)

In fact, Euler’s theorem states that given two coordinate systems, there is one invariant

4
 Figure 3: Vehicle-2 Frame.

Figure 4: Body Frame.

axis (namely, Euler axis), along which measurements are the same in both coordinate systems.
Euler’s theorem also shows that it is possible to move from one coordinate system to the other
through one rotation θ about that invariant axis. Quaternionic representation of the attitude
is based on Euler’s theorem. Given a unit vector e = (ex , ey , ez ) along the Euler axis, the
quaternion is defined to be  
cos(θ/2)
q= .
e sin(θ/2)

Besides the attitude orientation, quaternions have been widely applied to study life science,
physics and engineering. For instances, an interesting application is the description of protein
structure (see Figure 6) [33], neural networks [24], spatial rigid body transformation [15],
Frenet frames in differential geometry, fluid mechanics, quantum mechanics, and so on.

5
 Figure 5: Gimbal Lock.

Figure 6: The transformation that sends one peptide plane to the next is a screw motion.

1.1.1 Quaternion Frenet frames in differential geometry

The differential geometry of curves [16, 18, 19] traditionally begins with a vector ~x(s) that
describes the curve parametrically as a function of s that is at least three times differentiable.
Then the tangent vector T~ (s) is well-defined at every point ~x(s) and one may choose two
additional orthogonal vectors in the plane perpendicular to T~ (s) to form a complete local
orientation frame. Provided the curvature of ~x(s) vanishes nowhere, we can choose this local
coordinate system to be the Frenet frame consisting of the tangent T~ (s), the binormal B(s),~
and the principal normal N~ (s) (see Figure 7), which are given in terms of the curve itself by

6
these expressions:
~x0 (s)
T~ (s) = ,
k~x0 (s)k
~ ~x0 (s) × ~x00 (s)
B(s) = ,
k~x0 (s) × ~x00 (s)k
~ (s) = B(s)
N ~ × T~ (s),
The standard frame configuration is illustrated in Figure 7. The Frenet frames obeys the

Figure 7: Frenet frame for a curve.

following differential equations.

 ~0    ~ 
T (s) 0 κ(s) 0 T (s)
~ 0 (s)  = v(s)  −κ(s)
 B 0 ~
τ (s)   B(s)  (1.1)
~ 0 (s)
N 0 −τ (s) 0 ~ (s)
N

Here v(s) = k~x0 (s)k is the scalar magnitude of the curve derivative, and the intrinsic geometry
of the curve is embodied in the curvature κ(s) and the torsion τ (s), which may be written in
terms of the curve itself as
k~x0 (s) × ~x00 (s)k
κ(s)(s) = ,
k~x0 (s)k3
(~x0 (s) × ~x00 (s)) · ~x000 (s)
τ (s)(s) = .
k~x0 (s) × ~x00 (s)k2

Handson and Ma [30] pointed out that all 3D coordinate frame can be express in the form
of quaternion. Then q 0 (t) = (q00 (t), q10 (t), q20 (t), q30 (t)) can be written as a first-order quaternion
differential equation
1
q 0 (t) = vKq(t)
2
with  
0 −τ 0 −κ
 τ 0 κ 0 
K=  0 −κ 0
.
τ 
κ 0 −τ 0

7
1.1.2 QDEs appears in kinematic modelling and attitude dynamics

Global Positioning System (GPS) have developed rapidly in recent years. It is compulsory for
the navigation and guidance. In fact, GPS is based on attitude determination. And attitude
determination results in the measurement of vehicle orientation with respect to a reference
frame. Wertz [40] derived the general attitude motion equations for the quaternions assuming
a constant rotation rate over an infinitesimal time ∆t such that
qk+1 = qk + q̇∆t
where q = (q1 , q2 , q3 , q4 )T , and the quaternion differentiated with respect to time is approxi-
mated by
1
q̇ = Ω · qk
2
with the skew symmetric form of the body rotations about the reference frame as
 
0 ωz −ωy ωx
 −ωz 0 ωx ωy 
Ω=  ωy −ωx

0 ωz 
−ωx −ωy −ωz 0

and ω = (ωx , ωy , ωz )T is the angular velocity of body rotation.

Another attitude dynamic is the orientation tracking for humans and robots using inertial
sensors. In [32], the authors designed a Kalman filter for estimating orientation. They devel-
oped a process model of a rigid body under rotational motion. The state vector of a process
model will consists of the angular rate ω and parameters for characterizing orientation n. For
human body motions, it is shown that quaternion rates are related to angular rates through
the following first order quaternion differential equation:
1
ṅ = n ω,
2
where ω is treated as a quaternion with zero scalar part.

Moreover, a general formulation for rigid body rotational dynamics is developed using
quaternions. Udwadia et al. [39] provided a simple and direct route for obtaining Lagrange’s
equation describing rigid body rotational motion in terms of quaternions. He presented a
second order quaternionic differntial equation as follows.
4E T JE ü + 8Ė T JE u̇ + 4J0 N (u̇)u = 2E T ΓB − 2(ω T Jω)u,
where u = (u0 , u1 , u2 , u3 )T is the unit quaternion with uT u = 1 and E is the orthogonal matrix
given by  
u0 u1 u2 u3
 −u1 u0 u3 −u2 
E=  −u2 −u3 u0
,
u1 
−u3 u2 −u1 u0

8
ΓB is the four-vector containing the components of the body torque about the body-fixed
axes. ω = (0, ω1 , ω2 , ω3 )T , ω1 , ω2 , ω3 are components of angular velocity. The inertia matrix
of the rigid body is represented by Jˆ = diag(J1 , J2 , J3 )T and the 4 × 4 diagonal matrix
J = diag(J0 , J1 , J2 , J3 )T , where J0 is any positive number.

1.1.3 QDE appears in fluid mechanics

Recently Gibbon et al. [20, 21] used quaternion algebra to study the fluid mechanics. They
proved that the three-dimensional incompressible Euler equations have a natural quaternionic
Riccati structure in the dependent variable. He proved that the 4-vector ζ for the three-
dimensional incompressible Euler equations evolves according to the quaternionic Riccati equa-
tion
Dζ
+ ζ 2 + ζp = 0,
dt
where ζ = (α, χ) and ζp = (αp , χp )T are a 4-vectors, α is a scalar and χ is a 3-vector,
T

where αp , χp are defined in terms of the Hessian matrix P [20]. To convince the reader that
this structure is robust and no accident, it is shown that the more complicated equations
for ideal magneto-hydrodynamics (MHD) can also be written in an equivalent quaternionic
form. Roubtsov and Roulstone [37, 38] also expressed the nearly geostrophic flow in terms of
quaternions.

1.1.4 QDE appears in quantum mechanics

Quaternion has been used to study the quantum mechanics for a very long time, and the
quternionic quantum mechanic theory is well developed (see, for instance, Refs. [2, 3, 27, 28]
and the references therein). In the standard formulation of nonrelativistic quantum mechanics,
the complex wave function Φ(x, t), describing a particle without spin subjected to the influence
of a real potential V (x, t), satisfies the Schrödinger equation
i h ~2 2 i
∂t Φ(x, t) = ∇ − V (x, t) Φ(x, t).
~ 2m
Alder generalized the classical mechanics to quaternionic quantum mechanics [3], the anti-self-
adjoint operator
i h ~2 2 i
AV (x, t) = ∇ − V (x, t) ,
~ 2m
can be generalized by introducing the complex potential W (x, t) = |W (x, t)| exp[iθ(x, t)],

V,W i h ~2 2 i j
A (x, t) = ∇ − V (x, t) + W (x, t),
~ 2m ~
where the quaternionic imaginary units i, j and k satisfy the following associative but non-
commutative algebra:
i2 = j2 = k2 = ijk = −1.

9
As pointed out in [27], this generalization for the anti-self-adjoint Hamiltonian operator, the
quaternionic wave function Φ(x, t) satisfies the following equation:
n i h ~2 i j o
∂t Φ(x, t) = ∇2 − V (x, t) + W (x, t) Φ(x, t).
~ 2m ~
Later, [28] gave further study on the quaternionic quantum mechanics. The complex wave
function Φ(x, t), describing a particle without spin subjected to the influence of a potential
V (x, t), satisfies the Schrödinger equation

~2
i Φxx (x, t) − (iV1 + jV2 + kV3 )Φ(x, t) = ~Φt (x, t).
2m
This partial differential equation can be reduced by using the well-known separation of vari-
ables,
Φ(x, t) = φ(x) exp{−iEt/~},
to the following second order quaternion differential equation

~2
i φ̈(x) − (iV1 + jV2 + kV3 )φ(x) = −φ(x)Ei.
2m

1.2 History and motivation

Though the above mentioned works show the existence of the quaternionic differential equation
structure appearing in differential geometry, fluid mechanics, attitude dynamics, quantum
mechanics, etc. for the author knowledge, there is no research on the systematic study for
the theory of the quaternion-valued differential equations. Upon the non-commutativity of the
quaternion algebra, the study on the quaternion differential equations becomes a challenge and
much involved. The studies are rare in this area. However, if you would like to found more
behavior of a complicated dynamical system, it is necessary to make further mathematical
analysis. By using the real matrix representation of left/right acting quaternionic operators,
Leo and Ducati [29] tried to solve some simple second order quaternionic differential equations.
Later, Campos and Mawhin [7] studied the existence of periodic solutions of one-dimensional
first order periodic quaternion differential equation. Wilczynski [41] continued this study and
payed more attention on the existence of two periodic solutions of quaternion Riccati equations.
Gasull et al [22] studied a one-dimensional quaternion autonomous homogeneous differential
equation
q̇ = aq n .
Their interesting results describe the phase portraits of the homogeneous quaternion differen-
tial equation. For n = 2, 3, the existence of periodic orbits, homoclinic loops, invariant tori
and the first integrals are discussed in [22]. Zhang [42] studied the global structure of the
one-dimensional quaternion Bernoulli equations

q̇ = aq + aq n .

10
By using the Liouvillian theorem of integrability and the topological characterization of 2-
dimensional torus: orientable compact connected surface of genus one, Zhang proved that
the quaternion Bernoulli equations may have invariant tori, which possesses a full Lebesgue
measure subset of H. Moreover, if n = 2 all the invariant tori are full of periodic orbits; if
n = 3 there are infinitely many invariant tori fulfilling periodic orbits and also infinitely many
invariant ones fulfilling dense orbits.

As we know, the second order quaternionic differential equation can be easily transformed
to a two-dimensional QDEs. In fact, for a second order quaternionic differential equation

d2 x dx
+ q 1 (t) + q2 (t)x = 0,
dt2 dt
by the transformation 
x1 = x
x2 = ẋ
then we have 
ẋ1 = x1
ẋ2 = −q2 (t)x1 − q1 (t)x2 .
Most of above mentioned works focused on the one-dimensional QDE and they mainly dis-
cussed the qualitative properties of the QDE (e.g. the existence of periodic orbits, homoclinic
loops and invariant tori, integrability, the existence of periodic solutions and so on). They
did not provide the algorithm to compute the exact solutions to linear QDEs. The most
basic problem is to solve the QDEs. Since quaternions are a non-commutative extension of
the complex numbers, a lot of concepts and theory of ODEs are not valid for QDEs. Up till
now, there is no systematic theory on the algorithm of solving QDEs. What is superposition
principle? Can the solutions of QDEs span a linear vector space? How is the structure of
the general solutions? How to determine the relations of two independent solutions? How
to define the Wronskian determinant? How to construct a fundamental matrix of the linear
QDEs? There is no algorithm to compute the fundamental matrix of the linear QDEs yet. In
particular, there is no algorithm to compute the fundamental matrix by using the eigenvalue
and eigenvector theory.

In this paper, we try to establish a basic framework for the nonautonomous linear homoge-
nous two-dimensional QDEs taking form of

ẋ1 = a11 (t)x1 + a12 (t)x2
ẋ2 = a21 (t)x1 + a22 (t)x2 ,

where all the coefficients a11 , a12 , a21 , a22 are quaternionic functions. A systematic basic theory
on the solution to the higher dimensional linear QDEs is given in this paper. We estabilsh
the superposition principle and the structure of the general solutions. Moreover, we provide
two algorithms to compute the fundamental matrix of the linear QDEs. For the sake of
understanding, we focus on the 2-dimensional system. However, in fact, our basic results can
be easily extended to any arbitrary n-dimensional QDEs.

11
 In the following, we will show that there are four profound differences between QDEs and
ODEs.

1. On the non-commutativity of the quaternion algebra, the algebraic structure of the

solutions to QDEs is not a linear vector space. It is actually a right free module.
Consequently, to determine if two solutions of the QDEs are linearly independent or depen-
dent, we have to define the concept of left dependence (independence) and right dependence
(independence). While for classical ODEs, left dependence and right dependence are the same.

2. Wronskian of ODEs is defined by Caley determinant. However, since Caley determinant

depends on the expansion of i-th row and j-th column of quaternion, different expansions can
lead to different results. Wronskian of ODEs can not be extended to QDEs. It is necessary to
define the Wronskian of QDEs by a novel method (see Section 4). We use double determinants
1 1
ddetM = rddet(MM+ ) or rddet(M+ M) to define Wronskian of QDEs.
2 2
3. Liouville formula for QDEs and ODEs are different.

4. For QDEs, it is necessary to treat the eigenvalue problems with left and right, separately.
This is a large-difference between QDEs and ODEs.

1.3 Outline of the paper

In next section, quaternion algebra is introduced. In section 3, existence and uniqueness of the
solutions to QDEs are proved. Section 4 devotes to defining the Wronskian determinant for
QDEs and proving the Liouville formulas. Superposition principle and structure of the general
solution are given. In section 5, fundamental matrix and solution to linear QDEs are given.
Section 6 presents two algorithms for computing fundamental matrix. Some examples are
given to show the feasibility of the obtained algorithms. Finally, a conclusion and discussion
ends the paper.

2 Quaternion algebra

The quaternions were discovered in 1843 by Sir William Rowan Hamilton [13]. We adopt the
standard notation in [9, 13]. We denote as usual a quaternion q = (q0 , q1 , q2 , q3 )T ∈ R4 by

q = q0 + q1 i + q2 j + q3 k,

where q0 , q1 , q2 , q3 are real numbers and i, j, k symbols satisfying the multiplication table formed
by
i2 = j2 = k2 = ijk = −1.
We denote the set of quaternion by H.

12
 We can define for the quaternion q and h = h0 + h1 i + h2 j + h3 k, the inner product, and
the modulus, respectively, by
< q, h >= q0 h0 + q1 h1 + q2 h2 + q3 h3 ,
p q
|q| = < q, h > = q02 + q12 + q22 + q32 .
<q = q0 , =q = q1 i + q2 j + q3 k.
And the conjugate is
q = q0 − q1 i − q2 j − q3 k.
Thus, we have
q
qh = hq, and q −1 = , q 6= 0.
|q|2
For any q, p ∈ H, we have
|q + p| ≤ |q| + |p|, |q · p| ≤ |q| · |p|.
Let ψ : R → H be a quaternion-valued function defined on R (t ∈ R is a real variable).
We denote
  the set of such quaternion-valued functions by H ⊗ R. Denote H2 = {Ψ Ψ =
ψ1
}. Then two dimensional quaternion-valued functions defined on real domain, Ψ(t) =
 ψ2 
ψ1 (t)
∈ H2 ⊗ R, where ψi (t) = ψi0 (t) + ψi1 (t)i + ψi2 (t)j + ψi3 (t)k. And the first derivative
ψ2 (t)
of two dimensional quaternionic functions with respect to the real variable t will be denoted
by  
ψ̇1 (t) dψi dψi0 dψi1 dψi2 dψi2
Ψ̇(t) = , ψ̇i := = + i+ j+ k.
ψ̇2 (t) dt dt dt dt dt
For any A(t) = (aij (t))n×n ∈ Hn×n ⊗ R,
 
dA(t) daij (t)
= .
dt dt
For any quaternion-valued functions f (t) defined on R, f (t) = f0 (t) + f1 (t)i + f2 (t)j + f3 (t)k
on the interval [a, b] ⊂ R, if fl (t) (l = 0, 1, 2, 3) is integrable, then f (t) is integrable and
Z b Z b Z b Z b Z b
f (t)dt = f0 (t)dt + f1 (t)dti + f2 (t)dtj + f3 (t)dtk, (a ≤ t ≤ b).
a a a a a

Moreover, Z b Z b 
A(t)dt = aij (t)dt .
a a
We define the norms for the quaternionic vector and matrix. For any A = (aij )n×n ∈ Hn×n
and x = (x1 , x2 , · · · , xn )T ∈ Hn , we define
n
X n
X
kAk = |aij |, kxk = |xi |.
i,j=1 i=1

13
For A, B ∈ H2×2 and x, y ∈ H2 , we see that
kABk ≤ kAk · kBk,
kAxk ≤ kAk · kxk,
kA + Bk ≤ kAk + kBk,
kx + yk ≤ kxk + kyk.

3 Existence and Uniqueness of Solution to QDEs

Consider the linear system

    
ψ̇1 (t) a11 (t) a12 (t) ψ1 (t)
Ψ̇(t) = A(t)Ψ(t) or = , (3.1)
ψ̇2 (t) a21 (t) a22 (t) ψ2 (t)

where Ψ(t) ∈ H2 ⊗ R, A(t) ∈ H2×2 ⊗ R is 2 × 2 continuous quaternion-valued function matrix

on the interval I = [a, b], I ⊂ R.

System (3.1) is associated with the following initial value problem

Ψ(t0 ) = ξ, ξ ∈ H2 . (3.2)

To discuss the existence and uniqueness of the solutions to QDEs (3.1), we need introduce
the limit of the quaternion-valued sequence and convergence of the quaternion-valued sequence.
Let {xk } ∈ Hn be a sequence of quaternionic vector, xk = (xk1 , xk2 , · · · , xkn )T , where xki =
xki0 + xki1 i + xki2 j + xki3 k, xkil , (i = 1, 2, · · · , n; l = 0, 1, 2, 3) are real numbers. If the sequences of
real numbers {xkil } are convergent, then the sequence of quaternonic vectors {xk } is said to be
convergent. And we define

lim xki = lim xki0 + lim xki1 i + lim xki2 j + lim xki3 k.
k→∞ k→∞ k→∞ k→∞ k→∞

Let {xk (t)} ∈ Hn ⊗ R be a sequence of quaternion-valued vector functions defined on R.

xk (t) = (xk1 (t), xk2 (t), · · · , xkn (t))T , where xki (t) = xki0 (t) + xki1 (t)i + xki2 (t)j + xki3 (t)k, xkil (t) ∈
Rn ⊗ R, (i = 1, 2, · · · , n; l = 0, 1, 2, 3) are real functions. If the sequences of real-valued
functions {xkil (t)} are convergent (uniform convergent), then the sequence of quaternion-valued
vector functions {xk (t)} is said to be convergent (uniform convergent). Moreover,

lim xki (t) = lim xki0 (t) + lim xki1 (t)i + lim xki2 (t)j + lim xki3 (t)k.
k→∞ k→∞ k→∞ k→∞ k→∞

n
xk (t) is convergent (uniform convergent), then
P
If the sequence of vector functions Sn (t) =
k=1
∞
k
P
the series x (t) is said to be convergent (uniform convergent). It is not difficult to prove
k=1
that if there exist positive real constants δk such that

kxk (t)k ≤ δk , (a ≤ t ≤ b)

14
 ∞ ∞
xk (t) is absolutely
P P
and the real series δk is convergent, then the quaternion-valued series
k=1 k=1
convergent on the interval [a, b]. Certainly, it is uniform convergent on the interval [a, b].
Moreover, if the sequence of quaternion-valued vector functions {xk (t)} is uniform convergent
on the interval [a, b], then
Z b Z b
k
lim x (t)dt = lim xk (t)dt.
k→∞ a a k→∞

Based on the above concepts on the convergence of quaternion-valued series, by similar argu-
ments to Theorem 3.1 in Zhang et al [44], we have

Theorem 3.1. The initial value problem ( 3.1)-( 3.2) has exactly a unique solution.

4 Wronskian and structure of general solution to QDEs

To study the superposition principle and structure of the general solution to QDEs, we should
introduce the concept of right- module. On the non-commutativity of the quaternion algebra,
we will see that the algebraic structure of the solutions to QDEs is not a linear vector
space. It is actually a right- module. This is the largest difference between QDEs and ODEs.

Now we introduce some definitions and known results on groups, rings, modules in the
abstract algebraic theory (e.g. [25]). An abelian group is a set, A, together with an operation
· that combines any two elements a and b to form another element denoted a · b. The symbol
· is a general placeholder for a concretely given operation. The set and operation, (A, ·), must
satisfy five requirements known as the abelian group axioms:

(i) Closure: For all a, b in A, the result of the operation a · b is also in A.

(ii) Associativity: For all a, b and c in A, the equation (a · b) · c = a · (b · c) holds.

(iii) Identity element: There exists an element e in A, such that for all elements a in A, the
equation e · a = a · e = a holds.

(iv) Inverse element: For each a in A, there exists an element b in A such that a · b = b · a = e,
where e is the identity element.

(v) Commutativity: For all a, b in A, a · b = b · a.

More compactly, an abelian group is a commutative group.

A ring is a set (R, +, ·) equipped with binary operations + and · satisfying the following
three sets of axioms, called the ring axioms:

15
(i) R is an abelian group under addition, meaning that
• (a + b) + c = a + (b + c) for all a, b, c in R (+ is associative).
• a + b = b + a for all a, b in R (+ is commutative).
• There is an element 0 in R such that a+0 = a for all a in R (0 is the additive identity).
• For each a in R there exists .a in R such that a + (−a) = 0 (−a is the additive inverse
of a).

(ii) R is a monoid under multiplication, meaning that:

• (a · b) · c = a · (b · c) for all a, b, c in R (· is associative).
• There is an element 1 in R such that a · 1 = a and 1 · a = a for all a in R (1 is the
multiplicative identity).

(iii) Multiplication is distributive with respect to addition:

• a · (b + c) = (a · b) + (a · c) for all a, b, c in R (left distributivity).
• (b + c) · a = (b · a) + (c · a) for all a, b, c in R (right distributivity).

(F, +, ·) is said to be a field if (F, +, ·) is a ring and F ∗ = F/{0} with respect to mul-
tiplication · is a commutative group. The most commonly used fields are the field of real
numbers, the field of complex numbers. A ring in which division is possible but commuta-
tivity is not assumed (such as the quaternions) is called a division ring. In fact, a field is a
nonzero commutative division ring.
Remark 4.1. Different from real or complex numbers, quaternion is a division ring, not a
field due to the non-commutativity. Thus, the vector space H2 over quaternion is not a linear
vector space. We need introduce the concept of right- module.

Module over a ring R is a generalization of abelian group. You may view an R-mod as
a “vector space over R”. A right R-module is an additive abelian group A together with a
function A × R → A (by (a, r) → a · r ) such that for all a, b ∈ A and r, s ∈ R:

(i) (a + b) · r = a · r + b · r.

(ii) a · (r + s) = a · r + a · s.

(iii) (a · r) · s = a · (rs).

Similarly, we can define the left R-module.

The multiplication xr, where x ∈ A and r ∈ R, is called the right action of R on the
abelian group A. We will write AR R to indicate that A is a right R-module (over ring R).
Similarly, we can define the left R-module. And we denote the left R-module by ALR (over
ring R).

16
 Right Submodules. Let AR R R
R be a right R-module, a subset N ⊂ AR is called a right
R-submodule, if the following conditions are satisfied:

(i) NR is a subgroup of the (additive, abelian) group AR .

(ii) xr is in NR for all r ∈ R and x ∈ NR .

Similarly, we can define the concept left submodules.

Direct sums. If N1 , · · · , Nm are submodules of a module AR over a ring R, their sum

N1 + · · · + Nm is defined to be the set of all sums of elements of the Ni , namely,

N1 + · · · + Nm = {x1 + · · · + xm |xi ∈ Ni f or each i}

this is a submodule of A. We say N1 + · · · + Nm is a direct sum if this representation is unique

for each x ∈ N1 + · · · + Nm ; that is, given two representations

x1 + · · · + xm = x = y1 + · · · + ym

where xi ∈ Ni and yi ∈ Ni for each i, then necessarily xi = yi for each i.

Lemma 4.2. The following conditions are equivalent for a set N1 , · · · , Nm of submodules of
a module A:

(i) The sum N1 + · · · + Nm is direct.

P
(ii) ( Ni ) ∩ Nk = 0 for each k.
i6=k

(iii) N1 + · · · + Nk−1 ∩ Nk = 0 for each k ≥ 2.

(iv) If x1 + · · · + xm = 0 where xi ∈ Ni for each i, then each xi = 0

If N1 + · · · + Nm is direct sum, we denote it by N1 ⊕ · · · ⊕ Nm .

Given elements x1 , · · · , xk in a right- module AR R over ring R, a sum x1 r1 + · · · + xk rk ,

ri ∈ R, is called right- linear combination of the xi , with coefficients ri . The sum of all the
right- submodules xi R, (i = 1, 2, · · · , k) is a right- submodule

x1 R + · · · + xk R = {x1 r1 + · · · + xk rk |ri ∈ R}

consisting of all such right linear combinations. If AR R = x1 R + · · · + xk R, we say that

{x1 , · · · , xk } is a generating set for in a right- module AR
R . Similarly, the sum of all the left
submodules Rxi , (i = 1, 2, · · · , k) is a left submodule

Rx1 + · · · + Rxk = {r1 x1 + · · · + rk xk |ri ∈ R}

17
consisting of all such right linear combinations. If ALR = Rx1 + · · · + Rxk , we say that
{x1 , · · · , xk } is a generating set for in a left module ALR .

The elements x1 , · · · , xk ∈ AR
R are called independent if

x1 r1 + · · · + xk rk = 0, ri ∈ R implies that r1 = · · · = rk = 0.
For the sake of convenience, we call it right independence. A subset {x1 , · · · , xk } is called a
basis of a right module AR R
R if it is right independent and generates AR . Similarly, the elements
L
x1 , · · · , xk ∈ AR are called left independent if
r1 x1 + · · · + rk xk = 0, ri ∈ R implies that r1 = · · · = rk = 0.
A subset {x1 , · · · , xk } is called a basis of a left module ALR if it is left independent and generates
ALR . A module that has a finite basis is called a free module.

From the above discussion, we know that the quaternionic vector space Hn over the division
ring H is a right- or left- module. Thus, the quaternionic vectors x1 , x2 , · · · , xn , xi ∈ Hn are
right independent if
x1 r1 + · · · + xn rn = 0, ri ∈ H implies that r1 = · · · = rk = 0.
On the contrary, if there is a nonzero ri , ı = 1, 2, · · · , n such that x1 r1 + · · · + xn rn = 0,
then x1 , x2 , · · · , xn are rihgt dependent. Similarly, we can define the left independence or
dependence for the quaternionic vector Hn over the division ring H.

For the sake of easily understanding, we focus on the two dimensional system (3.1).
Definition 4.3. For two quaternion-valued functions x1 (t) and x2 (t) defined on the real in-
terval I, if there are two quaternionic constants q1 , q2 ∈ H (not both zero) such that
x1 (t)q1 + x2 (t)q2 = 0, f or any t ∈ I,
then x1 (t) and x2 (t) are said to be right linearly dependent on I. On the contrary, x1 (t) and
x2 (t) are said to be right linearly independent on I if the algebric equation
x1 (t)q1 + x2 (t)q2 = 0, f or any t ∈ I,
can only be satisfied by
q1 = q2 = 0.
If there are two quaternionic constants q1 , q2 ∈ H (not both zero) such that
q1 x1 (t) + q2 x2 (t) = 0, f or any t ∈ I,
then x1 (t) and x2 (t) are said to be left linearly dependent on I. x1 (t) and x2 (t) are said to be
left linearly independent on I if the algebric equation
q1 x1 (t) + q2 x2 (t) = 0, f or any t ∈ I,
can only be satisfied by q1 = q2 = 0.

18
 For sake of convenience, in this paper, we adopt Cayley determinant over the quaternion
[8]. For a 2 × 2 determinant we could use a11 a22 − a12 a21 (expanding along the first row). That
is,  
a11 a12
rdet = a11 a22 − a12 a21 .
a21 a22
Certainly, one can also use a11 a22 − a21 a12 (expanding along the first column) as follows.
 
a11 a12
cdet = a11 a22 − a21 a12 .
a21 a22

There are other definitions (expanding along the second column) or (expanding along the
second row). Without loss of generality, we employ rdet to proceed our study in this paper.
In fact, one can prove the main results in this paper by replacing the Cayley determinant by
other definitions.

To study the solution of Eq. (3.1), firstly, we should define the concept of Wronskian for
QDEs. Consider any two solutions of (3.1),
   
x11 (t) x12 (t)
x1 (t) = , x2 (t) = ∈ H2 ⊗ R.
x21 (t) x22 (t)

If we define Wronskian as standard consideration of ODEs,

 
x11 (t) x12 (t)
WODE (t) = rdet = x11 (t)x22 (t) − x12 (t)x21 (t). (4.1)
x21 (t) x22 (t)

This Wronskian definition cannot be extended to QDEs. In fact, let us consider two right
linearly dependent solutions of (3.1),
   
x11 (t) x12 (t)
x1 (t) = x2 (t)η, or = η, x1 (t), x2 (t) ∈ H2 ⊗ R, η ∈ H, (4.2)
x21 (t) x22 (t)

which implies
x12 (t) = x11 (t)η −1 , and x21 (t) = x22 (t)η.
Substituting these inequalities into (4.1), we have

WODE (t) = x11 (t)x22 (t) − x12 (t)x21 (t) = x11 (t)x22 (t) − x11 (t)η −1 x22 (t)η 6= 0.

The reason lies in η is a quaternion, which can not communicate. Therefore, we should change
the standard definition of Wronskian. Now we give the concept of Wronskian for QDEs as
follows.

Definition 4.4. Let x1 (t), x2 (t) be two solutions of Eq. ( 3.1). Denote
 
x11 (t) x12 (t)
M (t) = .
x21 (t) x22 (t)

19
The Wronskian of QDEs is defined by
1
WQDE (t) = ddetM(t) := rdet M(t)M+ (t)


(4.3)
2
where ddetM := rdet MM+ is called double determinant, M + is the conjugate transpose of

M , that is,  
+ x11 (t) x21 (t)
M (t) = .
x12 (t) x22 (t)
Remark 4.5. From Definition 4.4, we have

WQDE (t) = {|x11 (t)|2 |x22 (t)|2 + |x12 (t)|2 |x21 (t)|2
(4.4)
−x12 (t)x22 (t)x21 (t)x11 (t) − x11 (t)x21 (t)x22 (t)x12 (t)}.

We see that
x12 (t)x22 (t)x21 (t)x11 (t) = x11 (t)x21 (t)x22 (t)x12 (t),
which, combining with ( 4.4), implies that WQDE (t) is a real number.

For quaternionic matrix, it should be pointed out that

rdet M (t)M + (t) 6= rdetM (t) · rdetM + (t).

Proposition 4.6. If x1 (t) and x2 (t) are right linearly dependent on I then WQDE (t) = 0.

PROOF. If x1 (t) and x2 (t) are right linearly dependent on I then Eq. (4.2) holds. Then
by (4.4),

WQDE (t) = {|x12 (t)|2 |η|2 |x22 (t)|2 + |x11 (t)|2 |η|2 |x22 (t)|2
−x12 (t)x22 (t)x22 (t)|η|2 x12 (t) − x12 (t)|η|2 x22 (t)x22 (t)x12 (t)} = 0.

That is, WQDE (t) = 0.

Theorem 4.7. (Liouville formula) The Wronskian WQDE (t) of Eq. ( 3.1) satisfies the follow-
ing quaternionic Liouville formula.
Z t 
+
WQDE (t) = exp tr[A(s) + A (s)]ds WQDE (t0 ), (4.5)
t0

where trA(t) is the trace of the coefficient matrix A(t), i.e. trA(t) = a11 (t) + a22 (t). Moreover,
if WQDE (t) = 0 at some t0 in I then WQDE (t) = 0 on I.

20
 PROOF. Let us consider Eq. (4.8). By calculating the first derivative of the left-hand-side
and right-hand-side terms, by using (4.4), we obtain

d
WQDE (t)
dt
d
rdet M (t)M + (t)


=
dt h
d
= x11 (t)x11 (t)x22 (t)x22 (t) + x12 (t)x12 (t)x21 (t)x21 (t)
dt i
−x12 (t)x22 (t)x21 (t)x11 (t) − x11 (t)x21 (t)x22 (t)x12 (t)
d d
= x11 (t)x11 (t)x22 (t)x22 (t) + x11 (t) x11 (t)x22 (t)x22 (t)
dt dt
d d
+x11 (t)x11 (t) x22 (t)x22 (t) + x11 (t)x11 (t)x22 (t) x22 (t)
dt dt
d d
+ x12 (t)x12 (t)x21 (t)x21 (t) + x12 (t) x12 (t)x21 (t)x21 (t)
dt dt
d d
+x12 (t)x12 (t) x21 (t)x21 (t) + x12 (t)x12 (t)x21 (t) x21 (t)
dt dt
d d
− x12 (t)x22 (t)x21 (t)x11 (t) − x12 (t) x22 (t)x21 (t)x11 (t)
dt dt
d d
−x12 (t)x22 (t) x21 (t)x11 (t) − x12 (t)x22 (t)x21 (t) x11 (t)
dt dt
d d
− x11 (t)x21 (t)x22 (t)x12 (t) − x11 (t) x21 (t)x22 (t)x12 (t)
dt dt
d d
−x11 (t)x21 (t) x22 (t)x12 (t) − x11 (t)x21 (t)x22 (t) x12 (t)
dt dt
= [a11 (t)x11 (t) + a12 (t)x21 (t)]x11 (t)x22 (t)x22 (t)
+x11 (t)[x11 (t)a11 (t) + x21 (t)a12 (t)]x22 (t)x22 (t)
+x11 (t)x11 (t)[a21 (t)x12 (t) + a22 (t)x22 (t)]x22 (t)
+x11 (t)x11 (t)x22 (t)[x12 (t)a21 (t) + x22 (t)a22 (t)]
+[a11 (t)x12 (t) + a12 (t)x22 (t)]x12 (t)x21 (t)x21 (t)
+x12 (t)[x12 (t)a11 (t) + x22 (t)a12 (t)]x21 (t)x21 (t)
+x12 (t)x12 (t)[a21 (t)x11 (t) + a22 (t)x21 (t)]x21 (t)
+x12 (t)x12 (t)x21 (t)[x11 (t)a21 (t) + x21 (t)a22 (t)]
−[a11 (t)x12 (t) + a12 (t)x22 (t)]x22 (t)x21 (t)x11 (t)
−x12 (t)[x12 (t)a21 (t) + x22 (t)a22 (t)]x21 (t)x11 (t)
−x12 (t)x22 (t)[a21 (t)x11 (t) + a22 (t)x21 (t)]x11 (t)
−x12 (t)x22 (t)x21 (t)[x11 (t)a11 (t) + x21 (t)a12 (t)]
−[a11 (t)x11 (t) + a12 (t)x21 (t)]x21 (t)x22 (t)x12 (t)
−x11 (t)[x11 (t)a21 (t) + x21 (t)a22 (t)]x22 (t)x12 (t)
−x11 (t)x21 (t)[a21 (t)x12 (t) + a22 (t)x22 (t)]x12 (t)
−x11 (t)x21 (t)x22 (t)[x12 (t)a11 (t) + x22 (t)a12 (t)]

21
 n
= a11 (t)|x11 (t)|2 |x22 (t)|2 + a12 (t)x21 (t)x11 (t)|x22 (t)|2
+|x11 (t)|2 a11 (t)|x22 (t)|2 + x11 (t)x21 (t)a12 (t)]x22 (t)x22 (t)
+|x11 (t)|2 a21 (t)x12 (t)x22 (t) + |x11 (t)|2 a22 (t)|x22 (t)|2
+|x11 (t)|2 x22 (t)x12 (t)a21 (t) + |x11 (t)|2 |x22 (t)|2 a22 (t)]
+a11 (t)|x12 (t)|2 |x21 (t)|2 + a12 (t)x22 (t)x12 (t)|x21 (t)|2
+|x12 (t)|2 a11 (t)|x21 (t)|2 + x12 (t)x22 (t)a12 (t)|x21 (t)|2
+|x12 (t)|2 a21 (t)x11 (t)x21 (t) + |x12 (t)|2 a22 (t)|x21 (t)|2
+|x12 (t)|2 x21 (t)x11 (t)a21 (t) + |x12 (t)|2 |x21 (t)|2 a22 (t)]
−a11 (t)x12 (t)x22 (t)x21 (t)x11 (t) − a12 (t)|x22 (t)|2 x21 (t)x11 (t)
−|x12 (t)|2 a21 (t)x21 (t)x11 (t) − x12 (t)x22 (t)a22 (t)]x21 (t)x11 (t)
−x12 (t)x22 (t)a21 (t)|x11 (t)|2 − x12 (t)x22 (t)a22 (t)x21 (t)]x11 (t)
−x12 (t)x22 (t)x21 (t)x11 (t)a11 (t) − x12 (t)x22 (t)|x21 (t)|2 a12 (t)]
−a11 (t)x11 (t)x21 (t)x22 (t)x12 (t) − a12 (t)|x12 (t)|2 x22 (t)x12 (t)
−|x11 (t)|2 a21 (t)x22 (t)x12 (t) − x11 (t)x21 (t)a22 (t)]x22 (t)x12 (t)
−x11 (t)x21 (t)a21 (t)|x12 (t)|2 − x11 (t)x21 (t)a22 (t)x22 (t)]x12 (t)o
−x11 (t)x21 (t)x22 (t)x12 (t)a11 (t) − x11 (t)x21 (t)|x22 (t)|2 a12 (t)
n
a11 (t) + a22 (t) + a11 (t) + a22 (t) |x11 (t)|2 |x22 (t)|2

=
 
+ a11 (t) + a22 (t) + a11 (t) + a22 (t) |x12 (t)|2 |x21 (t)|2 
−a 11 (t) x 12 (t)x 22 (t)x 21 (t)x 11 (t) + x 11 (t)x 21 (t)x 22 (t)x
 12 (t)
− x12 (t)x22 (t)x  21 (t)x11 (t) + x11 (t)x21 (t)x22 (t)x12 (t) a11 (t)
−x12 (t)x22 (t)a22 (t) + a22 (t)x21 (t)x11 (t)
(4.6)
−x11 (t)x21 (t) a22 (t) + a22 (t) x22 (t)x12 (t)
+|x11 (t)|2 a21 (t)x12 (t)x22 (t) + |x11 (t)|2 x22 (t)x12 (t)a21 (t)
−|x11 (t)|2 a21 (t)x22 (t)x12 (t) − a21 (t)x12 (t)x22 (t)|x11 (t)|2
+|x12 (t)|2 a21 (t)x11 (t)x21 (t) + |x12 (t)|2 x21 (t)x11 (t)a21 (t) o
−|x12 (t)|2 a21 (t)x21 (t)x11 (t) − a21 (t)x11 (t)x21 (t)|x12 (t)|2 .
Set
Q1 = x12 (t)x22 (t)x21 (t)x11 (t) + x11 (t)x21 (t)x22 (t)x12 (t)
It is easy to see that
Q1 = Q1 ,
which implies that Q1 is a real number. Thus,
 
−a
 11 (t) x 12 (t)x 22 (t)x 21 (t)x 11 (t) + x 11 (t)x 21 (t)x 22 (t)x
 12 (t)
−x12 (t)x22 (t)x21 (t)x11 (t) + x11 (t)x21 (t)x22 (t)x12 (t) a11 (t) 
= − a11 (t) + a11 (t) x12 (t)x22 (t)x21 (t)x11 (t) + x11 (t)x21 (t)x22 (t)x12 (t) .
It is obvious that
|x11 (t)|2 a21 (t)x12 (t)x22 (t) + |x11 (t)|2 x22 (t)x12 (t)a21 (t) = 2|x11 (t)|2 <{a21 (t)x12 (t)x22 (t)},
and
|x11 (t)|2 a21 (t)x22 (t)x12 (t) + x12 (t)x22 (t)a21 (t)|x11 (t)|2 = 2|x11 (t)|2 <{a21 (t)x22 (t)x12 (t)}.

22
Now we claim that
<{a21 (t)x12 (t)x22 (t)} = <{a21 (t)x12 (t)x22 (t)}. (4.7)
In fact, for any a, b ∈ H, it is easy to check that

<{ab} = <{ab}.

Thus, it follows from (4.7) that

|x11 (t)|2 a21 (t)x12 (t)x22 (t) + |x11 (t)|2 x22 (t)x12 (t)a21 (t)
−|x11 (t)|2 a21 (t)x22 (t)x12 (t) − x12 (t)x22 (t)a21 (t)|x11 (t)|2 = 0

Similarly, we have

|x12 (t)|2 a21 (t)x11 (t)x21 (t) + |x12 (t)|2 x21 (t)x11 (t)a21 (t)
−|x12 (t)|2 a21 (t)x21 (t)x11 (t) − a21 (t)x11 (t)x21 (t)|x12 (t)|2 = 0.

Substituting above equalities into (4.6), we have

d n
a11 (t) + a22 (t) + a11 (t) + a22 (t) |x11 (t)|2 |x22 (t)|2

WQDE (t) =
dt  
+a11 (t) + a22 (t) + a11 (t) + a22 (t)|x12 (t)|2 |x21 (t)|2
− a11 (t) + a22 (t) + a11 (t) + a22 (t) x12 (t)x22 (t)x21 (t)x11 (t) o (4.8)
  
− a11 (t) + a22 (t) + a11 (t) + a22 (t) x11 (t)x21 (t)x22 (t)x12 (t)
 
= a11 (t) + a22 (t) + a11 (t) + a22 (t) WQDE (t)
= trA(t) + trA+ (t) WQDE (t).

Integrating (4.8) over [t0 , t], Liouville formula (4.5) follows. Consequently, if WQDE (t) = 0 at
some t0 in I then WQDE (t) = 0 on I.

Remark 4.8. It should be noted that the Wronskian of QDEs can also be defined by

WQDE (t) = ddetM(t) := rdet M+ (t)M(t) .

By this definition, Proposition 4.6 and Theorem 4.7 can be similarly proved.

We need a lemma from (Theorem 2.10 [26]).

Lemma 4.9. The necessary and sufficient condition of invertibility of quaternionic matrix M
is ddetM 6= 0 (or WQDE (t) 6= 0).

Proposition 4.10. If WQDE (t) = 0 at some t0 in I then x1 (t) and x2 (t) are right linearly
dependent on I.

PROOF. From Liouville formula (Theorem 4.7), we have

WQDE (t0 ) = 0, implies WQDE (t) = 0, for any t ∈ I,

23
which implies that the quaternionic matrix M is not invertible on I, in view of Lemma 4.1.
Hence, the linear system

M q = 0, or (x1 (t), x2 (t))q = 0, q = (q1 , q2 )T ∈ H,

has a non-zero solution. Consequently, the two solution x1 (t) and x2 (t) are right linearly
dependent on I.

From Proposition 4.6, Proposition 4.10 and Theorem 4.7, we immediately have

Theorem 4.11. Let A(t) in Eq. ( 3.1) be continuous functions of t on an interval I. Two
solutions x1 (t), x2 (t) of Eq. ( 3.1) on I are right linearly dependent on I if and only if the
absolute value of the Wronskian, WQDE (t), is zero at some t0 in I.

Now, we present an important result on the superposition principle and structure of the
general solution.

Theorem 4.12. (Structure of the general solution) If x1 (t) and x2 (t) are two solutions of Eq.
( 3.1), then
x(t) = x1 (t)q1 + x2 (t)q2 ,
is also a solution of ( 3.1), where q1 , q2 are quaternionic constants. Let SH denote the set
which contains of all the solutions of ( 3.1). Then the set SH is a right free module (the
quaternionic vector space H2 over a division ring H).

Remark 4.13. We consider the one-dimensional nonautonomous system

ẋ = a(t)x.

We see that it has a fundamental solution denoted by x̄(t). The set of solutions S = {x̄(t)q|q ∈
H} is a 1-dimensional right free module. In fact, qx̄(t) does not satisfy the equation, that is
qx̄(t) 6∈ S. Simple computation shows:

d
(qx̄(t)) = q¯˙x(t) = qa(t)x̄(t) 6= a(t)qx̄(t).
dt
Consequently, qx̄(t) is not a solution of ẋ = a(t)x, i.e. qx̄(t) 6∈ S.

In particular, we consider the one-dimensional autonomous system

ẋ = ax.

We see that it has a fundamental solution denoted by x̄(t) = exp{at}. The set of solutions
S = {exp{at}q|q ∈ H} is a 1-dimensional right free module.

Suppose that x1 (t), x2 (t) are two fundamental solutions (they are independent in classical
sense, i.e., the quaternionic vector space H2 over a field C) of (3.1). It is easy to see that

24
the solution set SC = {x(t) is a solution of (3.1) x(t) = c1 x1 (t) + c2 x2 (t), c1 , c2 ∈ C} (the
quaternionic vector space H2 over a field C) is a linear vector space. However, we claim that
the set SC is not complete. The following counterexample shows this fact.
Couterexample: Consider the following QDEs
 
i 0
ẋ = x, x = (x1 , x2 )T . (4.9)
0 j
It is easy to see that x1 (t) = (eit , 0)T , x2 (t) = (0, ejt )T are two fundamental solutions of (4.9).
Also, it is easy to verify that x(t) = (eit k, 0)T is a solution of (4.9). Considering the solution
set SC = {x(t) is a solution of (3.1) x(t) = c1 x1 (t) + c2 x2 (t), c1 , c2 ∈ C}, we claim that the
set SC is not complete. If this is not the case, there exists the complex numbers a, b such that
x(t) = x1 (t)a + x2 (t)b, i.e., (eit k, 0)T = (eit , 0)T a + (0, ejt )T b. Then,
 it
e k = eit a,
(4.10)
0 = ejt b.
Consequently, we see that −ak = 1, b = 0. Therefore, a = k, b = 0, which contradicts to that
a is a complex number.

PROOF of THEOREM 4.12. The first assertion is easily obtained by direct differen-
tiation. Now we prove the second assertion. We proceed this with two steps. First, we will
prove the existence of two linearly independent solutions of Eq. (3.1). We can easily choose
two linearly independent quaternionic vector x01 = (x011 , x021 )T ∈ H2 and x02 = (x012 , x022 )T ∈ H2
(e.g. the natural basis e1 = (1, 0)T , e2 = (0, 1)T ). According to Theorem 3.1, for any k = 1, 2
and t0 ∈ I, Eq. (3.1) exists a unique solution xk (t), k = 1, 2 satisfying xk (t0 ) = x0k . If there
are two quaternionic constants q1 , q2 ∈ H such that
x1 (t)q1 + x2 (t)q2 = 0, f or any t ∈ I.
In particularly, taking t = t0 , we have
x1 (t0 )q1 + x2 (t0 )q2 = 0, i.e., x01 q1 + x02 q2 = 0.
Noting that x01 and x02 are two right linearly independent quaternionic vector, it follows that
q1 = q2 = 0. Thus, x1 (t), x2 (t) of Eq. (3.1) on I are right linearly dependent on I.

Secondly, we prove that any solution x(t) of Eq. (3.1) can be represented by the combina-
tion of above two right linearly independent solutions x1 (t), x2 (t)
x(t) = x1 (t)q1 + x2 (t)q2 , f or any t ∈ I, (4.11)
where q1 , q2 ∈ H. On one hand, x01 and x02 are two right linearly independent quaternionic
vector, which are a basis of two dimensional vector space H2 . Thus, there are two quaternionic
constants q1 , q2 ∈ H such that
x(t0 ) = x01 q1 + x02 q2 = x1 (t0 )q1 + x2 (t0 )q2 .
On the other hand, from the first assertion of this theorem, we know that x1 (t)q1 + x2 (t)q2 is
also a solution of (3.1) satisfying the IVP x(t0 ). Therefore, x(t) and x1 (t)q1 + x2 (t)q2 are two
solutions of Eq. (3.1) satisfying the same IVP x(t0 ). By the uniqueness theorem (Theorem
3.1), the equality (4.11) holds.

25
5 Fundamental matrix and solution to QDEs

To study the solutions of Eq. (3.1), it is necessary to introduce the quaternion-valued expo-
nential function. If q ∈ H, the exponential of q is defined by
∞
X qn
exp q =
n=0
n!

(the series converges absolutely and uniformly on compact subsets). If r ∈ H is such that
qr = rq, the exponential satisfies the addition theorem

exp{q} exp{r} = exp{q + r}.

=q
exp{q} = (exp <q)[(cos k=qk)e + (sin k=qk) ].
k=qk
Thus, exp q = 1 if and only if

<q = 0, and k=qk = 0, (mod2π).

Since [q(t)r(t)]0 = q 0 (t)r(t) + q(t)r0 (t), ( here 0 = d

dt
),
n−1
X
0
n
[q (t)] = q j (t)q 0 (t)q n−1−j (t). (5.1)
j=0

Remark 5.1. It is known that if q(t) reduces to the real-valued function,

[q n (t)]0 = nq n−1 (t), (5.2)

which is different from ( 5.1). This is the difference between quaternion-valued function and
real function (complex function). But if q is a quaternion function and q(t)q 0 (t) = q 0 (t)q(t),
then ( 5.2) holds.

Proposition 5.2. If q(t) is differentiable and if q(t)q 0 (t) = q 0 (t)q(t), it follows from that

[exp q(t)]0 = [exp q(t)]q 0 (t). (5.3)

Remark 5.3. If q(t) ≡ q is a constant quaternion, then q(t)q 0 (t) = q 0 (t)q(t) is always satisfied.
We see the difference of the derivative between the quaternionic exponential function and
complex (real) exponential function. For quaternion-value functions, the condition q(t)q 0 (t) =
q 0 (t)q(t) is needed to guarantee that the equality ( 5.3) hold. For the real or complex function,
q(t)q 0 (t) = q 0 (t)q(t) is always satisfied.

We prove Proposition 5.2 by two methods.

26
 PROOF 1: The way of Taylor expansion. Since [q(t)r(t)]0 = q 0 (t)r(t) + q(t)r0 (t), by using
(5.1),
∞ Xn−1 j
0
P∞ qn (t) 0 X q (t)q 0 (t)q n−1−j (t)
[exp q(t)] = [ n=0 n! ] =
n=0 j=0
n!
∞
X nq n−1 (t) ∞
X q n−1 (t)
= q 0 (t) = q 0 (t) = [exp q(t)]q 0 (t).
n=0
n! n=0
(n − 1)!

PROOF 2: The way of Euler formula. We rewrite q(t) as

q(t) = q0 (t) + q1 (t)i + q2 (t)j + q3 (t)k = q0 (t) + =q(t),

where
=q(t)
=q(t) = q1 (t)i + q2 (t)j + q3 (t)k = |=q(t)|.
|=q(t)|
Then we can rewrite exp{q(t)} as

eq(t) = eq0 (t)+=q(t)

= eq0 (t) e=q(t)
=q(t)
= eq0 (t) e |=q(t)| |=q(t)|
=q(t)
= eq0 (t) [cos(|=q(t)|) + |=q(t)|
sin(|=q(t)|)]

If q(t) = 0, obviously, the equality (5.3) holds. If q(t) 6= 0, differentiate exp{q(t)} with respect
to t, we have
d q(t)
e
dt n
d q0 (t) =q(t) o
= e [cos(|=q(t)|) + sin(|=q(t)|)]
dt |=q(t)|
=q(t)
= q00 (t)eq0 (t) [cos(|=q(t)|) + sin(|=q(t)|)]
|=q(t)|
d =q(t)
+eq0 (t) [cos(|=q(t)|) + sin(|=q(t)|)]
dt |=q(t)|
h =q(t) i
= q00 (t)eq0 (t) cos(|=q(t)|) + sin(|=q(t)|)
|=q(t)|
q0 (t)
h d =q 0 (t)|=q(t)| − =q(t) dt
d
|=q(t)|
+e − sin(|=q(t)|) |=q(t)| + 2
sin(|=q(t)|)
dt |=q(t)|
=q(t) d i
+ cos(|=q(t)|) |=q(t)|
|=q(t)| dt

27
 d
=q 0 (t)|=q(t)|−=q(t) dt |=q(t)|
Since q(t)q 0 (t) = q 0 (t)q(t), |=q(t)| 2 = 0, it follows that
d q(t)
e
dt h =q(t) i
= q00 (t)eq0 (t) cos(|=q(t)|) + sin(|=q(t)|)
|=q(t)|
q0 (t)
h d =q(t) d i
+e − sin(|=q(t)|) |=q(t)| + cos(|=q(t)|) |=q(t)|
dt |=q(t)| dt
0 q0 (t)
h =q(t) i
= q0 (t)e cos(|=q(t)|) + sin(|=q(t)|)
|=q(t)|
q0 (t) =q(t) d
h i
−1
+e |=q(t)| − =q (t)|=q(t)| sin(|=q(t)|) + cos(|=q(t)|) .
|=q(t)| dt
Note that
=q(t)
=q −1 (t) = ,
|=q(t)|2
and −=q(t) = =q(t). Thus, we have
d q(t)
e
dt h =q(t) i
= q00 (t)eq0 (t) cos(|=q(t)|) + sin(|=q(t)|) (5.4)
|=q(t)|
d =q(t) q0 (t) h =q(t) i
+ |=q(t)| e cos(|=q(t)|) + sin(|=q(t)|) .
dt |=q(t)| |=q(t)|
From q(t)q 0 (t) = q 0 (t)q(t), we conclude that =q(t)=q 0 (t) = =q 0 (t)=q(t). Substituting =q(t) by
q1 (t)i + q2 (t)j + q3 (t)k, we obtain
q1 (t)q20 (t) = q10 (t)q2 (t),
q1 (t)q30 (t) = q10 (t)q3 (t),
q2 (t)q30 (t) = q20 (t)q3 (t).
By these equalities, we have
d =q(t) q1 (t)q10 (t) + q2 (t)q20 (t) + q3 (t)q30 (t) q1 (t)i + q2 (t)j + q3 (t)k
|=q(t)| = p p
dt |=q(t)| q12 (t) + q22 (t) + q22 (t) q12 (t) + q22 (t) + q22 (t)
[q1 (t) + q2 (t) + q2 (t)][q1 (t)i + q2 (t)j + q30 (t)k]
2 2 2 0 0
=
q12 (t) + q22 (t) + q22 (t)
0
= =q (t).
Consequently, it follows from (5.4) that
d q(t)
e
dt
0 q0 (t)
h =q(t) i
= q0 (t)e cos(|=q(t)|) + sin(|=q(t)|)
|=q(t)|
h =q(t) i
+=q 0 (t)eq0 (t) cos(|=q(t)|) + sin(|=q(t)|)
|=q(t)|
= (q0 (t) + =q 0 (t))eq(t)
= q 0 (t)eq(t) .

28
This completes the proof of Proposition 5.2.

If a matrix A ∈ H2×2 , the exponential of A is defined by

∞
X An A A2
exp A = =E+ + + ··· .
n=0
n! 1! 2!
Moreover, for t ∈ R,
∞
X An tn At A2 t2
exp(At) = =E+ + + ··· .
n=0
n! 1! 2!

Proposition 5.4. The sires exp A = ∞ An

P
n=0 n!
is absolutely convergent. The sires exp(At) =
P ∞ An
n=0 n! is uniformly convergent on any finite interval I.

PROOF. First, we show that the sires exp A = ∞ An

P
n=0 n! is absolutely convergent. In fact,
for any positive integer k, we have
Ak kAkk
k k≤ .
k! k!
For any quaternionic matrix A, kAk is a real number. Then it is easy to see that the real
kAkn
series ∞
P P∞ A n
n=0 n! is convergent. Thus, the sires exp A = n=0 n! is absolutely convergent.

Now we prove the second conclusion. In fact, for any finite interval I, without loss of
generality, we assume that t ≤ µ. Then we have
(At)k kAkk ktkk kAkk kµkk
k k≤ ≤ .
k! k! k!
k k
Note that kAk k!kµk is convergent. Thus, the sires ∞ An
P
n=0 n! is uniformly convergent on any
finite interval I.
Proposition 5.5. For two matrix A, B ∈ H2×2 , if AB = BA, then
exp(A + B) = exp A exp B.

From Proposition 5.2, we have

Proposition 5.6. If A(t) is differentiable and if A(t)A0 (t) = A0 (t)A(t), then
[exp A(t)]0 = [exp A(t)]A0 (t).
Definition 5.7. Let x1 (t), x2 (t) be any two solutions of Eq. ( 3.1) on I. Then we call
 
T x11 (t) x12 (t)
M (t) = (x1 (t), x2 (t)) =
x21 (t) x22 (t)
as a solution matrix of Eq. ( 3.1). Moreover, if x1 (t), x2 (t) be two right linearly independent
solutions of Eq. ( 3.1) on I, the solution matrix of Eq. ( 3.1) is said to be a fundamental matrix
of Eq. ( 3.1). Further, if M (t0 ) = E(identity), M (t) is said to be a normal fundamental
matrix.

29
Remark 5.8. Note that fundamental matrix is not unique. From Definition 5.7, it is easy to
know that if M (t) is a solution matrix of Eq. ( 3.1). M (t) also satisfies Eq. ( 3.1), that is,

Ṁ (t) = A(t)M (t).

From Theorem 4.11, we know that if M (t) is a fundamental matrix of Eq. (3.1), the
Wronskian determinant WQDE (t) 6= 0. By Theorem 4.12, we have

Theorem 5.9. Let M (t) be a fundamental matrix of Eq. ( 3.1). Any solution x(t) of Eq.
( 3.1) can be represented by
x(t) = M (t)q,
where q is a constant quaternionic vector. Moreover, for given IVP x(t0 ) = x0 ,

x(t) = M (t)M −1 (t0 )x0 .

By Theorem 4.11, we also have

Theorem 5.10. A solution matrix M (t) of Eq. ( 3.1) on I is a fundamental matrix if and only
if ddetM(t) 6= 0 (or WQDE (t) 6= 0) on I. Moreover, for some t0 in I such that ddetM(t0 ) 6= 0
(or WQDE (t0 ) 6= 0 ), then ddetM(t) 6= 0 (or WQDE (t) 6= 0 ).

Example 1 Show that

Rt
e−it · ejt dt
 
T eit eit t0
M (t) = (φ1 (t), φ2 (t)) =
0 ejt

is a fundamental matrix of the QDES

    
ẋ1 i 1 x1
= . (5.5)
ẋ2 0 j x2

 Proof.
it
 First, we show that M (t) is a solution matrix of Eq.(5.5). In fact, let φ1 (t) =
e
, in view of (it)0 (it) = (it)(it)0 , then
0
 it     it   
ie i 1 e i 1
φ̇1 (t) = = = φ1 (t),
0 0 j 0 0 j
Rt
e−it · ejt dt
 
eit t0
which implies that φ1 (t) is a solution of Eq.(5.5). Similarly, let φ2 (t) = ,
ejt
then
Rt Rt
e−it · ejt + ejt dt e−it · ejt dt
      
ieit t0
i 1 eit t0
i 1
φ̇2 (t) = = = φ2 (t),
jejt 0 j ejt 0 j

30
which implies that φ2 (t) is another solution of Eq.(5.5). Therefore, M (t) = (φ1 (t), φ2 (t))> is
a solution matrix of Eq.(5.5).

Secondly, by Theorem 5.10, and taking t0 = 0,

ddetM(t0 ) = rdet(M(t0 )M+ (t0 )) 6= 0,

then ddetM(t) 6= 0. Therefore, M (t) is a fundamental matrix of Eq. (5.5).

From Theorem 5.9 and Theorem 5.10, the following corollaries follows immediately.

Corollary 5.11. Let M (t) be a fundamental matrix of Eq. ( 5.5) on I. Q ∈ H2×2 is an

invertible constant quaternionic matrix. Then M (t)Q is also a fundamental matrix of Eq.
( 5.5) on I.

PROOF. As pointed out in Remark 5.8, any solution matrix M (t) satisfies Eq. (3.1),
that is,
Ṁ (t) = A(t)M (t).
The inverse also holds. Taking Ψ(t) = M (t)Q, differentiating it, and noticing that M (t) is a
fundamental matrix satisfying Eq. (3.1), we have

Ψ0 (t) = M 0 (t)Q = A(t)M (t)Q = A(t)Ψ(t),

which implies that Ψ(t) is a solution matrix of Eq. (3.1). Moreover, since Q is invertible, we
have
ddetΨ(t) = ddetM(t)Q = ddetM(t) · ddetQ 6= 0.
Therefore, Ψ(t) = M (t)Q is also a fundamental matrix of Eq. (3.1).

Corollary 5.12. If M (t), Ψ(t) are two fundamental matrices of Eq. ( 3.1) on I, there exists
an invertible constant quaternionic matrix Q ∈ H2×2 such that

Ψ(t) = M (t)Q, t ∈ I.

PROOF. Since M (t) is a fundamental matrix, M (t) is invertible and denote the inverse
by M −1 (t). Set
Ψ(t) = M (t)X(t), t ∈ I. (5.6)
M (t), Ψ(t) are fundamental matrices, so ddetM (t) 6= 0 and ddetΨ(t) 6= 0. Thus, ddetX(t) 6= 0.
Differentiation on (5.6) gives

Ψ0 (t) = M 0 (t)X(t) + M (t)X 0 (t) = A(t)M (t)X(t) + M (t)X 0 (t) = A(t)Ψ(t) + M (t)X 0 (t)

On the other hand, Ψ(t) is a fundamental matrix, which implies that

Ψ0 (t) = A(t)Ψ(t).

31
Compare above two equalities, we have

M (t)X 0 (t) = 0.

Noticing that M (t) is invertible, we obtain that X 0 (t) = 0, which implies that X(t) is a
constant matrix. Letting X(t) = Q. The proof is complete.

Now consider the following quaternion-valued equations with constant quaternionic coeffi-
cients
ẋ(t) = Ax(t), (5.7)
where A ∈ H2×2 is a constant quaternion matrix. Then we have
Theorem 5.13. M (t) = exp{At} is a fundamental matrix of ( 5.7). Moreover, any solution
x(t) of ( 5.7) can be represented by

x(t) = exp{At}q,

where q ∈ H2 is an arbitrary constant quaternion. For the IVP x(t0 ) = x0 , any solution x(t)
of ( 5.7) can be represented by
x(t) = exp{A(t − s)}x0 .

Proof. Note that A ∈ H2×2 is a constant quaternion matrix, by using At·(At)0 = (At)0 ·At
we have
M 0 (t) = [exp{At}]0 = A exp{At} = AM (t),
which implies exp{At} is a solution matrix of (5.7). Moreover, ddetM(0) = ddet exp(A · 0) =
ddetE = 1 6= 0. Consequently, ddetM(t) 6= 0. Therefore, exp{At} is a fundamental matrix of
(5.7).

Now consider the diagonal homogenous system.

    
ẋ1 (t) a1 (t) 0 x1 (t)
= (5.8)
ẋ2 (t) 0 a2 (t) x2 (t)

Theorem 5.14. Assume that

Z t Z t
ai (t) ai (s)ds = ai (s)dsai (t) (5.9)
t0 t0

holds. The fundamental matrix can be chose by

Rt !
exp{ t0 a1 (s)ds} 0
M (t) = Rt .
0 exp{ t0 a2 (s)ds}

Then the solution of the diagonal system ( 5.8) with the initial value x(t0 ) = x0 is given by

x(t) = M (t)x0 .

32
6 Algorithm for computing fundamental matrix

Since the fundamental matrix plays great role in solving QDEs, in this section, we provide an
algorithm for computing fundamental matrix of linear QDEs with constant coefficients.

6.1 Method 1: using expansion of exp{At}

From Theorem 5.13, we know that exp{At} is a fundamental matrix of (5.7). So if the
coefficient matrix is not very complicate, we can use the definition of exp(At) to compute
fundamental matrix of linear QDEs with constant coefficients.
Theorem 6.1. If A = diag(λ1 , λ2 ) ∈ H2×2 is a diagonal matrix, then
 
exp{λ1 t} 0
exp{At} = .
0 exp{λ2 t}

Proof. By the expansion,

   2 2  
λ1 0 t λ1 0 t exp{λ1 t} 0
exp{At} = E + + + ··· = .
0 λ2 1! 0 λ2 2! 0 exp{λ2 t}

By Proposition 5.5, we can divide the matrix to some simple ones and use the expansion
to compute the fundamental matrix.
A = diagA + N,
where N is a nilpotent matrix. That is, N n = 0 and n is a finite number.

Example 2 Find a fundamental matrix of the following QDES

 
k 1
ẋ = x, x = (x1 , x2 )T .
0 k

Answer. We see that

     
k 1 k 0 0 1
= + .
0 k 0 k 0 0
     
k 0 0 1 0 1 k 0
Noticing that = , by Proposition 5.3 and Theorem
0 k 0 0 0 0 0 k
6.1, we have
   
k 0 0 1
exp(At) = exp exp
0 k 0 0
 kt     2
e 0 0 1 t 0 1 t2

= E + + + · · · .
0 ekt 0 0 1! 0 0 2!

33
Note that  2    3  
0 1 0 0 0 1 0 0
= , = ,··· ,
0 0 0 0 0 0 0 0
Then the fundamental matrix
    
ekt 0 1 t ekt tekt
exp{At} = = .
0 ekt 0 1 0 ekt
 
λ 1
For the matrix , λ ∈ H, similar computation shows that
0 λ
Theorem 6.2. For any λ ∈ H, we have
    λt  
λ 1 e 0 1 t
exp t = .
0 λ 0 eλt 0 1

If we generalize this result to n-dimensional case, we have the following formula.

Theorem 6.3. For any λ ∈ H, we have
 
t2 k−2
t tk−1
  

 λ 1 0···0 0 
 1 t 2!
· · · (k−2)! (k−1)!
1···0 tk−3 tk−2
 0 λ 0 
   

 
 

  0 1 t · · · (k−3)! (k−2)!

exp 
 ... ... 
 t = diag(exp λt, · · · , exp λt) 
 .. ..

.
     . . 




 0 0 0···λ 1  




 0 0 0···1 t


 0 0 0···0 λ 
0 0 0···0 1

We remark that this method can not be applied extensively when the coefficient matrix A
is complicate. Usually, the two divisions of A can not be communicate. For example,
     
i k i 0 0 k
= + .
0 j 0 j 0 0
We see that      
i 0 0 k 0 k i 0
6=
0 j 0 0 0 0 0 j
 
i 1
In this case, this method can not be used. Even, for a simple matrix , you can
0 j
not use this method. So we should find more effective method to compute the fundamental
matrix. In what follows, we try to use the eigenvalue and eigenvector theory to compute it.

6.2 Method 2: eigenvalue and eigenvector theory

Since quaternion is a noncommunicative algebra, eigenvalue of a matrix A ∈ H2×2 should be

defined by left eigenvalue and right eigenvalue, respectively. A quaternion λ is said to be a
right eigenvalue of A if
Ax = xλ

34
for some nonzero (column) vector x with quaternion components. Similarly, a quaternion λ is
a left eigenvalue of A if
Ax = λx
for some nonzero (column) vector x with quaternion components. Right and left eigenvalues
are in general unrelated. Usually, right and left eigenvalues are different. It should be noted
that there are many differences in the eigenvalue problem for complex (real) and quaternion
matrices. Eigenvalues of complex matrices satisfy Brauer’s theorem [5] for the inclusion of
the eigenvalues. While right eigenvalues of quaternion matrices do not have this property.
Another contrast in the eigenvalue problems for complex and quaternion matrices is that a
complex n×n matrix cannot have more than n complex eigenvalues, while it can have infinitely
many quaternion left eigenvalues (see Theorem 2.1 in [17]). Finally, it is well-known that any
eigenvalue of an n × n complex matrix must be a root of its characteristic polynomial. Left
eigenvalues and right eigenvalues of quaternion matrices do not share this property of complex
matrices.

In this paper, we focus on the right eigenvalue and eigenvector. Because we emphasize on
finding the solution taking the form

x = qeλt , q = (q1 , q2 )T ,

where λ is a quaternionic constant and q is a constant quaternionic vector. Substituting it

into Eq. (5.7), we have
qλeλt = Aqeλt .
Because eλt 6= 0, it follows that
qλ = Aq,
or
Aq = qλ. (6.1)
So if we find such eigenvalue λ and eigenvector q in (6.1), we will find the solution of Eq. (5.7).
We also say that the eigenvalue λ and eigenvector q are the right eigenvalue and eigenvector
of the QDEs (5.7).

Notice also that if 0 6= α ∈ H, then

Aq = qλ ⇒ Aqα = qα(α−1 λα),

so we can sensibly talk about the eigenline spanned by an eigenvector q, even though there
may be many associated eigenvalues! In what follows, if

θ = α−1 λα

we call it that θ is similar to λ. If there exists a nonsingular matrix T such that

A = T −1 BT,

we call it that A is similar to B, denoted by A ∼ B.

35
Remark 6.4. If λ is a characteristic root of A, then so is α−1 λα.

From the definition of fundamental matrix, we have

Theorem 6.5. If the matrix A has two right linearly independent eigenvectors q1 and q1 ,
corresponding to the eigenvalues λ1 and λ2 (λ1 and λ2 can be conjugate), then

M (t) = (q1 eλ1 t , q2 eλ2 t )

is a fundamental matrix of Eq. ( 5.7).

Proof. From above discussion, we know q1 eλ1 t , q2 eλ2 t are two solution of Eq. (5.7). Thus,
M (t) is a solution matrix of Eq. (5.7). Moreover, by using the linear independence of q1 and
q2 , we have
ddetM (0) = ddet(q1 , q2 ) 6= 0.
Therefore, M (t) is a fundamental matrix of Eq. (5.7).

Now we need a lemma from ([4] Proposition 2.4).

Lemma 6.6. Suppose that λ1 , λ2 , · · · , λr are distinct eigenvalues for A, no two of which
are similar, and let q1 , · · · , qr be corresponding eigenvectors. Then q1 , · · · , qr are linearly
independent.

Then by Lemma 6.6, we have the corollary

Corollary 6.7. If the matrix A has two distinct eigenvalues λ1 and λ2 (λ1 and λ2 can not be
conjugate), then
M (t) = (q1 eλ1 t , q2 eλ2 t )
is a fundamental matrix of Eq. ( 5.7).

Now we need some results from ([6] Theorem 11, Theorem 12).

Lemma 6.8. If A is in triangular form, then every diagonal element is a characteristic root.

Lemma 6.9. Let a matrix of quaternion be in triangular form. Then the only characteristic
roots are the diagonal elements (and the numbers similar to them).

Lemma 6.10. Similar matrices have the same characteristic roots.

Lemma 6.11. ([6] Theorem 2) Every matrix of quaternion can be transformed into triangular
form by a unitary matrix.

36
Example 3 Find a fundamental matrix of the following QDEs
 
i j
ẋ = x, x = (x1 , x2 )T . (6.2)
0 i+j

Answer: From Lemma 6.9 and Lemma 6.10, we see that λ1 = i and λ2 = i + j. To find
the eigenvector of λ1 = i, we consider the following equation

Aq = qλ1 ,

that is 
iq1 + jq2 = q1 i
(6.3)
(i + j)q2 = q2 i
From the second equation of (6.3), we can take q2 = 0. Substituting it into the first equation
of (6.3), we obtain q1 = 1. So we can take one eigenvector as
   
q1 1
ν1 = =
q2 0

To find the eigenvector of λ2 = i + j, we consider the following equation

Aq = qλ2 ,

that is 
iq1 + jq2 = q1 (i + j)
(6.4)
(i + j)q2 = q2 (i + j)
We can take one eigenvector as
   
q1 1
ν2 = =
q2 1

Since   h 1 1  1 0 i
1 1
ddet(ν1 , ν2 ) = ddet = rdet 6= 0,
0 1 0 1 1 1
the eigenvectors ν1 and ν2 is linearly independent. Taking
 it (i+j)t 
λ1 t λ2 t e e
M (t) = (ν1 e , ν2 e ) = ,
0 e(i+j)t

from Theorem 6.1, M (t) is a fundamental matrix. In fact, by the definition of fundamental
matrix, we can verify that M (t) is a fundamental matrix of Eq. (6.2), which is consistent with
Theorem 6.1. Now we verify the fundamental matrix as follows. First, we show that M (t) is
a solution matrix of Eq. (6.2). Let φ1 (t) = ν1 eλ1 t and φ2 (t) = ν2 eλ2 t , then
      
1 it i j 1 it i j
φ̇1 (t) = ie = e = φ1 (t),
0 0 i+j 0 0 i+j

37
  
e(i+j)t
which implies that φ1 (t) is a solution of Eq. (6.2). Similarly, let φ2 (t) = , then
e(i+j)t
      
1 (i+j)t i j 1 (i+j)t i j
φ̇2 (t) = (i + j)e = e = φ2 (t),
1 0 i+j 1 0 i+j

which implies that φ2 (t) is another solution of Eq. (6.2). Therefore, M (t) = (φ1 (t), φ2 (t))T is
a solution matrix of Eq. (6.2).

Secondly, by Theorem 5.10, and taking t0 = 0,

ddetM(t0 ) = ddetM(0) = ddet(ν1 , ν2 ) 6= 0,

then ddetM(t) 6= 0. Therefore, M (t) is a fundamental matrix of Eq. (6.2).

7 Conclusion and Discussion

Conclusion

This paper established a systematic basic theory for two-dimensional linear QDEs, which
have many applications in the real world such as quantum mechanics, differential geometry,
Kalman filter design, attitude dynamics, fluid mechanics, and so on. For the sake of un-
derstanding, we focus on the 2-dimensional system. However, in fact, our basic results can
be easily extended to any arbitrary n-dimensional QDEs. Due to the non-commutativity of
quaternion algebra, there are four profound differences between QDEs and ODEs.

1. Due to the non-commutativity of the quaternion algebra, the algebraic structure of the
solutions to QDEs is not a linear vector space. It is actually a right-free module.

2. Wronskian of ODEs is defined by Caley determinant. However, since Caley determinant

depends on the expansion of i-th row and j-th column of quaternion, different expansions can
lead to different results. Wronskian of ODEs can not be extended to QDEs. It is necessary to
define the Wronskian of QDEs by a novel method (see Section 4). We use double determinants
ddetM = rddet(MM+ ) to define Wronskian of QDEs. We remark that the Wronskian of QDEs
can also be defined by ddetM = rddet(M+ M).

3. Liouville formula for QDEs and ODEs are different.

4. For QDEs, it is necessary to treat the eigenvalue problems with left and right, separately.
This is a large-difference between QDEs and ODEs.

38
8 Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this
article.

9 Acknowledgement

The authors would express their great gratitude to Professor Weinian Zhang for his valuable
comments and discussion when they prepared this paper.

References
[1] V. I. Arnold, Ordinary Differential Equations, Cambridge:M. I. T. Press, (1973).

[2] S. L. Alder, Quaternionic quantum field theory, Commun. Math. Phys. 104, 611-656,
(1986).

[3] S. L. Adler, Quaternionic Quantum Mechanics and Quantum Fields, Oxford University
Press, New York, 1995.

[4] A. Baker, Right eigenvalues for quaternionic matrices: A topological approach, Linear
Algebra Appl., 286 (1999), pp. 303-309.

[5] A. Brauer, Limits for the characteristic roots of matrices II, Duke Math. J. 14 (1947), pp.
21-26.

[6] J.L. Brenner, Matrices of quaternions, Pac. J. Math. 1 (1951) 329-335.

[7] J. Campos and J. Mawhin, Periodic solutions of quaternionic-values ordinary differential

equations, Annali di Mathmatica, 185 (2006), S109-S127.

[8] A. Cayley, On certain results relating to quaternions, Philos. Mag. 26 (1845), 141-145;
reprinted in The Collected Mathematical Papers Vol. 1, Cambridge: Cambridge 39. and
invariant real University Press, 1989, pp. 123-126.

[9] L. Chen, Definition of determinant and Cramer solution over the quaternion field, Acta
Math. Sinica (N.S.) 7 (1991), 171-180.

[10] L. Chen, Inverse matrix and properties of double determinant over quaternion field, Sci.
China Ser. A, 34 (1991), 528-540.

[11] C. Chicone, Ordinary Differential Equations, Texts in Applied Mathematics, Volume 34

2006, Springer.

39
[12] Chow, S. N. and J. K. Hale, Methods of Bifurcation Theory. New York: Springer-
Verlag(1982)..

[13] R. Dimitrid and B. Goldsmith, Sir William Rowan Hamilton, Math. Intelligencer, 11
(1989), no. 2, 29-30

[14] T. Ding, C. Li, Ordinary Differential Equations, Higer Education Press, Beijing, 2004.

[15] J. R. Dooley, J. M. McCarthy, Spatial rigid body dynamics using dual quaternion conipo-
nents, Proceedings of the 1991 IEEE Intemational Conference on Robotics and Automa-
tion Sacramento, Califomia- April 1991.

[16] L. P. Eisenhart. A Treatise on the Differential Geometry of Curves and Surfaces. Dover,
New York, 1960. Originally published in 1909.

[17] F. Farid, Q. Wang, F. Zhang, On the eigenvalues of quaternion matrices, Linear and
Multilinear Algebra, 59:4, 451-473.

[18] Harley Flanders. Differential Forms with Applications to Physical Sci- ences. Academic
Press, New York, 1963.

[19] Alfred Gray. Modern Differential Geometry of Curves and Surfaces. CRC Press, Inc.,
Boca Raton, FL, 1993.

[20] J.D. Gibbon, A quaternionic structure in the three-dimensional Euler and ideal magne-
tohydrodynamics equation, Physica D 166 (2002), 17-28.

[21] J.D. Gibbon, D.D. Holm, R.M. Kerr and I. Roulstone, Quaternions and particle dynamics
in the Euler fluid equations, Nonlinearity, 19 (2006), 1969-1983.

[22] A. Gasull. J. Llibre, X. Zhang, One-dimensional quaternion homogeneous polynomial

differential equations, J. Math. Phys., 50 (2009), 082705.

[23] J. K. Hale, Ordinary Differential Equations, 2nd ed. Malabar: R. E. Krieger Pub. Co.,
(1980).

[24] T. Isokawa, T. Kusakabe, N. Matsui, F. Peper, Quaternion neural network and its appli-
cation, Lecture Notes in Computer Science, vol 2774, 318-324.

[25] N. Jacobson, Lecture in Abstract Algebra, Springer-Verlag, New York, 1951.

[26] I.I. Kyrchei, Cramer’s rule for some quaternion matrix equations, Appl. Math. Comput.,
217 (2010) 2024-2030.

[27] S. Leo, G. Ducati, C. Nishi, Quaternionic potentials in non-relativistic quantum mechan-

ics, J. Phys. A., 11/2001; 35(26).

[28] S. Leo, G. Ducati, Quaternionic wave packets, J. Math. Phys. 48, 052111 (2007).

40
[29] S. Leo, G. Ducati, Solving simple quaternionic differential equations, J. Math. Phys.
07/2003; 44(5).

[30] A. Handson, H. Hui, Quaternion frame approach to streamline visualization, IEEE Tran.
Vis. Comput. Grap., 1(1995), 164-172.

[31] P. Hartman, Ordinary Differential Equations. New York: John Wiley & Sons, Inc., (1964).

[32] J. Marins, X. Yun, E. Bachmann, R. McGhee, and M. Zyda, An Extended Kalman

Filter for Quaternion-Based Orientation Estimation Using MARG Sensors, Proceedings
of the 2001 IEEE/RSJ International Conference on Intelligent Robots and Systems Maui,
Hawaii, USA, Oct. 29- Nov. 03, 2001.

[33] J. Proškovǎ, Description of protein secondary structure using dual quaternions, Journal
of Molecular Structure, 1076 (2014) 89-93.

[34] CH Robotic Lab, Understanding Euler Angles, available online at:

http://www.chrobotics.com/library/understanding-euler-angles.

[35] CH Robotic Lab, Understanding Quaternions, available online at:

http://www.chrobotics.com/library/understanding-quaternions.

[36] Marsden, J. E, Sirovich, L, Golubitsky, M. and Jäger, W., Ordinary Differential Equations
with Applications, Springer-Verlag, New York, 1999.

[37] Rubtsov, V. N. and Roulstone, I., Examples of quaternionic and Keller structures in
Hamiltonian models of nearly geostrophic flow, J. Phys. A 30, L63,(1997).

[38] Roubtsov, V. N. and Roulstone, I., Holomorphic structures in hydrodynamical models of

nearly geostrophic flow, Proc. R. Soc. London, Ser. A 457, 1519, (2001).

[39] F. Udwadia, A. Schttle, An Alternative Derivation of the Quaternion Equations of Motion

for Rigid-Body Rotational Dynamics, J. Appl. Mech.77:4(2010).

[40] J.R. Wertz, Spacecraft Attitude Determination and Control, Kluwer Academic Publish-
ers, The Netherlands(1978).

[41] P. Wilczynski, Quaternionic-valued ordinary differential equations. The Riccati equation,

J. Differential Equations 247 (2009), 2163-2187.

[42] X. Zhang, Global structure of quaternion polynomial differential equations, Commun.

Math. Phys., 2011,303(2), pp 301-316.

[43] F. Zhang, Y. Wei., Jordan cannonical form of a oartitional complex matrix and its appli-
cation to real quaternion matrix, Commun. Algebra, 29(6) (2001), 2363-2375.

[44] W. Zhang, Z. Du, B. Xu., Ordinary differential equations, Higher Education Press, Bei-
jing, 2006 (in Chinese).

41