Simplify 3rd order ODE to v'''+Z v = F general formula​
​by Rauan Kaldybaev
A third-order linear ordinary differential equation in general has the form
u'''+au''+bu'+cu=f
, where
u=u[t]
is the unknown function and
a,b,c,f
are known functions of
t
. These equations are typically very difficult to solve explicitly, so instead, people try to at least understand their behavior through simplifications and tricks. One such trick is presented in this essay: how to convert a general 3rd order linear ODE to the (relatively) simple form
v'''+Zv=F
via a change of variables - in other words, how to set the coefficients before
v''
and
v'
to zero. The computation is fairly nontrivial: it involves solving another ODE
w''[t]=
2
p
[t]w[t]
and taking a couple of difficult integrals; however, more often than not it is doable and thus makes a valuable addition to the mathematical toolbox.

Deriving the formula

A general third-order linear ordinary differential equation has the form
u'''[t]+a[t]u''[t]+b[t]u'[t]+c[t]u[t]=f[t]
(
1
)
The equation is currently written in terms of
u[t]
, however, we can make a change of variables and re-write it in terms of some other function
v[t]
. There are many ways in which a change of variables could be done. Consider the following transformation:
u[t]=A[t]v[ϕ[t]]
(
2
)
Effectively, this change of variables replaces
t
with
ϕ
and multiplies the expression by some scaling function
A[t]
. Substituting this into (1) yields
A
3
(ϕ')
3
d
v
3
dϕ
+(aA
2
(ϕ')
+3A'
2
(ϕ')
+3Aϕ'ϕ'')
2
d
v
2
dϕ
​​+(Abϕ'+2aA'ϕ'+3ϕ'A''+aAϕ''+3A'ϕ''+Aϕ''')
dv
dϕ
+(Ac+bA'+aA''+A''')v=f[t]
(
3
)
That is how the ODE is written in terms of the variables
v
and
ϕ
. Effectively, it is in the same form as the initial ODE (1), the only difference being that the coefficients have changed. Here,
d
dϕ
denotes a
ϕ
derivative and
'
is reserved for the time derivative. At this point, the functions
A
and
ϕ
could be anything, but for the change of variables to be useful, let us try to find such
A
and
ϕ
that the coefficients before
2
d
v
2
dϕ
and
dv
dϕ
are zero:
aA
2
k
+3A'
2
k
+3Akk'=0​​Abk+2aA'k+3kA''+aAk'+3A'k'+Ak''=0
(
4
)
Here,
k≡ϕ'[t]
is a short-hand notation for
dϕ/dt
. If we solve these equations, the ODE (3) will have the desired simple form. Dividing the first equation by
A
2
k
yields
a+3
A'
A
+3
k'
k
=0
. Integrating this with respect to time, we get
t
∫
a[
t
1
]
dt
1
+3Log[A]+3Log[k]=
∼

. Solving for
A[t]
yields
A=

1
Exp-
1
3
t
∫
a[
t
1
]
dt
1

1
k
(
5
)
Plugging this in to the second equation in (4), we get
3
2
k'
k
-2
k''
k
=
1
3
2
a
+a'-b
. The equation is nonlinear which makes it difficult to solve; however, it can be made linear by introducing
w[t]≡
1
k[t]
(
6
)
The equation
3
2
k'
k
-2
k''
k
=
1
3
2
a
+a'-b
is written in terms of
w
as
w''[t]=
2
p
[t]w[t]
(
7
)
Where
2
p
[t]≡
2
a[t]
12
+
a'[t]-b
4
. Unlike the previous equation, (7) is linear, which makes it tremendously easier to solve, both exactly and approximately. From
w[t]≡1
k[t]
(6),
k
can then be computed as
k=
1
2
w
. The definition
w[t]≡1
k[t]
requires that
k!=0
. Coming back to the definition
k[t]≡ϕ'[t]
and requiring
ϕ[t]
to be smooth, we see that the condition
k!=0
simply demands a one-to-one relationship between
ϕ
and
t
. This is very natural because
ϕ
in the change of variables (2) is the “transformed time”, and it only makes sense if there is a one-to-one relationship between time and “transformed time”. Substituting this into (5) and setting

1
=1
, we get the final formula for
A[t]
:
A[t]=Exp-
1
3
t
∫
a[
t
1
]
dt
1

2
w
[t]
(
8
)
Making

1
equal to any other nonzero number (for instance, 2) wouldn’t really change much because it would just correspond to rescaling
A
by some constant. Since
k[t]≡ϕ'[t]
,
ϕ[t]
is the integral of
k[t]=1/
2
w
:
ϕ[t]=
t
∫
dt
1
2
w
[
t
1
]
(
9
)
And this is it! The problem is solved. If we substitute
u[t]=A[t]v[ϕ[t]]
, where
A
and
ϕ
are given by (8) and (9), into the ODE (3), we would see that the coefficients before
v''
and
v'
are indeed zero. The ODE can be written in the form
3
d
v
3
dϕ
+Zv=F
(
10
)
Where
Z
and
F
are given by
Z[ϕ]≡
1
54
6
w
[t](4
3
a[t]
-18a[t](b[t]-
′
a
[t])+9(6c[t]-3
′
b
[t]+
′′
a
[t]))​​F[ϕ]=
6
w
[t]
A[t]
f[t]
(
11
)

Code used to do the calculations


Computer implementation

The function below executes the algorithm described above. It first computes
2
p
and tries to find
w[t]
. If no solutions were found, it just returns the general formulas for
ϕ,A,Z,F
. If one solution was found, it substitutes the formula for
w
into (8) and (9) and outputs concrete formulas for
A,ϕ,Z
. If multiple solutions were found, it outputs the general formulas for
ϕ,A,Z,F
and gives the user the list of possible
w
.
In[]:=
ClearAll[findϕAZF];​​findϕAZF[a_Function, b_Function, c_Function, f_Function] := Block​​ ​​ t, p2, w, possiblew, wIndicator, Aval​​ ,​​ p2 =
2
a[t]
12
+
a'[t] - b[t]
4
;​​ possiblew = DSolve[w''[t]  p2 w[t], w, t];​​ ​​ Which[​​ SameQ[Head[possiblew], DSolve], wIndicator = Evaluate[w''[t]  p2 w[t]], (* if couldn't solve *)​​ Length[possiblew]  1, w = possiblew[[1, 1, 2]]; wIndicator = Nothing, (* if found exactly one solution *)​​ Length[possiblew] > 1, wIndicator = (w ∈ possiblew[[All, 1, 2]]) (* if found several solutions *)​​ ];​​ ​​ Aval = Exp-
1
3
Integrate[a[t], t]
2
w[t]
;​​ Simplify@ϕ  Integrate
1
2
w[t]
, t, wIndicator, A  Aval,​​ Z 
1
54
6
w[t]
(4
3
a[t]
-18 a[t] (b[t]-
′
a
[t])+9 (6 c[t]-3
′
b
[t]+
′′
a
[t])),​​ F ->
6
w[t]
Aval
f[t]​​;
In[]:=
ClearAll[getZFform];​​getZFform[{ϕ -> ϕval_, A -> Aval_, Z -> Zval_, F -> Fval_}] := Block[{texpr},​​ texpr = Solve[ϕval == ϕ, t];​​ Which[​​ SameQ[Head[texpr], Solve], Return[$Failed],​​ Length[texpr] >= 1, texpr = texpr[[1, 1, 2]];​​ ​​ u[t] == Aval v[ϕ],​​ t == texpr,​​ Z -> Zval /. t -> texpr,​​ F -> Fval /. t -> texpr​​ ​​ ]​​];
Here is how to use this method:
1
.
Get the ODE to the form
u'''[t]+a[t]u''[t]+b[t]u'[t]+c[t]u[t]=f[t]
; identify the functions
a[t],b[t],c[t],f[t]
.
2
.
Run findϕAZF to find
ϕ[t],A[t],Z[t],F[t]
.
3
.
The expressions for
ϕ[t],A[t],Z[t],F[t]
are going to contain two arbitrary constants. Choose values for the constants, preferably to make the expressions simpler.
4
.
Run getZFform to find
Z
and
F
expressed in terms of
ϕ
.
Simple! Of course, this isn’t going to work with each and every ODE. Often, the program is unable to find
ϕ[t]
, or invert
ϕ[t]
to find
t[ϕ]
, or to compute some of the integrals. In those cases, one can resort to manual tuning and computation, which would often help lift the stalemate. It is also possible to use approximate methods such as Frobenius series. This would require some work, but even then the process isn’t too laborious.

Simple example:
u'''+u''+u'+u=Log[6+
2
t
]

The equation is very simple, which makes it a great initial example to show that the formula works. On practice we could just solve it explicitly without bothering to get it to the
v'''+Zv=F
form, but it is just interesting to see how things would go. Step 1: We can see that the coefficients in this ODE are
a=b=c=1
, and the external force is
f[t]=Log[6+
2
t
]
. Step 2: Run the function findAϕZF:
In[]:=
findϕAZF[1&,1&,1&,Function[t,Log[6+
2
t
]]]
Out[]=
ϕ
6
Sin
t
6

2

1
Cos
t
6
+

1

2
Sin
t
6
,A
-t/3

2

1
Cos
t
6
+

2
Sin
t
6

,Z
20
27
6

1
Cos
t
6
+

2
Sin
t
6

,F
t/3

Log[6+
2
t
]
4

1
Cos
t
6
+

2
Sin
t
6


Step 3: The choice of constants

1
,

2
is arbitrary. Let us take

1
=1,

2
=0
:
In[]:=
%/.{C[1]1,C[2]0}
Out[]=
ϕ
6
Tan
t
6
,A
-t/3

2
Cos
t
6

,Z
20
27
6
Cos
t
6

,F
t/3

4
Cos
t
6

Log[6+
2
t
]
Step 4: Run getZFform:
In[]:=
getZFform[%]
Out[]=
u[t]
-t/3

2
Cos
t
6

v[ϕ],
t
6
ArcTan
ϕ
6
+π

1
if

1
∈
,Z
20
27
6
CosArcTan
ϕ
6
+π

1

if

1
∈
,F
2
3
ArcTan
ϕ
6
+π

1

4
CosArcTan
ϕ
6
+π

1

Log6+6
2
ArcTan
ϕ
6
+π

1
 if

1
∈

The expression contains an arbitrary constant

1
due to the translation symmetry of trigonometric functions;
t[ϕ]
can in general contain arbitrary constants. Take

1
=0
:
In[]:=
%/.C[1]->0
Out[]=
u[t]
-t/3

2
Cos
t
6

v[ϕ],t
6
ArcTan
ϕ
6
,Z
20
27
3
1+
2
ϕ
6
,F
2
3
ArcTan
ϕ
6


Log6+6
2
ArcTan
ϕ
6


2
1+
2
ϕ
6

So, the ODE can be written as
3
d
v
3
dϕ
+
160
3
(6+
2
ϕ
)
v=36
2
3
ArcTan
ϕ
6


Log61+
2
ArcTan
ϕ
6


2
(6+
2
ϕ
)
If we now substitute
u[t]=A[t]v[ϕ]
into the original ODE
u'''[t]+u''[t]+u'[t]+u[t]=Log[6+
2
t
]
, we see that it indeed holds:
In[]:=
Simplifyu'''[t]+u''[t]+u'[t]+u[t]/.u->Function[t,A[t]v[ϕ[t]]]/.ϕFunctiont,
6
Tan
t
6
,AFunctiont,
-t/3

2
Cos
t
6

/.v'''[ϕ_]->-
160
3
(6+
2
ϕ
)
v[ϕ]+36
2
3
ArcTan
ϕ
6


Log61+
2
ArcTan
ϕ
6


2
(6+
2
ϕ
)
,Assumptions->{-π/2<t<π/2}
Out[]=
Log[6+
2
t
]
Notice that the result only makes sense for
-π/2<t<π/2
because
ϕ=
6
Tan
t
6

goes to infinity at
t=±π/2
. This is one of the drawbacks of the formula: it can often introduce limitations out of nowhere.

A more complicated example:
u'''+tu''+u'+u=Sin[t]

The ODE
u'''+tu''+u'+u=Sin[t]
is very difficult to solve exactly, and this is one of those cases when the method presented in this essay may be of interest. The function
ϕ[t]
could not be found explicitly for this ODE, and instead, an indefinite integral is returned:
In[]:=
findϕAZF[#&,1&,1&,Function[t,Sin[t]]]
Out[]=
ϕ1
2

2
ParabolicCylinderD-
1
2
,
t
1/4
3
+

1
ParabolicCylinderD-
1
2
,
t
1/4
3

t,A
-
2
t
6

2

2
ParabolicCylinderD-
1
2
,
t
1/4
3
+

1
ParabolicCylinderD-
1
2
,
t
1/4
3

,Z
1
27
(27+2
3
t
)
6

2
ParabolicCylinderD-
1
2
,
t
1/4
3
+

1
ParabolicCylinderD-
1
2
,
t
1/4
3

,F
2
t
6

4

2
ParabolicCylinderD-
1
2
,
t
1/4
3
+

1
ParabolicCylinderD-
1
2
,
t
1/4
3

Sin[t]
The integral can be evaluated approximately using methods such as Taylor series or numerically; for simplicity, the example will go with the latter option. Choose

1
=1
,

2
=0
:
In[]:=
%/.{C[1]->1,C[2]->0}
Out[]=
ϕ1
2
ParabolicCylinderD-
1
2
,
t
1/4
3

t,A
-
2
t
6

2
ParabolicCylinderD-
1
2
,
t
1/4
3

,Z
1
27
(27+2
3
t
)
6
ParabolicCylinderD-
1
2
,
t
1/4
3

,F
2
t
6

4
ParabolicCylinderD-
1
2
,
t
1/4
3

Sin[t]
Now, the equation can be written as
3
d
v
3
dϕ
+
1
27
(27+2
3
t
)
6
ParabolicCylinderD-
1
2
,
t[ϕ]
1/4
3

v=
2
t
6

4
ParabolicCylinderD-
1
2
,
t
1/4
3

Sin[t]
The function
ϕ[t]
could not be found, not to mention inverted to find
t[ϕ]
. However, it’s easy to compute a numerical solution. If
ϕ=
1
2
ParabolicCylinderD-
1
2
,
t
1/4
3

t
, then
dϕ
dt
=
1
2
ParabolicCylinderD-
1
2
,
t
1/4
3

, and
dt
dϕ
=
2
ParabolicCylinderD-
1
2
,
t
1/4
3

. This is a simple differential equation that can be solved numerically, or, if more accuracy is needed, using a series solution. For simplicity, this example will go with the former option:
In[]:=
myt=NDSolvet'[ϕ]==
2
ParabolicCylinderD-
1
2
,
t[ϕ]
1/4
3

,t[0]==0,t,{ϕ,0,3}[[1,1,2]];
In[]:=
Plot[myt[ϕ],{ϕ,0,3},AxesLabel->{"ϕ","t"}]
Out[]=
To find the functions
Z
and
F
, we just substitute the numerical result into the formulas:
In[]:=
Plot[#[[2]],{ϕ,0,3},AxesLabel->{"ϕ",#[[1]]},ImageSize->Medium]&/@"Z",
1
27
(27+2
3
t
)
6
ParabolicCylinderD-
1
2
,
t
1/4
3

,"F",
2
t
6

4
ParabolicCylinderD-
1
2
,
t
1/4
3

Sin[t]/.t->myt[ϕ]
Out[]=

,

There is no closed formula for
Z
and
F
, however, the functions can be computed to any specified accuracy. Alternatively, if a Taylor series solution for
t[ϕ]
is found, it can also be substituted into
Z
and
F
to obtain a series expansion for those. The important thing is that the ODE can be converted to the
v'''+Zv=F
form, which paves the way for many different ways of analyzing it further. As you can see, the simplification method also works on the more difficult ODEs.