Ordinary Differential Equations (ODEs) has a broad range of applications in other subjects, such as Physics and Differential Geometry. Under specific circumstances, the problems encountered in these related areas are extremely difficult to solve using real numbers as parameters. To solve higher-dimensional cases, we need to use constant matrices as parameters (H. Cartan, 1981) . Ordinary differential equations (ODEs) are fundamental for the study for differential equations with multiple variables (Sheldon Axler, 2015) . Although linear ODEs have a comparatively easy form, they are effective in solving certain physical and geometrical problems (Zoltan Vizvari, 2020) . There exist various researches combining matrix with linear ODEs, yet most are purely theoretical and very few of them aim to discover the application of such a computational method in simplifying calculations in other related subjects.
In this article, we will explore the applications of the matrix method in linear ordinary differential equations (linear ODEs) in Physics and other branches of Mathematics. We will begin by introducing fundamental knowledge in Linear Algebra and proving the existence and uniqueness of solution for ODEs. Then, we will concentrate on finding the solutions for ODEs and introducing the matrix method for solving linear ODEs. Eventually, we will apply the conclusions we’ve gathered from the previous parts into solving problems concerning Physics and differential curve. Our ultimate aim is to be able to solve any given linear ordinary equations using constant matrices as parameters, and apply such a computational method into solving the evolution function of the movement of a free electron in a constant three-dimensional electric field and proving the Frenet formula in calculating the torsion of a 3-D curve (Fage M.K., 1974) .
Let us begin with some basic examples of how ODEs can be used to solve physical problems.
Example 1 (Radioactive chains of decay)
The differential equation for the number N of radioactive Nucleus is
where k is the decay constant:
1) Give the evolution function with initial condition .
2) We now look at a chain, where the original nucleus decays into another radioactive nucleon. Let α, β, γbe three types of nucleus. αtransfers to βwith decay constant , whereas βtransfers to γwith decay constant . Assume that there are N αin the beginning, find the evolution functions of the number of nucleus α, βand γ.
3) If , show that after a long term evolution, the ratio of the number of αand that of βtends to
Example 2 (Motion in the liquid)
In this question, we consider a ball sinking in a liquid under the influence of gravity. There are three forces horizontally acts on this body: the gravity, the buoyancy and the viscous force.
1) Write the differential equation for this motion and solve it. One may dispute whether the gravity is greater than the buoyancy.
2) Find the ultimate velocity.
The first two examples can be solved mainly through using separation of variable and the resolving kernel for linear ODE. The third example is, however, difficult to solve because the magnetic field B is a three-dimensional matrix instead of a constant. To find the answer for a problem like this, we need to introduce the matrix method for solving linear ODE.
Example 3 (Motion of a charged particle in an electromagnetic field)
A charged particle with mass m and charge q moves in electromagnetic field with a constant magnetic field B and a changing electronic field . We make the following assumption that the rate of the change of is small enough that the magnetic field it generates is negligible with respect toB. The evolution equation is
where m signifies the mass of the particle, v the velocity, q the charge,B the magnetic field and the electronic field depending with time, with the initial condition .
1) What is the resolving kernel of this linear differential equation?
2) What is the evolution function of the velocity, i.e., find .
3) Find the trajectory of this particle,i.e., find .
From the above examples, especially the last one, we see that even if the differential equation is relatively easy, we do not know in a general way how to solve it specifically if we cannot integrate the equation directly. Our goal of this article is to be able to solve the differential equation of the type of example 3.
2. Linear Algebra
2.1. Vector Spaces and Linear Maps
Definition 1 (Vector spaces) Aset V containingtwooperations that satisfy the eight axioms listed belowis considered a vector space (math world) :
Scalar Multiplication :
To be a vector space,all elements X,Y,Z inVand any scalars r,s inR must satisfy the following axioms:
3) For allX, .
4)For anyX,there exists a −X such that .
Definition 2 (Basis) We consider vectors a basis of if:
1) generate .
2) are linearly independent.
Definition 3 A map of vector spaces is a linear map,if , , .
2.2. Relations of Matrices and Linear Maps
Having given the definition of linear maps and matrices, we can show that the information in a linear map can be completely depicted through a matrix given a chosen basis.
What we would like to find out is , given . First, we find a basis , and write X under this basis:
Then, since for every i, to find , all we need to know is a . When is given for every , every is found, so that we can depict for every i through matrix:
Therefore, can be expressed as:
which tells us that given a basis, a linear map can be completely depicted through a matrix.
Additionally, linear map and matrix satisfy the following relations:
1) Linear map can be represented by matrix
2) Linear maps can be represented by matrix
3) Linear maps can be represented by matrix
4) Linear maps can be represented by matrix
2.3. Exponential Map of Matrices
In this subsection, we are going to give the definition of the exponential map of matrices, which is a generalization of exponential of real numbers. We will also present some basic properties of exponential maps of matrices.
Firstly, we know that for a complex number , ex can be expressed using the sum of an infinite Taylor series:
Then, we present some basic properties of the exponential map :
Having got these properties, we now substitute “x” in with a matrix A, so that . After the substitution, the exponential map for
matrix still satisfies most of the properties mentioned above, yet specifically, we need to notice here that for most matrices A and B, since the multiplication of matrices does not fit the commutative law. Our aim now is to be able to calculate using a relatively easy approach.
Given a matrix, its exponential is generally hard, even impossible to calculate just by following its definition. A procedure of calculating the exponential map will be given after the introduction of Jordan normal forms (JNF).
2.4. Jordan Normal Forms
We are going to admit the following theorem of the existence of Jordan normal forms.
Theorem 1 For any n × n matrix A,there exists a matrix
such that A is similar to . Here
and are not necessarily distinct.
For a proof, see for example .
In this subsection, we will present how to find the Jordan normal form and the transition matrix for a given matrix. The finding of JNF and transition matrix follows the following basic procedure:
1) Calculate the characteristic polynomial of the given matrix.
2) Let the characteristic polynomial equal to zero, and calculate the eigenvalues,
3) Put the eigenvalues on the main diagonal, and decide the numbers of “1s” on the up-right corner of every small Jordan matrix through comparing the number of set of eigenvectors that the potential solutions have.
4) Having found JNF, use to calculate the transition matrix, then find its inverse.
To clarify the process, we should now analyze an example: to calculate the Jordan normal form of matrix and transition matrix Q of : Firstly, calculate the characteristic polynomial, and let it equal zero :
Then solve the eigenvalues, we have . Therefore, the Jordan normal form of A is either or
Then, since Av = 2v has only one solution, matrix A has only one set of eigenvector. The matrix , however, has two eigenvectors. Therefore, only can be the JNF of A.
Having got JNF, we use to calculate the transition matrix. Such an equation can be solved using method of undetermined coefficients. Let Q be , and we can easily get a = 1; b = 0; c = 1; d = 1. Therefore, the transition matrix is .
2.5. Calculation of the Exponential Map of Matrices
Besides the properties of exponential maps mentioned above, we would like to further prove that ; where P is a matrix,Q its transition matrix, and the inverse of Q (Chantladze T., 2000) .
Having got such a property, we could then calculate the given example: , where
To simplify our calculation, we have to first find the Jordan normal form of A as well as the transition matrix and its inverse matrix. As has been calculated in 2.4,
According to the definition of matrix map,
Therefore, can be represented as:
Basically, the calculation of follows the following procedure:
Step 1: Calculate the JNF of A and find the transition matrix Q.
Step 2: Calculate
The example presented here is relatively easy, but it presents clearly the procedure of calculating the exponential maps of matrices. In the following section, we will utilize the exponential map of matrices to develop the matrix method for solving linear ordinary differential equations.
3. Ordinary Differential Equations
In this section, we are going to prove the central results of ODE theory, the Picard-Lindelöf theorem. It is on the existence and uniqueness of the initial value problems:
where F is a Lipschitz function with respect to vector x. For the proof reference, see . This theorem is quite useful in theory. For example, in Section 6, we are going to present a direct application of this theorem to differential geometry.
According to the definition, is continuous in a bounded region , and satisfies in U the Lipschitz condition in , where k is independent of . Find such that we have , then let R be defined as a rectangle such that . What we would like to prove now is that such that , the 1st order differential equation has one unique solution for which .
First, we need to define the following functions:
and we want to prove that is the solution we are looking for.
To begin with, we need to show that is bounded by the rectangleR, where . Assume that is bounded by R, we have:
which tells us that . Therefore, since is clearly bounded byR, and such a conclusion holds for , we can conclude that is bounded by R, through induction.
Next, we will prove
Suppose that this inequality holds for , we will have . According to the Lipschitz condition, we also got . Therefore, we have
which tells us that
When n = 1, such a conclusion holds true:
so according induction, the inequality
Then, we will prove that for , sequence converges. Based on , we have
Since converges absolutely , we know that
is converges uniformly for . Therefore, is continuous and converges uniformly for .
And now we can prove the existence of the solution by showing satisfying . Because converges uniformly on and , we know that tends uniformly to . Let n approach infinity for
we can get
Since is continuous, it has as its derivative and .
Having proved the existence of the solution , we will further prove its uniqueness given . Let us assume that there exists another different solution . When , let , and we have
hence . We can successively obtain the upper bound of on by repeating such argument, which eventually gives us
which approaches zero. Therefore, on interval , which contradicts the assumption that is different from , and this proves the uniqueness of solution .
4. Matrix Method for Linear ODE
The precedent section concerning the existence and uniqueness of solutions of ODEs with given boundary condition is quite useful theoretically. It has the disadvantage of not being able to give an explicit expression of the solution, though, which is demanded in many physical problems. In this section, we are going to focus on a special kind of ODEs: the linear ODEs and give an explicit expression of solutions using the “resolving kernel” (Halas Zdenek, 2005) .
4.1. Linear ODEs
An nth degree ordinary differential equation is called a nth degree linear ordinary differential equation if the variable function x and all of its nth order derivative are of first degree. Its normal for can be written as
where and are continuous functions on their domains. The linear ordinary differential equation is called homogeneous if , and a homogeneous linear ODE has the following 2 properties:
1) If are solutions for a homogeneous linear ODE, then ( is a constant) is also a solution.
2) If is a solution for such a homogeneous differential equation, and , then .
The linear ODE can also be written in the following form:
where is a continuous matrix, a continuous function of x, and an initial value. Such an ODE, as have been proven in section three, has a unique solution. The approach to find its solution will be thoroughly discussed in section 4.2 and 4.3.
4.2. Resolving Kernel
Definition 4 (Resolving kernel) Let be the unique solution of thefollowing linear ODE:
is called the Resolving Kernel of the linear ODE:
Property 1: if
Proof: Given , , by the definition of , is the unique solution of:
By the definition of is the unique solution of:
Therefore, based on the differential equations above, .
Property 2: is reversible
Proof: We try to inverse :
Therefore, is reversible. The inverse is .
Property 3: , where and . Det() means matrix determinant.
Proof: Let be a basis of , then let , where , :
Since , and :
4.3. Solving Linear ODEs
After showing the definition of Resolving Kernel and several of its properties, we are now going to apply it into solving certain linear ODEs:
How to solve:
X can be written as , and can be written as , which equals to:
Let be the resolving kernel of :
And since , where can be expressed as
, as , and so on.
where can be expressed as
4.4. Linear ODEs with Constant Coefficients
Given , where matrix is constant, we will prove the following theorem:
Theorem 2 Let be the resolving kernel of the ODE:
can then be expressed as
5. Applications in Physics
Having introduced basic knowledge about linear algebra and ODE and the matrix method for linear ODE, we are now going to apply them in solving Physics problems presented in the introduction:
1) (Radioactive chains of decay) The differential equation for the number Nof radioactive Nucleus is
where k is the decay constant. This is a problem that can be solved through ordinary approaches. To solve N, all we have to do is to apply integration after separating the variables, and we can get: .
Then, we will discuss some more complex circumstances. Now look at a chain, where the original nucleus decays into another radioactive nucleon. Let be three types of nucleus. transfers to with decay constant whereas transfers to with decay constant . Assume that there are N in the beginning.
Based on the information given, we can list the following differential equations which characterize the rate of change each kind of nucleus:
The first equation can be solved easily through the approach mentioned at the very beginning. The result is:
The second equation can be solved using
where , , , . Such an approach is usually utilized in solving linear ODE and will be thoroughly introduced in the third section in this article:
Having got can be easily calculated through integration:
so far, the three chains are all solved. Then, to evaluate the ratio of the number of α and that of β after a longtime of evolution, we need to calculate:
Let time t approach infinity, we get:
So far, this example is thoroughly solved.
2) (Linear motion of a particle in liquids) In this question, we consider a ball sinking in a liquid under the influence of gravity. There are three forces horizontally acts on this body: the gravity, the buoyancy and the viscous force.
The first thing we would like to do is to find the differential equation that characterizes the motion. The viscous force, as is mentioned in the question stem, is in direct proportion to the velocity v. The gravitational force points downwards, and the viscous force and buoyancy point upwards. Therefore, according to
newton’s second law, the net force can be expressed as ,
which equals to , where is the mass, rho the density of the liquid, V the volume of the ball, g the gravitational constant, v the velocity, and a constant. Substituting with , where is a constant, we can get:
Then, since the equation is a first order linear ODE, it can be solved using the approach is Example 1:
where , and represents the initial velocity. This differential equation can be solved in the following way:
So far, the differential equation which characterizes the motion of the ball has been solved. To find the ultimate velocity of the ball, let time t approach infinity, and we can get:
3) (Harmonic oscillations) A harmonic oscillation is a linear movement (along the axis), where the resulting force is always directed against and proportional to the distance to the position of the equilibrium.
We will first assume that there is no friction force on the oscillating object. Since the force is always directed against and proportional to the distance to the position of the equilibrium, the force can be expressed as , where k is a constant and x represents the equation for displacement. The we will have the following equation:
where represents the initial position and represents the initial velocity.
We will then use the matrix method directly to solve this equation as practice, though this might make the calculation more difficult. First, let be the matrix representing the displacement and velocity:
Then we write the resolving kernel of this equation:
Since and we can solve using its series expansion:
Having got , we can calculate using :
Then, to consider a more complex situation, we assume that there exists a friction which is also proportional to the velocity and directed opposite to the velocity. We can rewrite the differential equation that characterizes the motion as
write so that , where
Since now becomes more complex, the approach previously used in the non-frictional circumstance is no more applicable. To calculate there solving kernel, we need to use the method mentioned in section 2, to first find the Jordan Normal Form (JNF) and the transition matrix of matrix A. The first step is
to solve :
and we can get the solutions . Now that can either be or . Since the previous equation has two solutions, we know that should be the first one, since the first matrix has two eigenvalues.
Having got , we can then calculate the transition matrix Q by using . The value of Q is thus . Suggest that , can be calculated using , can we can get the result
Then, can be expressed as
so that the displacement equation , where
Now that there are two kinds for situations we need to discuss. If , ,
However, if , the equals
Since the final solution to can be written as
Having solved the three questions above, we’re now going to take a look at the one concerning electromagnetism, which cannot be solved without utilizing the matrix method for linear ODE:
Example 4 (Motion of a charged particle in an electromagnetic field)
A charged particle with mass m and charge q moves in an electromagnetic field with a constant magnetic fieldB and a changing electronic field . We make the following assumption that the rate of the change of is small enough that the magnetic field it generates is negligible with respect to B. The evolution equation is
where m signifies the mass of the particle, v the velocity, q the charge, B the magnetic field and the electronic field depending with time, with the initial condition . Here, we assume that the Matrix B is
To solve such a differential equation, we need to first solve the Jordan normal form of B and the transition matrix as well as its inverse matrix. The first step is to calculate the eigenvalues of B: we have to first calculate the determinant of , where λ suggests each eigenvalue, I the identity:
Let , we get , where . Therefore, the eigenvalues of B are 0, and Since three eigenvalues of B are different, the matrix B is a diagonalizable matrix. Therefore, the Jordan normal form of B will have 0, , and on its main diagonal, and zeros on the rest of its space:
Then we need to calculate the transition matrix Q:
When . Solving the equation above, we get .
When . Solving the equation above,
we get . When . Solving the equation above,
we get . Therefore, the eigenvectors are , , and .
To get the transition matrix, all we have to do is to put together the three eigenvectors we have got:
The next step is to calculate the inverse matrix of Q. Since the process of calculation is very complicated, here we will directly give the final result:
Then, we can calculate the resolving kernel through using the matrix method for solving ODEs mentioned before:
Having gotten the resolving kernel , the evolution function of the velocity can be expressed as where .
Then with can be found by the equation where .
6. Applications in Curves
6.1. Inner Product in R3
Assume that there are two vectors, , in R3. We define the inner product of u and v to be . The modulus of u is . The calculation of the inner product has the following properties:
2) , then if ,u is then perpendicular to v
3) , where is a constant
Let and be two differential maps of an open interval , we have
6.2. Regular Curve
Definition 5 (Regular curve) A curve is called a regular curve if
In this case, we can find another parameter such that under this parameter, the velocity is constantly 1 for a:
Proof: We try to find such that is of velocity 1.
That is, .
Therefore, . Then, to find , all we need to do is to calculate an integral:
Example 5 Given , find such that B and a are the same curve, but .
First, we calculate the derivative and the modulus of the derivative of :
Then, find :
t therefore equals to and
6.3. Vector Products in R3
Definition 6 (Vector Product) LetV be a three dimensional vector space, and be vectors in such a vector space,the vector product, , is defined tobe the vector in satisfying .
The vector product has the following properties:
3) If , the u and v are parallel.
Definition 7 (Curvature) Let be a regular curve, be the reparametrization of such that ,the curvature of ,k,is defined to be .
Property 1: If is a straight line, the curvature
Property 2: If is a circle of radius R, the curvature
Example 6 Find the curvature of curve :
We have already found the reparameterization of
To find the curvature k, we just have to calculate , which equals to
From now on, we assume that is of arc-length parameter (i.e, ).
Definition 8 The normal vector is defined to be .
Theorem 3 For a given curve of arc-length parameter, is perpendicular to .
Proof: Since . Then
Therefore, , and is perpendicular to .
Definition 9 The binormal vector of curve a at s is defined to be .
The modulus of is constantly 1, since
We now present the definition of torsion, a concept that will be further explored based on our previous introduction concerning solving linear ODEs.
Definition 10 (Torsion) For a regular curve of arc-length parameter, whose binormal vector is , defined to be the torsion ofthe curve.
Proof: can be expressed as
so , and since is perpendicular to , we can then conclude that is perpendicular to . Therefore, , as is also perpendicular to .
Now that assume we are given the curvature as well as the torsion of an arclength curve , we need to calculate the binormal vector, the normal vector, and the velocity vector. Our first step is to find the way to represent , , and using , , as well as and .
As defined, . Since is perpendicular to both and , we can conclude that . As , can be expressed as . Finally, since
Therefore, we have the following equations:
and the differential equation concerning can be written in the following form:
To find , , and , we need to introduce and prove the Fundamental Theorem of The Local Theory Of Curves: Given functions , differentiable, we can find a space curve such that , . Furthermore, if are two curves satisfying this condition, then there exists a rigid transformation A such that .
According to Frenet Formula (Reinhard Schäfke, Dieter Schmidt, 2016) , is clearly an ODE. According to our linear ODE theory, given , there exists a unique solution to the equation . Therefore, we can find a group of satisfying the given
equation. Then, accordingly, we can find a such that , where satisfies , . Then, assume vector is among the vectors satisfying this linear ODE, and vector , where A is a transformation, also satisfies this ODE, then according to the uniqueness of the solution to ODEs, must equal to . Till now, we have proved the existence of the unique solution and the existence and uniqueness of .
Finally, we have the following equations:
Therefore, we proved that if are two curves satisfying this condition, then there exists a rigid transformation A such that .
Based on the six sections above, we see that linear ODEs has a vast range of application in Physics and Geometry. Only through using them can we characterize certain evolution functions and solve equations relating the rate of change of a function and its value. Also, through applying the matrix method for solving linear ODEs, we are able to simplify and solve equations of higher dimension. Overall, the matrix method is quite useful and effective in calculating multiple variables at the same time.
Some physics background: A theoretical expression for the viscous force on a ball at low speed is given by Stoke: , where r represents the radius of the ball, v the speed and the viscosity of the fluid. In the question, we may write directly .
 Vizvari, Z., Sari, Z., Klincsik, M. and Odry, P. (2020) Exact Schemes for Second-Order Linear Differential Equations in Self-Adjoint Cases. Advances in Difference Equations, 2020.
 Zdeněk, H. (2005) Continuous Dependence of Inverse Fundamental Matrices of Generalized Linear Ordinary Differential Equations on a Parameter. Acta Universitatis Palackianae Olomucensis. Facultas Rerum Naturalium. Mathematica, 44, 39-48.
 Schäfke, R. and Schmidt, D. (2006) The Connection Problem for General Linear Ordinary Differential Equations at Two Regular Singular Points with Applications in the Theory of Special Functions. SIAM Journal on Mathematical Analysis, 11, 848-862.