Chebyshev Polynomial-Based Analytic Solution Algorithm with Efficiency, Stability and Sensitivity for Classic Vibrational Constant Coefficient Homogeneous IVPs with Derivative Orders n, n-1, n-2

David P. Stapleton

Department of Mathematics and Statistics, University of Central Oklahoma, Edmond, OK, USA.

DOI: 10.4236/ajcm.2022.124023                

Abstract

The Chebyshev polynomials are harnessed as functions of the one parameter of the nondimensionalized differential equation for trinomial homogeneous linear differential equations of arbitrary order n that have constant coefficients and exhibit vibration. The use of the Chebyshev polynomials allows calculation of the analytic solutions for arbitrary n in terms of the orthogonal Chebyshev polynomials to provide a more stable solution form and natural sensitivity analysis in terms of one parameter and the initial conditions in 6n + 7 arithmetic operations and one square root.

Keywords

Differential Equation, Stability, Sensitivity Analysis, Chebyshev Polynomials, Coefficient Formula

1. Introduction

The basic homogeneous linear IVP (Initial Value Problem) with constant coefficients and negative discriminant

$\begin{array}{l}a\frac{{\text{d}}^{n}y}{\text{d}{t}^n}+b\frac{{\text{d}}^{n-1}y}{\text{d}{t}^{n-1}}+c\frac{{\text{d}}^{n-2}y}{\text{d}{t}^{n-2}}=0\left(b^2-4ac<0\right);\\ y\left(0\right)=y_o,y^{\prime }\left(0\right)={y^{\prime }}_o,\cdots ,y^{\left(n-1\right)}\left(0\right)=y_o^{\left(n-1\right)}\end{array}$ (1)

is easily solved by elementary methods, but when the problem is of a high order, calculation of the coefficients and presentation of the solution through common means is prone to both computational and presentation instabilities. Some common applications involving higher order equations appear in control systems theory [1], the theory of vibrating beams [2], studies of consistency and stability [3], and theory of the behavior of electronic networks [4]. This paper addresses these issues for solving higher order problems in the special case of Equation (1) by using the simpler nondimensionalized form of the equation and efficiently giving its solution in terms of orthogonal polynomials.

2. The Classic Problem

Upon division of (1) by a and the change of coefficients

$T=\sqrt{\frac{a}{c}}$, $\xi =-\frac{b}{2\sqrt{ac}}$ (which is b/a = −2ξ/T, c/a = 1/T²) (2)

the problem becomes

$\begin{array}{l}\frac{{\text{d}}^{n}y}{\text{d}{t}^n}-\frac{2\xi }{T}\frac{{\text{d}}^{n-1}y}{\text{d}{t}^{n-1}}+\frac{1}{T^2}\frac{{\text{d}}^{n-2}y}{\text{d}{t}^{n-2}}=0;\\ y\left(0\right)=y_o,y^{\prime }\left(0\right)={y^{\prime }}_o,\cdots ,y^{\left(n-1\right)}\left(0\right)=y_o^{\left(n-1\right)}\end{array}$ (3)

with $-\text{1}<\xi <\text{1}$. The substitution $t=T\tau$, $\psi \left(\tau \right)=y\left(T\tau \right)=y\left(t\right)$, so that

$\frac{\text{d}y}{\text{d}t}=\frac{1}{T}\frac{\text{d}\psi }{\text{d}\tau },\cdots ,\frac{{\text{d}}^{n}y}{\text{d}{t}^n}=\frac{1}{T^n}\frac{{\text{d}}^n\psi }{\text{d}{\tau }^n}$ (4)

provides initial values ${\psi }_o^{\left(0\right)}=y_o$, ${\psi }_o^{\left(1\right)}=T{y^{\prime }}_o$, $\cdots$, ${\psi }_o^{\left(n-1\right)}=T^{n-1}{y}_o^{\left(n-1\right)}$ and the equivalent nondimensionalized IVP¹ for ( $-\text{1}<\xi <\text{1}$ )

$\begin{array}{l}\frac{{\text{d}}^n\psi }{\text{d}{\tau }^n}-2\xi \frac{{\text{d}}^{n-1}\psi }{\text{d}{\tau }^{n-1}}+\frac{{\text{d}}^{n-2}\psi }{\text{d}{\tau }^{n-2}}=0;\\ \psi \left(0\right)=y_o,{\psi }^{\prime }\left(0\right)=T{y^{\prime }}_o,\cdots ,{\psi }^{\left(n-1\right)}\left(0\right)=T^{n-1}{y}_o^{\left(n-1\right)}\end{array}$ (5)

It is this form of the problem that suggests the utility of a Chebyshev polynomial approach to solutions.

3. Solution in Terms of Chebyshev Polynomials of the Second Kind

The Laplace transform of (5), in terms of the Laplace transform $\Psi \left(s\right)=\mathcal{L}\left(\psi \left(t\right)\right)$, is

$\begin{array}{l}{s}^n\Psi \left(s\right)-\left(s^{n-1}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-1\right)}\right)-2\xi \left(s^{n-1}\Psi \left(s\right)-\left(s^{n-2}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-2\right)}\right)\right)\\ +\left(s^{n-2}\Psi \left(s\right)-\left(s^{n-3}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-3\right)}\right)\right)=0\end{array}$ (6)

and can be algebraically solved for Ψ(s) to get

$\Psi \left(s\right)=\frac{\left(s^{n-1}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-1\right)}\right)-2\xi \left(s^{n-2}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-2\right)}\right)+\left(s^{n-3}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-3\right)}\right)}{s^n+2\xi s^{n-1}+s^{n-2}}$ (7)

or in terms of $\omega =\sqrt{1-{\xi }^2}$

$\Psi \left(s\right)=\frac{\left(s^{n-1}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-1\right)}\right)-2\xi \left(s^{n-2}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-2\right)}\right)+\left(s^{n-3}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-3\right)}\right)}{\left({\left(s+\xi \right)}^2+{\omega }^2\right)s^{n-2}}$ (8)

By the method of partial fractions, there exist constants $c_0,\cdots ,c_{n-1}$ such that

$\Psi \left(s\right)=\frac{\left(c_0-2\xi c_1\right)\left(s-\xi \right)+\left(\xi c_0+\left(1-2{\xi }^2\right)c_1\right)}{{\left(s-\xi \right)}^2+{\omega }^2}+\frac{c_2}{s}+\cdots +\frac{c_{n-1}}{s^{n-2}}$ (9)

and the Laplace transform formulas

$\begin{array}{l}{\mathcal{L}}^{-1}\left(\frac{\omega }{{\left(s-a\right)}^2+{\omega }^2}\right)={\text{e}}^{a\tau }\sin \left(\omega \tau \right),\\ {\mathcal{L}}^{-1}\left(\frac{s-a}{{\left(s-a\right)}^2+{\omega }^2}\right)={\text{e}}^{a\tau }\cos \left(\omega \tau \right),\\ {\mathcal{L}}^{-1}\left(\frac{1}{s^n}\right)=\frac{{\tau }^{n-1}}{\left(n-1\right)!}\end{array}$ (10)

applied to (9) show that the form of the general solution of the DE in (5) is expressible as

$\boxed{\psi =\left(c_0-2\xi c_1\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)+\frac{\xi c_0+\left(1-2{\xi }^2\right)c_1}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)+{\displaystyle \underset{k=2}{\overset{n-1}{\sum }}\frac{c_k}{\left(k-2\right)!}{\tau }^{k-2}}}$ (11)

This leaves the problem of determining the constants. By (8) and (9) they satisfy

$\begin{array}{l}\frac{\left(c_0-2\xi c_1\right)\left(s-\xi \right)+\left(\xi c_0+\left(1-2{\xi }^2\right)c_1\right)}{{\left(s-\xi \right)}^2+{\omega }^2}+\frac{c_2}{s}+\cdots +\frac{c_{n-1}}{s^{n-2}}\\ =\frac{\left(s^{n-1}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-1\right)}\right)-2\xi \left(s^{n-2}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-2\right)}\right)+\left(s^{n-3}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-3\right)}\right)}{\left({\left(s-\xi \right)}^2+{\omega }^2\right)s^{n-2}}\end{array}$ (12)

and multiplication by $\left({\left(s-\xi \right)}^2+{\omega }^2\right)s^{n-2}$ shows that

$\begin{array}{l}\left(\left(c_0-2\xi c_1\right)\left(s-\xi \right)+\left(\xi c_0+\left(1-2{\xi }^2\right)c_1\right)\right)s^{n-2}+c_2\left({\left(s-\xi \right)}^2+{\omega }^2\right)s^{n-3}\\ +\cdots +c_{n-1}\left({\left(s-\xi \right)}^2+{\omega }^2\right)\\ =\left(s^{n-1}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-1\right)}\right)-2\xi \left(s^{n-2}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-2\right)}\right)+\left(s^{n-3}{\psi }_o^{\left(0\right)}+\cdots +{\psi }_o^{\left(n-3\right)}\right)\end{array}$ (13)

Since this holds for all sufficiently large Re(s) the coefficients of powers of s match as

$c_0-2\xi c_1+c_2={\psi }_o^{(\; 0\; )}$

$c_1-2\xi c_2+c_3={\psi }_o^{\left(1\right)}-2\xi {\psi }_o$

$\left({\xi }^2+{\omega }^2\right)c_k-2\xi c_{k+1}+c_{k+2}={\psi }_o^{\left(k+1\right)}-2\xi {\psi }_o^{\left(k\right)}+{\psi }_o^{\left(k-1\right)}$ ( $2\le k\le n-4$ ) (14)

$-2\xi c_{n-2}+\left({\xi }^2+{\omega }^2\right)c_{n-1}={\psi }_o^{\left(n-2\right)}-2\xi {\psi }_o^{\left(n-3\right)}+{\psi }_o^{\left(n-4\right)}$

$\left({\xi }^2+{\omega }^2\right)c_{n-1}={\psi }_o^{\left(n-1\right)}-2\xi {\psi }_o^{\left(n-2\right)}+{\psi }_o^{\left(n-3\right)}$

Also, ${\xi }^2+{\omega }^2=1$, so in matrix form this is a system $U^{-\text{1}}c={\left(U^{-\text{1}}\right)}^{\text{T}}{\psi }_o$ in terms of

$c={\left[\begin{array}{cccc}{c}_0& c_1& \cdots & c_{n-1}\end{array}\right]}^{\text{T}}$, ${\psi }_o={\left[\begin{array}{cccc}{\psi }_o^{\left(0\right)}& {\psi }_o^{\left(1\right)}& \cdots & {\psi }_o^{\left(n-1\right)}\end{array}\right]}^{\text{T}}$ (15)

in which U⁻¹ is an n × n, upper triangular Toeplitz matrix with unit diagonal.

$\left[\begin{array}{cccccc}1& -2\xi & 1& 0& \cdots & 0\\ 0& 1& -2\xi & 1& \ddots & 0\\ 0& 0& 1& -2\xi & \ddots & 0\\ 0& 0& 0& 1& \ddots & 1\\ ⋮& & \ddots & \ddots & \ddots & -2\xi \\ 0& 0& \cdots & 0& 0& 1\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ ⋮\\ c_{n-2}\\ c_{n-1}\end{array}\right]=\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ -2\xi & 1& 0& 0& \cdots & 0\\ 1& -2\xi & 1& 0& \cdots & 0\\ 0& 1& -2\xi & 1& \ddots & ⋮\\ & \ddots & \ddots & \ddots & \ddots & 0\\ 0& 0& 0& 1& -2\xi & 1\end{array}\right]\left[\begin{array}{c}{\psi }_o^{\left(0\right)}\\ {\psi }_o^{\left(1\right)}\\ {\psi }_o^{\left(2\right)}\\ {\psi }_o^{\left(3\right)}\\ ⋮\\ {\psi }_o^{\left(n-1\right)}\end{array}\right]$

To solve this, consider that the values U₀(ξ) = 1, U₁(ξ) = 2ξ of the Chebyshev polynomials of the second kind [5] U_k(ξ) ( $k=0,\text{1},\cdots$ ) imply that

$-2\xi U_0\left(\xi \right)+U_1\left(\xi \right)=1$ (16)

and the recurrence relation $U_{k+1}\left(\xi \right)=2\xi U_k\left(\xi \right)-U_{k-1}\left(\xi \right)$ ( $k=\text{1},\text{2},\cdots$ ) that defines the polynomials for larger values of k is also

$U_{k+1}\left(\xi \right)-2\xi U_k\left(\xi \right)+U_{k-1}\left(\xi \right)=0\left(k=1,2,\cdots \right)$ (17)

These relations show (checked by matrix multiplication) that the inverse of U⁻¹ is

$U=\left[\begin{array}{ccccc}{U}_0\left(\xi \right)& U_1\left(\xi \right)& U_2\left(\xi \right)& \cdots & U_{n-1}\left(\xi \right)\\ 0& U_0\left(\xi \right)& U_1\left(\xi \right)& \cdots & U_{n-2}\left(\xi \right)\\ 0& 0& U_0\left(\xi \right)& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & U_1\left(\xi \right)\\ 0& \cdots & 0& 0& U_0\left(\xi \right)\end{array}\right]$ (18)

The coefficients $c=U{\left(U^{-\text{1}}\right)}^{\text{T}}{\psi }_o$ in (11) are therefore

$\left[\begin{array}{c}{c}_0\\ c_1\\ ⋮\\ c_{n-2}\\ c_{n-1}\end{array}\right]=\left[\begin{array}{ccccc}{U}_0\left(\xi \right)& U_1\left(\xi \right)& U_2\left(\xi \right)& \cdots & U_{n-1}\left(\xi \right)\\ 0& U_0\left(\xi \right)& U_1\left(\xi \right)& \cdots & U_{n-2}\left(\xi \right)\\ 0& 0& U_0\left(\xi \right)& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & U_1\left(\xi \right)\\ 0& \cdots & 0& 0& U_0\left(\xi \right)\end{array}\right]\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ -2\xi & 1& 0& 0& \cdots & 0\\ 1& -2\xi & 1& 0& \cdots & 0\\ 0& 1& -2\xi & 1& \ddots & ⋮\\ & \ddots & \ddots & \ddots & \ddots & 0\\ 0& 0& 0& 1& -2\xi & 1\end{array}\right]\left[\begin{array}{c}{\psi }_o^{\left(0\right)}\\ {\psi }_o^{\left(1\right)}\\ {\psi }_o^{\left(2\right)}\\ {\psi }_o^{\left(3\right)}\\ ⋮\\ {\psi }_o^{\left(n-1\right)}\end{array}\right]$ (19)

Let A = U(U⁻¹)^T denote the n × n matrix product in this expression (so that c = Aψ_o). Then

A has the following structure. The zeros in the matrices force a_jk = 0 for j ≥ k + 3. By multiplication a_k+2,k = U_o(ξ) = 1 for 1 ≤ k ≤ n − 2, and a_k+1,k = −2ξU_o(ξ) + U₁(ξ). = 0 (by (16)) for 1 ≤ k ≤ n − 1. By multiplication a_jk = U_k+1(ξ) − 2ξU_k(ξ) + U_k−1(ξ) = 0 (by (17)) for j ≤ k ≤ n − 2. This leaves only the last two columns of A. By multiplication and (17) column n − 1 contains a_j,n−1 = U_n−j−1(ξ) − 2ξU_n−j(ξ) = −U_n−j+1(ξ) for j ≤ n − 1 and a_n,n−1 = −2ξU₀(ξ) = −2ξ = −U₁(ξ). By multiplication column n contains a_j,n = U_n−j(ξ). Hence A is the n × n matrix

$A=\left[\begin{array}{cccccc}0& 0& \cdots & 0& -U_n\left(\xi \right)& U_{n-1}\left(\xi \right)\\ 0& 0& \cdots & 0& -U_{n-1}\left(\xi \right)& U_{n-2}\left(\xi \right)\\ 1& 0& \ddots & ⋮& ⋮& ⋮\\ 0& 1& \ddots & 0& -U_3\left(\xi \right)& U_2\left(\xi \right)\\ ⋮& \ddots & \ddots & 0& -U_2\left(\xi \right)& U_1\left(\xi \right)\\ 0& \cdots & 0& 1& -U_1\left(\xi \right)& U_0\left(\xi \right)\end{array}\right]$ (20)

Upon defining ${\psi }_o^{\left(-1\right)}=0$, ${\psi }_o^{\left(-2\right)}=0$, the coefficients c = Aψ_o can therefore be written as²

$c_k={\psi }_o^{\left(k-2\right)}-U_{n-k}\left(\xi \right){\psi }_o^{\left(n-2\right)}+U_{n-k-1}\left(\xi \right){\psi }_o^{\left(n-1\right)}$

or

$\boxed{c_k={\psi }_o^{\left(k-2\right)}+U_{n-k-1}\left(\xi \right){\psi }_o^{\left(n-1\right)}-U_{n-k}\left(\xi \right){\psi }_o^{\left(n-2\right)}\left(k=0,\cdots ,n-1\right)}$ (21)

Theorem: The unique solution of IVP (5) for n ≥ 2 with $\omega =\sqrt{1-{\xi }^2}$, ${\psi }_o^{\left(-1\right)}=0$, ${\psi }_o^{\left(-2\right)}=0$ defined can be expressed as

$\psi =\left(c_0-2\xi c_1\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)+\frac{\xi c_0+\left(1-2{\xi }^2\right)c_1}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)+{\displaystyle \underset{k=2}{\overset{n-1}{\sum }}{c}_k\frac{{\tau }^{k-2}}{\left(k-2\right)!}}$ (22)

in which the coefficients c_k are given by (21) and U_k(ξ) denotes the kth Chebyshev polynomial of the second kind.³

For given values of ξ and initial conditions, the coefficients in the analytic solution (9) of (5) can be calculated in O(n) time as follows.⁴ Computation of $\omega =\sqrt{1-{\xi }^2}$ uses two arithmetic operations and one square root. The Chebyshev polynomials of zeroth through nth degree can be evaluated at the given value of ξ with one multiplication to produce U₁(ξ) and two operations each to compute the subsequent U_k(ξ) ( $k=\text{2},\text{3},\cdots ,n$ ) by means of the recurrence relation (17) (with 2ξ available), for a total of 2n − 1 flops. The coefficients $c_0,\cdots ,c_{n-1}$ in (9) may each be computed by the four flops appearing in (8) for 4n total flops, or in 4n − 2 flops by omitting adds of ${\psi }_o^{\left(-1\right)}=0$, ${\psi }_o^{\left(-2\right)}=0$. The coefficients of ${\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)$ and ${\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)$ can then be completed in 8 flops since 2ξ, c₀, c₁ have already been evaluated. Hence, $\omega =\sqrt{1-{\xi }^2}$ and the coefficients in expression (9) for given ξ and given initial conditions can be evaluated with one square root and 2 + (2n − 1) + (4n − 2) + 8 = 6n + 7 arithmetic operations.

Corollary 1: The value of $\omega =\sqrt{1-{\xi }^2}$ and the coefficients of the terms $\frac{{\tau }^{k-2}}{\left(k-2\right)!}$ ( $k=\text{2},\cdots ,n-1$ ), ${\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)$ and ${\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)$ in Formula (22) for the analytic solution (5) of (2) can be computed with 6n + 7 arithmetic operations and one square root.⁵

Elementary example: By Formula (22) find the unique solution of

$\frac{{\text{d}}^{3}y}{\text{d}{t}^3}+6\frac{{\text{d}}^{2}y}{\text{d}{t}^2}+25\frac{\text{d}y}{\text{d}t}=0;y\left(0\right)=3,y^{\prime }\left(0\right)=-3,y^{″}\left(0\right)=-7$ (23)

Use

$T=\sqrt{\frac{a}{c}}=\sqrt{\frac{1}{25}}=\frac{1}{5}$, $\xi =-\frac{b}{2\sqrt{ac}}=-\frac{6}{2\sqrt{25}}=-\frac{3}{5}$ (24)

with t = Tτ = τ/5 and the initial values for ψ(τ) = y(Tτ)

${\psi }_o^{\left(0\right)}=y_o=3$, ${\psi }_o^{\left(1\right)}=T{y^{\prime }}_o=-3/5$, ${\psi }_o^{\left(2\right)}=T^2{y^{″}}_o=-7/25$ (25)

Then, by (5) the given IVP has equivalent form

$\frac{{\text{d}}^3\psi }{\text{d}{\tau }^3}-2\left(-\frac{3}{5}\right)\frac{{\text{d}}^2\psi }{\text{d}{\tau }^2}+\frac{\text{d}\psi }{\text{d}\tau }=0;\psi \left(0\right)=3,{\psi }^{\prime }\left(0\right)=-3/5,{\psi }^{″}\left(0\right)=-7/25$ (26)

whose solution coefficients and ω are computable with one square root and 6 × 3 + 7 = 25 arithmetic operations. The Chebyshev values with 2ξ = −6/5 known are

$\begin{array}{l}{U}_0\left(-3/5\right)=1\\ U_1\left(-3/5\right)=-6/5\\ U_2\left(-3/5\right)=\left(-6/5\right)\cdot U_1\left(-3/5\right)-U_0\left(-3/5\right)=11/25\\ U_3\left(-3/5\right)=\left(-6/5\right)\cdot U_2\left(-3/5\right)-U_1\left(-3/5\right)=84/125\end{array}$ (27)

By (8), the numerators of the coefficients of the sine, cosine and powers of τ are

$\begin{array}{l}{c}_0={\psi }_o^{\left(-2\right)}+U_2\left(-\frac{3}{5}\right){\psi }_o^{\left(2\right)}-U_3\left(-\frac{3}{5}\right){\psi }_o^{\left(1\right)}=0+\left(\frac{11}{25}\right)\left(-\frac{7}{25}\right)-\left(\frac{84}{125}\right)\left(-\frac{3}{5}\right)=\frac{7}{25}\\ c_1={\psi }_o^{\left(-1\right)}+U_1\left(-\frac{3}{5}\right){\psi }_o^{\left(2\right)}-U_2\left(-\frac{3}{5}\right){\psi }_o^{\left(1\right)}=0+\left(-\frac{6}{5}\right)\left(-\frac{7}{25}\right)-\left(\frac{11}{25}\right)\left(-\frac{3}{5}\right)=\frac{3}{5}\\ c_2={\psi }_o^{\left(0\right)}+U_0\left(-\frac{3}{5}\right){\psi }_o^{\left(2\right)}-U_1\left(-\frac{3}{5}\right){\psi }_o^{\left(1\right)}=3+1\left(-\frac{7}{25}\right)-\left(-\frac{6}{5}\right)\left(-\frac{3}{5}\right)=2\end{array}$ (28)

so in (11) with $\omega =\sqrt{1-{\xi }^2}=\sqrt{1-{\left(-3/5\right)}^2}=4/5$

$\begin{array}{c}\psi =\left(\frac{7}{25}-2\left(-\frac{3}{5}\right)\left(\frac{3}{5}\right)\right){\text{e}}^{-3\tau /5}\cos \left(4\tau /5\right)\\ +\frac{-\frac{3}{5}\frac{7}{25}+\left(1-2{\left(-\frac{3}{5}\right)}^2\right)\frac{3}{5}}{\frac{4}{5}}{\text{e}}^{-3\tau /5}\sin \left(4\tau /5\right)+\frac{2}{0!}\\ ={\text{e}}^{-3\tau /5}\cos \left(4\tau /5\right)+2\end{array}$ (29)

Since τ = t/T = 5t the solution of the original IVP is

$y=\psi \left(5t\right)={\text{e}}^{-3t}\cos \left(4t\right)+2$ (30)

Programming example: The parameters and coefficients of the analytic solution to (1) can be obtained in decimal form by the simplistic MATLAB/Octave function (not optimized and not set up to avoid underflow or overflow) below.

Execution of this program with the above IVP example uses the command

[T,xi,omega,coefs] = IVP_Cheby(1,6,25,[3,−3,−7])

And the output is

T = 0.20000

xi = −0.60000

omega = 0.80000

coefs =

1.00000 0.00000 2.00000

4. Sensitivity Analysis

The derivative of each coefficient (21) with respect to an initial value ${\psi }_o^{\left(j\right)}$ in terms of the Kronecker delta function ${\delta }_k^j$ (augmented with ${\delta }_k^j=0$ for j < 0 or k < 0) by inspection is

$\frac{\partial c_k}{\partial {\psi }_o^{\left(j\right)}}=\{\begin{array}{l}{\delta }_{k-2}^j\left(j<n-2\right)\\ -U_{n-k}\left(\xi \right)\left(j=n-2\right)\\ U_{n-k-1}\left(\xi \right)\left(j=n-1\right)\end{array}$ (31)

so the derivative of the solution (22) of (5) with respect to an initial value ${\psi }_o^{\left(j\right)}$ is⁶

$\frac{\partial \psi }{\partial {\psi }_o^{\left(j\right)}}=\{\begin{array}{l}\left({\displaystyle \underset{k=2}{\overset{n-1}{\sum }}{\delta }_{k-2}^j\frac{{\tau }^{k-2}}{\left(k-2\right)!}}\right)\left(j<n-2\right)\\ \left(\frac{\partial c_0}{\partial {\psi }_o^{\left(n-2\right)}}-2\xi \frac{\partial c_1}{\partial {\psi }_o^{\left(n-2\right)}}\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)\\ +\frac{\xi \frac{\partial c_0}{\partial {\psi }_o^{\left(n-2\right)}}+\left(1-2{\xi }^2\right)\frac{\partial c_1}{\partial {\psi }_o^{\left(n-2\right)}}}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-2\right)\\ \left(\frac{\partial c_0}{\partial {\psi }_o^{\left(n-1\right)}}-2\xi \frac{\partial c_1}{\partial {\psi }_o^{\left(n-1\right)}}\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)\\ +\frac{\xi \frac{\partial c_0}{\partial {\psi }_o^{\left(n-1\right)}}+\left(1-2{\xi }^2\right)\frac{\partial c_1}{\partial {\psi }_o^{\left(n-1\right)}}}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-1\right)\end{array}$ (32)

or

$\frac{\partial \psi }{\partial {\psi }_o^{\left(j\right)}}=\{\begin{array}{l}\frac{{\tau }^j}{j!}\left(j<n-2\right)\\ \left(-U_n\left(\xi \right)+2\xi U_{n-1}\left(\xi \right)\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)-\frac{\xi U_n\left(\xi \right)+\left(1-2{\xi }^2\right)U_{n-1}\left(\xi \right)}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-2\right)\\ \left(U_{n-1}\left(\xi \right)-2\xi U_{n-2}\left(\xi \right)\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)+\frac{\xi U_{n-1}\left(\xi \right)+\left(1-2{\xi }^2\right)U_{n-2}\left(\xi \right)}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-1\right)\end{array}$ (33)

or, by (17) (that is, by $-U_n\left(\xi \right)+2\xi U_{n-1}\left(\xi \right)=U_{n-2}\left(\xi \right)$ and $\xi \left(U_n\left(\xi \right)-2\xi U_{n-1}\left(\xi \right)\right)=-\xi U_{n-2}(\; \xi \; ))$

$\frac{\partial \psi }{\partial {\psi }_o^{\left(j\right)}}=\{\begin{array}{l}\frac{{\tau }^j}{j!}\left(j<n-2\right)\\ U_{n-2}\left(\xi \right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)-\frac{U_{n-1}\left(\xi \right)-\xi U_{n-2}\left(\xi \right)}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-2\right)\\ -U_{n-3}\left(\xi \right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)+\frac{U_{n-2}\left(\xi \right)-\xi U_{n-3}\left(\xi \right)}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-1\right)\end{array}$ (34)

From the Chebyshev relation $U_k\left(\xi \right)-\xi U_{k-1}\left(\xi \right)=T_k\left(\xi \right)$, in which T_k(ξ) denotes the kth Chebyshev polynomial of the first kind, therefore⁷

$\boxed{\frac{\partial \psi }{\partial {\psi }_o^{\left(j\right)}}=\{\begin{array}{l}\frac{{\tau }^j}{j!}\left(j<n-2\right)\\ U_{n-2}\left(\xi \right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)-\frac{T_{n-1}\left(\xi \right)}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-2\right)\\ -U_{n-3}\left(\xi \right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)+\frac{T_{n-2}\left(\xi \right)}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)\left(j=n-1\right)\end{array}}$ (35)

Equations (35) yield the partial derivatives of ψ with respect to ${\psi }_o^{\left(n-2\right)}$ and ${\psi }_o^{\left(n-1\right)}$ (in particular) at τ = 0 in terms of Chebyshev polynomials.

Corollary 2: The partial derivatives of the solution ψ(τ) with respect to initial values ${\psi }_o^{\left(n-2\right)}$, ${\psi }_o^{\left(n-1\right)}$, at τ = 0 are

$\frac{\partial \psi \left(0\right)}{\partial {\psi }_o^{\left(n-1\right)}}=U_{n-2}\left(\xi \right)$, $\frac{\partial \psi \left(0\right)}{\partial {\psi }_o^{\left(n-1\right)}}=-U_{n-3}\left(\xi \right)$ (36)

Elementary example: Find the initial partial derivatives of the solution y(t) to

$\frac{{\text{d}}^{3}y}{\text{d}{t}^3}+6\frac{{\text{d}}^{2}y}{\text{d}{t}^2}+25\frac{\text{d}y}{\text{d}t}=0;y\left(0\right)=3,y^{\prime }\left(0\right)=-3,y^{″}\left(0\right)=-7$ (37)

(the last example) with respect to each of the initial values $y\left(0\right)=3$, $y^{\prime }\left(0\right)=-3$, $y^{″}\left(0\right)=-7$.

By (35) and (36) the partial derivatives of the solution y(t) = ψ(τ) at t = 0 with respect to the initial conditions are

$\begin{array}{l}\frac{\partial y\left(0\right)}{\partial y_o}=\frac{\partial \psi \left(0\right)}{\partial {\psi }_o^{\left(0\right)}}\cdot \frac{\text{d}{\psi }_o^{\left(0\right)}}{\text{d}{y}_o}=\frac{1}{0!}\cdot 1=1\\ \frac{\partial y\left(0\right)}{\partial {y^{\prime }}_o}=\frac{\partial \psi \left(0\right)}{\partial {\psi }_o^{\left(1\right)}}\cdot \frac{\text{d}{\psi }_o^{\left(1\right)}}{\text{d}{y^{\prime }}_o}=U_1\left(\xi \right)\cdot T=2\xi \cdot T=2\cdot \left(-\frac{3}{5}\right)\cdot \frac{1}{5}=-\frac{6}{25}\\ \frac{\partial y\left(0\right)}{\partial {y^{″}}_o}=\frac{\partial \psi \left(0\right)}{\partial {\psi }_o^{\left(2\right)}}\cdot \frac{\text{d}{\psi }_o^{\left(2\right)}}{\text{d}{y^{″}}_o}=-U_0\left(\xi \right)\cdot T^2=-1\cdot T^2=-1\cdot {\left(\frac{1}{5}\right)}^2=-\frac{1}{25}\end{array}$ (38)

5. Summary

The unique analytic solution to a nondimensionalized nth order homogeneous IVP (5) whose differential equation has derivatives of orders n, n − 1, n − 2 only and constant coefficient −1 < ξ < 1 has been expressed using Chebyshev polynomials of the second kind in 6n + 7 arithmetic operations and one square root. This formulation of the solution provides natural Chebyshev polynomial formulas for the partial derivatives of the solution with respect to initial values. By substitution, the solution and its sensitivity analysis can be applied to the common problem of solving trinomial homogeneous linear differential equations with constant coefficients involving derivatives of orders n, n − 1, n − 2, negative discriminant and n initial values.

NOTES

¹Actually y(t) must be divided by its units to make ψ(τ) dimensionless, but this scaling can be ignored since such scaling before solving the IVP only leads to corresponding rescaling after solving. IVP (2) is the n − 2 derivative of the universal oscillator problem when ξ is replaced by −ξ. It will be observed in (5) that the differential Equation (2) has the universal oscillator general solution plus an arbitrary polynomial of degree n − 2 and with ξ replaced by −ξ.

²The reader should recall that the elements c_k, ψ_k of vectors c, ψ are indexed with k running from zero to n − 1.

³In terms of Chebyshev polynomials of the first and second kinds T_k(ξ), U_k(ξ), Equation (9) is also: $\psi =\left(U_0\left(\xi \right)c_0-U_1\left(\xi \right)c_1\right){\text{e}}^{\xi \tau }\cos \left(\omega \tau \right)+\frac{T_0\left(\xi \right)c_0-T_1\left(\xi \right)c_1}{\omega }{\text{e}}^{\xi \tau }\sin \left(\omega \tau \right)+{\displaystyle \underset{k=2}{\overset{n-1}{\sum }}\frac{c_k}{\left(k-2\right)!}{\tau }^{k-2}}$.

⁴The operations count for symbolic representation in terms of arbitrary ξ is different.

⁵The factorials are not in this calculation since for high order problems their inclusion makes underflow likely.

⁶The derivative of (9) with respect to parameter ξ is also possible, but is not presented here.

⁷The derivative ψ’(τ) of (9) has similarly derived partial derivatives that are not presented here.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hendricks, E., Jannerup, O. and Hasse Sorenson, P. (2008) Linear Systems Control: Deterministic and Stochastic Methods. Springer-Verlag, Berlin.
[2]	Han, S., Benaroya, H. and Wei, T. (1999) Dynamics of Transversely Vibrating Beams Using Four Engineering Theories. Journal of Sound and Vibration, 225, 935-988. https://doi.org/10.1006/jsvi.1999.2257
[3]	Surana, K. and Alverio, E. (2020) Consistency and Validity of the Mathematical Models and the Solution Methods for BVPs and IVPs Based on Energy Methods and Principle of Virtual Work for Homogeneous Isotropic and Non-Homogeneous Non-Isotropic Solid Continua. Applied Mathematics, 11, 546-578. https://doi.org/10.4236/am.2020.117039
[4]	Cooman, A., Ferranti, F., Rolain, Y. and Vanderstein, G. (2016) Common-Denominator Modelling for Stability Analysis of Electronic Circuits. 11th European Microwave Integrated Circuits Conference, London, 3-4 October 2016, 361-364. https://doi.org/10.1109/EuMIC.2016.7777565
[5]	Mason, J. and Handscomb, D. (2002) Chebyshev Polynomials. Chapman and Hall/ CRC, London. https://doi.org/10.1201/9781420036114