On the Approximated Solution of a Special Type of Nonlinear Third-Order Matrix Ordinary Differential Problem

Abstract

Matrix differential equations are at the heart of many science and engineering problems. In this paper, a procedure based on higher-order matrix splines is proposed to provide the approximated numerical solution of special nonlinear third-order matrix differential equations, having the form $Y^{\left(3\right)}\left(x\right)=f(x,Y\left(x\right))$. Some numerical test problems are also included, whose solutions are computed by our method.

1. Introduction

In different areas of engineering and applied science, such as those described in [1,2], third-order ordinary differential equation arise in the form of:

$$y^{\left(3\right)}=f(x,y),x\in [a,b],$$

(1)

with initial conditions $y\left(a\right)=y_a,y^{\prime }\left(a\right)=y_a^{\prime }$, and $y^{″}\left(a\right)=y_a^{″}$. Here, first and second derivatives are not part of the equation to be solved.

For the resolution of this type of equation, various different numerical techniques have been developed and exist in the literature [3,4,5,6,7]. Nonetheless, they frequently depend on very well-known numerical methods, once the original equations have been converted into a system of first-order ordinary differential equations. However, this conversion requires an increase in the number of unknowns of the problem, with the unavoidable growth in the computational cost that this entails. On the other hand, it is possible to use direct-integration methods to solve higher-order differential equations, such as (1), which have been proven to be computationally efficient and offer results with a very satisfactory accuracy (see [8,9,10] and the references therein).

In this paper, we present a new algorithm for solving matrix third-order differential equations of the type shown in Equation (1). Regarding the mathematical notation to be used hereinafter, we decided to adopt the formulation commonly used in matrix calculus, which is also part of the papers [11,12], previously written by the authors. Accordingly, given a vector $x\in {\mathbb{C}}^n$, its Euclidean norm is defined as ${\parallel x\parallel }_2=({x^{\top }x)}^{1/2}$ and its 1-norm as ${\displaystyle {\parallel x\parallel }_1=\sum _{i=1}^s\left|x_i\right|}$. Likewise, for a rectangular $m\times n$ matrix $A\in {\mathbb{C}}^{m\times n}$, its 2-norm is:

$${∥A∥}_2=\underset{x\ne 0}{sup}\frac{{∥Ax∥}_2}{{∥x∥}_2}.$$

Suppose $A\in {\mathbb{C}}^{m\times n}$ and $B\in {\mathbb{C}}^{p\times q}$. Then, the Kronecker product $A\otimes B$ is the $pm\times qn$ block matrix:

$$A\otimes B=\left(\begin{array}{ccc}{a}_{11}B& \dots & a_{1n}B\\ ⋮& & ⋮\\ a_{m1}B& \dots & a_{mn}B\end{array}\right).$$

If $A\in {\mathbb{C}}^{m\times n}$, the column-vector operator on it is defined by:

$$\mathbf{vec}\left(A\right)=\left(\begin{array}{c}{\mathbf{A}}_{•1}\\ ⋮\\ {\mathbf{A}}_{•n}\end{array}\right),\mathrm{being}{\mathbf{A}}_{•k}=\left(\begin{array}{c}{a}_{1k}\\ ⋮\\ a_{mk}\end{array}\right).$$

The derivative of a matrix $A\in {\mathbb{C}}^{m\times n}$ with respect to a matrix $B\in {\mathbb{C}}^{p\times q}$ was defined by [13] as:

$$\frac{\partial A}{\partial B}=\left(\begin{array}{ccc}{\displaystyle \frac{\partial A}{\partial b_{11}}}& \dots & {\displaystyle \frac{\partial A}{\partial b_{1q}}}\\ ⋮& & ⋮\\ {\displaystyle \frac{\partial A}{\partial b_{p1}}}& \dots & {\displaystyle \frac{\partial A}{\partial b_{pq}}}\end{array}\right),\mathrm{where}\frac{\partial A}{\partial b_{ij}}=\left(\begin{array}{ccc}{\displaystyle \frac{\partial a_{11}}{\partial b_{ij}}}& \dots & {\displaystyle \frac{\partial a_{1n}}{\partial b_{ij}}}\\ ⋮& & ⋮\\ {\displaystyle \frac{\partial a_{m1}}{\partial b_{ij}}}& \dots & {\displaystyle \frac{\partial a_{mn}}{\partial b_{ij}}}\end{array}\right).$$

In addition, the derivative of the product of two matrices $A\in {\mathbb{C}}^{m\times n}$ and $B\in {\mathbb{C}}^{n\times s}$ with respect to another matrix $C\in {\mathbb{C}}^{p\times q}$ is given by:

$$\frac{\partial AB}{\partial C}=\frac{\partial A}{\partial C}\left[I_q\otimes B\right]+\left[I_p\otimes A\right]\frac{\partial B}{\partial C},$$

(2)

$I_q$ and $I_p$ being, respectively, identity matrices of orders q and p. Moreover, given the above matrices, even the chain rule is also applicable:

$$\frac{\partial C}{\partial A}=\left[\frac{\partial {\left[\mathbf{vec}\left(B\right)\right]}^{\top }}{\partial A}\otimes I_p\right]\left[I_n\otimes \frac{\partial C}{\partial \left[\mathbf{vec}\left(B\right)\right]}\right].$$

(3)

This paper is organized as follows. In Section 2, the proposed method is described in detail. Then, in Section 3, an algorithm for solving this type of equation and its corresponding implementation in MATLAB are included. Next, in Section 4, distinct problems are numerically solved via the suggested method. Finally, conclusions are given in Section 5.

2. Description of the Method

Let the following third-order matrix differential equation be:

$$\left.\begin{array}{ccc}\hfill Y^{\left(3\right)}\left(x\right)& =& f(x,Y(x\left)\right)\hfill \\ \hfill Y\left(a\right)& =& Y_a\hfill \\ \hfill Y^{\prime }\left(a\right)& =& Y_a^{\prime }\hfill \\ \hfill Y^{″}\left(a\right)& =& Y_a^{″}\hfill \end{array}\right\},a\le x\le b,$$

(4)

where $Y\left(x\right)\in {\mathbb{C}}^{p\times q}$ is the unknown matrix and $Y_a,Y_a^{\prime },Y_a^{″}\in {\mathbb{C}}^{p\times q}$ are matrices that collect the initial solution. $f:[a,b]\times {\mathbb{C}}^{p\times q}\to {\mathbb{C}}^{p\times q}$ is a matrix-valued function and, as the differentiability class, $f\in {\mathcal{C}}^s\left(T\right)$, being $s\ge 1$, with:

$$T=\left\{(x,Y);a\le x\le b,Y\in {\mathbb{C}}^{p\times q}\right\}.$$

(5)

To ensure the uniqueness and the existence of $Y\left(x\right)$ as the differentiable continuously solution of problem (4) (see [14], p. 99), f also complies with the global Lipschitz condition:

$$∥f\left(x,Y_1\right)-f\left(x,Y_2\right)∥\le L∥Y_1-Y_2∥,a\le x\le b,Y_1,Y_2\in {\mathbb{C}}^{p\times q},$$

(6)

where $L>0$ is the Lipschitz constant.

The initial interval $[a,b]$, in which the solution must be provided, is partitioned into n subintervals, each of size $h=(b-a)/n$, as follows:

$${\Delta }_{[a,b]}=\left\{a=x_0<x_1<\dots <x_n=b\right\},x_i=a+ih,i=0,1,\dots ,n.$$

(7)

In each of the n aforementioned intervals and taking s as the order of the differentiability class of the function f, a matrix spline $S\left(x\right)$ of order $m\in \mathbb{N}$, with $3\le m\le s$, is defined. Thus, the differential Equation (4) is approximately solved so that $S\left(x\right)\in {\mathcal{C}}^3\left([a,b]\right)$.

For the first subinterval, the spline is designed as:

$$\begin{array}{ccc}\hfill S_{\left|{}_{\left[a,a+h\right]}\right.}\left(x\right)& =& Y\left(a\right)+Y^{\prime }\left(a\right)(x-a)+\frac{1}{2!}{Y}^{″}\left(a\right){(x-a)}^2+\frac{1}{3!}{Y}^{\left(3\right)}\left(a\right){(x-a)}^3\hfill \\ & +& \frac{1}{4!}{Y}^{\left(4\right)}\left(a\right){(x-a)}^4+\dots +\frac{1}{(m-1)!}{Y}^{(m-1)}\left(a\right){(x-a)}^{m-1}\hfill \\ \hfill & +& \frac{1}{m!}{A}_0{(x-a)}^m.\hfill \end{array}$$

(8)

As can be easily verified, this spline fulfills Equation (4) at $x=a$, since:

$$S_{\left|{}_{\left[a,a+h\right]}\right.}\left(a\right)=Y\left(a\right)=Y_a,S_{\left|{}_{\left[a,a+h\right]}\right.}^{\prime }\left(a\right)=Y^{\prime }\left(a\right)=Y_a^{\prime },S_{\left|{}_{\left[a,a+h\right]}\right.}^{″}\left(a\right)=Y^{″}\left(a\right)=Y_a^{″}$$

and:

$$S_{\left|{}_{\left[a,a+h\right]}\right.}^{\left(3\right)}\left(a\right)=Y^{\left(3\right)}\left(a\right)=f(a,Y\left(a\right))=f(a,S_{\left|{}_{\left[a,a+h\right]}\right.}\left(a\right)).$$

Obviously, the values corresponding to $Y^{\left(4\right)}\left(a\right),Y^{\left(5\right)}\left(a\right),Y^{\left(6\right)}\left(a\right),\dots ,Y^{(m-1)}\left(a\right)$ and the matrix $A_0\in {\mathbb{C}}^{p\times q}$ must be previously determined to fully define the spline (8). Based on this, the procedure given in [12] can be followed to calculate the fourth-order derivative $Y^{\left(4\right)}\left(x\right)$, that is,

$$\begin{array}{ccc}{Y}^{\left(4\right)}\left(x\right)\hfill & =\hfill & \frac{\partial f(x,Y(x\left)\right)}{\partial x}+\left[{\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }\otimes I_p\right]\frac{\partial f(x,Y(x\left)\right)}{\partial \mathbf{vec}Y\left(x\right)}\hfill \\ \hfill & =\hfill & g_1\left(x,Y\left(x\right),Y^{\prime }\left(x\right)\right),\hfill \end{array}$$

(9)

where $g_1\in {\mathcal{C}}^{s-1}\left(T\right)$. As a result, $Y^{\left(4\right)}\left(a\right)$ can be estimated as:

$$Y^{\left(4\right)}\left(a\right)=g_1\left(a,Y\left(a\right),Y^{\prime }\left(a\right)\right)=g_1\left(a,Y_a,Y_a^{\prime }\right).$$

Let us assume that $f\in {\mathcal{C}}^s\left(T\right)$, for $s\ge 2$. Thus, the second partial derivatives of f exist and are also continuous. The fifth-order derivative $Y^{\left(5\right)}\left(x\right)$ can be expressed as:

$$\begin{array}{c}{Y}^{\left(5\right)}\left(x\right)=\frac{{\partial }^{2}f(x,Y(x\left)\right)}{\partial x^2}+\left({\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }\otimes I_p\right)\frac{\partial }{\partial x}\left(\frac{\partial f(x,Y(x\left)\right)}{\partial \mathbf{vec}Y\left(x\right)}\right)\hfill \\ +\left(\frac{\partial {\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }}{\partial x}\otimes I_p\right)\frac{\partial f(x,Y(x\left)\right)}{\partial \mathbf{vec}Y\left(x\right)}\hfill \\ +\left({\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }\otimes I_p\right)\frac{\partial }{\partial \mathbf{vec}Y\left(x\right)}\left(\frac{\partial f(x,Y(x\left)\right)}{\partial x}\right)\hfill \\ +\left({\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }\otimes I_p\right)\left(\frac{\partial {\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }}{\partial \mathbf{vec}Y\left(x\right)}\otimes I_p\right)\frac{\partial f(x,Y(x\left)\right)}{\partial \mathbf{vec}Y\left(x\right)}\hfill \\ +\left({\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }\otimes I_p\right)\left({\left[\mathbf{vec}f(x,Y(x\left)\right)\right]}^{\top }\otimes I_{p^{2}q}\right)\frac{{\partial }^{2}f(x,Y(x\left)\right)}{{(\partial \mathbf{vec}Y(x\left)\right)}^2}\hfill \\ =g_2\left(x,Y\left(x\right),Y^{\prime }\left(x\right),Y^{″}\left(x\right)\right)\in {\mathcal{C}}^{s-2}\left(T\right).\hfill \end{array}$$

(10)

Therefore, $Y^{\left(5\right)}\left(a\right)=g_2\left(a,Y\left(a\right),Y^{\prime }\left(a\right),Y^{″}\left(a\right)\right)=g_2\left(a,Y_a,Y_a^{\prime },Y_a^{″}\right)$ can be immediately evaluated. For the calculation of all the other higher-order derivatives, from $Y^{\left(6\right)}\left(x\right)$ to $Y^{(m-1)}\left(x\right)$, the same previously described procedure can be carried out. In this way, we have:   

$$\left.\begin{array}{ccc}\hfill Y^{\left(6\right)}\left(x\right)& =& g_3\left(x,Y\left(x\right),Y^{\prime }\left(x\right),Y^{″}\left(x\right),Y^{‴}\left(x\right)\right)\in {\mathcal{C}}^{s-3}\left(T\right)\hfill \\ \hfill & ⋮& \\ \hfill Y^{(m-1)}\left(x\right)& =& g_{m-4}\left(x,Y\left(x\right),Y^{\prime }\left(x\right),\dots ,Y^{(m-4)}\left(x\right)\right)\in {\mathcal{C}}^{s-(m-4)}\left(T\right)\hfill \end{array}\right\}.$$

(11)

At this point, $Y^{\left(6\right)}\left(a\right),\dots ,Y^{(m-1)}\left(a\right)$ can be determined by just replacing $x=a$ in (11).

Once all these derivatives have been worked out, the matrix parameter $A_0$ should be determined as well. For this purpose, suppose that (8) is the solution of Equation (4) at the point $x=a+h$. Then, one obtains:

$$S_{\left|{}_{\left[a,a+h\right]}\right.}^{\left(3\right)}(a+h)=f\left(a+h,S_{\left|{}_{\left[a,a+h\right]}\right.}(a+h)\right),$$

(12)

Precisely, it is from here that the following implicit matrix equation can be formulated, with $A_0$ being its only unknown:

$$\begin{array}{cc}\hfill & A_0=\frac{(m-3)!}{h^{m-3}}\left[f\left(a+h,Y\left(a\right)+Y^{\prime }\left(a\right)h+\dots +\frac{h^{m-1}}{(m-1)!}{Y}^{(m-1)}\left(a\right)+\frac{h^m}{m!}{A}_0\right)\right.\hfill \\ \hfill & -\left.Y^{\left(3\right)}\left(a\right)-Y^{\left(4\right)}\left(a\right)h-\dots -\frac{1}{(m-4)!}{Y}^{(m-1)}\left(a\right)h^{m-4}\right].\hfill \end{array}$$

(13)

Once this equation has been solved for $A_0$, under the assumption that this solution is unique, all the parameters that form the spline in the interval $[a,a+h]$ will have been determined, and its analytical expression can be obtained.

For the next interval $[a+h,a+2h]$, the matrix spline can be expressed as:

$${\displaystyle \begin{array}{cc}{S}_{\left|{}_{\left[a+h,a+2h\right]}\right.}\left(x\right)& ={\sum }_{j=0}^2\frac{S_{\left|{}_{\left[a,a+h\right]}\right.}^{\left(j\right)}(a+h)}{i!}{(x-(a+h))}^j+\sum _{j=3}^{m-1}\frac{\overline{Y^{\left(j\right)}(a+h)}}{j!}{(x-(a+h))}^j\\ & +\frac{A_1}{m!}{(x-(a+h))}^m,\end{array}}$$

(14)

where

$$\overline{Y^{\left(3\right)}(a+h)}=f\left(a+h,S_{\left|{}_{\left[a,a+h\right]}\right.}(a+h)\right).$$

(15)

After having evaluated the corresponding derivates of $Y\left(x\right)$ using $S_{\left|{}_{\left[a,a+h\right]}\right.}(a+h)$ in (9)–(11), the expressions for $\overline{Y^{\left(4\right)}(a+h)},\dots ,\overline{Y^{(m-1)}(a+h)}$ are obtained, and they can be rewritten in the form:

$$\left.\begin{array}{ccc}\overline{Y^{\left(4\right)}(a+h)}\hfill & =& g_1\left(a+h,S_{\left|{}_{\left[a,a+h\right]}\right.}(a+h),S_{\left|{}_{\left[a,a+h\right]}\right.}^{\prime }(a+h)\right)\hfill \\ \hfill & ⋮& \\ \overline{Y^{(m-1)}(a+h)}\hfill & =& g_{m-4}\left(a+h,S_{\left|{}_{\left[a,a+h\right]}\right.}(a+h),\dots ,S_{\left|{}_{\left[a,a+h\right]}\right.}^{(m-4)}(a+h)\right)\hfill \end{array}\right\}.$$

(16)

The matrix spline $S\left(x\right)$, defined by means of (8) and (14), is of differentiability class ${\mathcal{C}}^3\left([a,a+h]\cup [a+h,a+2h]\right)$, and therefore, it satisfies Equation (4) at the point $x=a+h$. Similar to the previous interval, all spline coefficients have already been identified, with the exception of $A_1\in {\mathbb{R}}^{p\times q}$. Nevertheless, its value can be determined if the mentioned spline is considered as a solution of (4) at $x=a+2h$:

$$S_{\left|{}_{\left[a+h,a+2h\right]}\right.}^{\left(3\right)}(a+2h)=f\left(a+2h,S_{\left|{}_{\left[a+h,a+2h\right]}\right.}(a+2h)\right),$$

which leads to the following matrix equation, whose only unknown is precisely $A_1$:

$$\begin{array}{cc}{A}_1& =\frac{(m-3)!}{h^{m-3}}\left[f\left(a+2h,{\sum }_{j=0}^2\frac{S_{\left|{}_{\left[a,a+h\right]}\right.}^{\left(j\right)}(a+h)}{j!}{h}^j+{\sum }_{j=3}^{m-1}\frac{\overline{Y^{\left(j\right)}(a+h)}}{j!}{h}^j\right.\right.\\ & +\left.\frac{A_1{h}^m}{m!}\right)\left.-\overline{Y^{\left(3\right)}(a+h)}-\overline{Y^{\left(4\right)}(a+h)}h-\dots -\frac{h^{m-4}}{(m-4)!}\overline{Y^{(m-1)}(a+h)}\right].\end{array}$$

(17)

Assuming once again the uniqueness of the solution of the above equation for $A_1$, the spline is completely specified in the interval $[a+h,a+2h]$.

Now, repeating this procedure iteratively, the interval $\left[a+ih,a+(i+1)h\right]$ is reached, whose resulting spline then would be defined in the form:

$$\begin{array}{ccc}\hfill S_{\left|{}_{\left[a+ih,a+(i+1)h\right]}\right.}\left(x\right)& =& {\displaystyle \sum _{j=0}^2\frac{S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}^{\left(j\right)}(a+ih)}{j!}{(x-(a+ih))}^j}\hfill \\ & +& \sum _{j=3}^{m-1}\frac{\overline{Y^{\left(j\right)}(a+ih)}}{j!}{(x-(a+ih))}^j\hfill \\ & +& \frac{A_i}{m!}{(x-(a+ih))}^m,\hfill \end{array}$$

(18)

where:

$$\overline{Y^{\left(3\right)}(a+ih)}=f\left(a+ih,S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}(a+ih)\right),$$

(19)

and, similarly to what has been described previously,

$$\left.\begin{array}{ccc}\overline{Y^{\left(4\right)}(a+ih)}\hfill & =& g_1\left(a+ih,S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}(a+ih),S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}^{\prime }(a+ih)\right)\hfill \\ \hfill & ⋮& \\ \overline{Y^{(m-1)}(a+ih)}\hfill & =& g_{m-4}\left(a+ih,S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}(a+ih),\dots ,S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}^{(m-4)}(a+ih)\right)\hfill \end{array}\right\}.$$

(20)

Under this premise, the matrix spline ${\displaystyle S\left(x\right)\in {\mathcal{C}}^3\left(\bigcup _{j=0}^i[a+jh,a+(j+1)h]\right)}$ satisfies the differential Equation (4) at the point $x=a+ih$. If we additionally consider that $S_{\left|{}_{\left[a+ih,a+(i+1)h\right]}\right.}\left(x\right)$ fulfills the mentioned equation at point $x=a+(i+1)h$, then:

$$S_{\left|{}_{\left[a+ih,a+(i+1)h\right]}\right.}^{\left(3\right)}(a+(i+1)h)=f\left(a+(i+1)h,S_{\left|{}_{\left[a+ih,a+(i+1)h\right]}\right.}(a+(i+1)h)\right).$$

After expanding this previous formulation, one arrives at:

$$\begin{array}{ccc}{A}_i\hfill & =\hfill & \frac{(m-3)!}{h^{m-3}}\left[f\left({\displaystyle a+(i+1)h,{\sum }_{j=0}^2\frac{S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}^{\left(j\right)}(a+ih)}{j!}{h}^j+\sum _{j=3}^{m-1}\frac{\overline{Y^{\left(j\right)}(a+ih)}}{j!}{h}^j}\right.\right.\hfill \\ \hfill & +\hfill & \left.\left.\frac{A_i}{m!}{h}^m\right)-\overline{Y^{\left(3\right)}(a+ih)}-\dots -\frac{h^{m-4}}{(m-4)!}\overline{Y^{(m-1)}(a+ih)}\right].\hfill \end{array}$$

(21)

Note that Equations (13) and (17) can be derived from Equation (21) when $i=0$ or $i=1$, respectively.

Let us now demonstrate the uniqueness of the solution for these equations by using a fixed-point argument. Thus, let the matrix-valued function $g:{\mathbb{C}}^{p\times q}\to {\mathbb{C}}^{p\times q}$ be defined as follows for given values of h and i:

$$\begin{array}{ccc}g\left(T\right)\hfill & =\hfill & \frac{(m-3)!}{h^{m-3}}\left[f\left({\displaystyle a+(i+1)h,{\sum }_{j=0}^2\frac{S_{\left|{}_{\left[a+(i-1)h,a+ih\right]}\right.}^{\left(j\right)}(a+ih)}{j!}{h}^j+\sum _{j=3}^{m-1}\frac{\overline{Y^{\left(j\right)}(a+ih)}}{j!}{h}^j}\right.\right.\hfill \\ \hfill & +\hfill & \left.\left.\frac{T}{m!}{h}^m\right)-\overline{Y^{\left(3\right)}(a+ih)}-\dots -\frac{h^{m-4}}{(m-4)!}\overline{Y^{(m-1)}(a+ih)}\right].\hfill \end{array}$$

(22)

Equation (21) is satisfied if and only if $A_i=g\left(A_i\right)$, i.e., if the matrix $A_i$ is a fixed point for function $g\left(T\right)$. Considering this definition for function $g\left(T\right)$ and the global Lipschitz condition (6) for f, we have:

$$∥g\left(T_1\right)-g\left(T_2\right)∥\le \frac{L{h}^3}{m(m-1)(m-2)}∥T_1-T_2∥.$$

If we take $h<\sqrt[3]{m(m-1)(m-2)/L}$, the matrix function $g\left(T\right)$ will be contractive, and the solutions $A_i$, $i=0,1,\dots ,n-1$, of Equation (21) will be unique and, as a consequence, entirely determine the matrix spline. We can therefore state that the following theorem has been proven.

Theorem 1.

Take the third-order matrix differential equation expressed in (4) and the Lipschitz constant L described in (6). If the interval $[a,b]$ is partitioned as in (7) by considering a step size h such as:

$$h<\sqrt[3]{m(m-1)(m-2)/L},$$

then the matrix spline $S\left(x\right)$ of order $m,3\le m\le s$, s being the order of the differentiability class of the function f, is of class ${\mathcal{C}}^3[a,b]$, and it can be defined as in (18) for each subinterval $\left[a+ih,a+(i+1)h\right]$, $i=0,1,\dots ,n-1$.

If an analysis similar to the one performed in [15] were carried out, we could conclude that the defined splines have a global error of $O\left(h^{m-1}\right)$.

3. Algorithm for the MATLAB Program

Consider the following third-order ordinary differential equation for the k-th step:

$$\left.\begin{array}{ccc}\hfill Y^{\left(3\right)}\left(x\right)& =& f(x,Y(x\left)\right)\hfill \\ \\ \hfill Y\left(x_k\right)& =& Y_0,Y^{\prime }\left(x_k\right)=Y_1,Y^{″}\left(x_k\right)=Y_2\hfill \end{array}\right\}{x}_k\le x\le x_k+h,$$

(23)

where h is the step size, $x_k=a+kh$ and $Y_0$, $Y_1$, and $Y_2$ are, respectively, the matrices Y, $Y^{\prime }$, and $Y^{″}$ obtained in the previous step at $x_k$. If we denote by $S_k\left(x\right)$ the spline of order m in the interval $[x_k,x_k+h]$, then:

$$\begin{array}{c}{S}_k(x_k+h)=\sum _{i=0}^{m-1}\frac{Y^{\left(i\right)}\left(x_k\right)h^i}{i!}+\frac{h^m}{m!}{A}_k\equiv B_k^{\left(0\right)}+\frac{h^m}{m!}{A}_k,\hfill \end{array}$$

(24)

$$\begin{array}{c}{S}_k^{{}^{\prime }}(x_k+h)=\sum _{k=1}^{m-1}\frac{Y^{\left(i\right)}\left(x_k\right)h^{i-1}}{(i-1)!}+\frac{h^{m-1}}{(m-1)!}{A}_k=B_k^{\left(1\right)}+\frac{h^{m-1}}{(m-1)!}{A}_k,\hfill \end{array}$$

(25)

$$\begin{array}{c}{S}_k^{{}^{″}}(x_k+h)=\sum _{k=2}^{m-1}\frac{Y^{\left(i\right)}\left(x_k\right)h^{i-2}}{(i-2)!}+\frac{h^{m-2}}{(m-2)!}{A}_k\equiv B_k^{\left(2\right)}+\frac{h^{m-2}}{(m-2)!}{A}_k,\hfill \end{array}$$

(26)

$$\begin{array}{c}{S}_k^{{}^{‴}}(x_k+h)=\sum _{k=3}^{m-1}\frac{Y^{\left(i\right)}\left(x_k\right)h^{i-3}}{(i-3)!}+\frac{h^{m-3}}{(m-3)!}{A}_k\equiv B_k^{\left(3\right)}+\frac{h^{m-3}}{(m-3)!}{A}_k.\end{array}$$

(27)

If we substitute expressions (24) and (27) into differential Equation (23), we obtain:

$$B_k^{\left(3\right)}+\frac{h^{m-3}}{(m-3)!}{A}_k=f\left(x,B_k^{\left(0\right)}+\frac{h^m}{m!}{A}_k\right).$$

Then, the matrix $A_k$ can be found by using the following fixed-point iteration:

$$A_k=\frac{(m-3)!}{h^{m-3}}\left[f\left(x,B_k^{\left(0\right)}+\frac{h^m}{m!}{A}_k\right)-B_k^{\left(3\right)}\right].$$

(28)

Hence, the approximated values for Y, $Y^{\prime }$, and $Y^{″}$ at $x_k+h$ are:

$$\begin{array}{c}Y(x_k+h)=B_k^{\left(0\right)}+\frac{h^m}{m!}{A}_k,\hfill \\ Y^{\prime }(x_k+h)=B_k^{\left(1\right)}+\frac{h^{m-1}}{(m-1)!}{A}_k,\hfill \\ Y^{″}(x_k+h)=B_k^{\left(2\right)}+\frac{h^{m-2}}{(m-2)!}{A}_k.\hfill \end{array}$$

Figure 1 shows the MATLAB code that approximately computes the solution matrices $Y\left(b\right)$, $Y^{\prime }\left(b\right)$, and $Y^{″}\left(b\right)$ of the third-order ordinary differential Equation (4). This code uses the cell-array data type for storing sets of matrices. In Line 31, $Y^{\left(i\right)}\left(x_k\right),i=0,1,\dots ,m$, are obtained and stored in the cell-array variable ym by invoking the f MATLAB function to return the matrices that appear in expressions (24)–(27). In Lines 11–16, the expressions $B_k^{\left(0\right)}$, $B_k^{\left(1\right)}$, $B_k^{\left(2\right)}$, and $B_k^{\left(3\right)}$ from Equations (24)–(27) are computed. In Lines 38–47, the matrix $A_k$ is worked out by using the fixed-point iteration. Finally, in Lines 48–50, $Y(x_k+h)$, $Y^{\prime }(x_k+h)$, and $Y^{″}(x_k+h)$ are computed approximately.

Figure 2 reproduces the MATLAB function f for the thin-film flow of problem (33). The dotted Line 24 indicates the completion of the list of all derivates until the m-th one, m being the spline order.

The memory requirements for this function are $(m+10)$ matrices, i.e., m matrices for the cell-array variable ym, 4 matrices for the cell-array variable Bk, 3 matrices for the variables Ak, Ak1, and Skm, and 3 matrices for the cell-array variable y.

4. Numerical Test Examples

This section discusses the approximation for solutions of seven distinct third-order differential equations. For all of them, special Mathematica and MATLAB software packages were employed to obtain the solutions of the algebraic equations that emerge after applying the algorithm described in the preceding section.

4.1. A Scalar Inhomogeneous Non-Linear Problem

Consider this simple test problem from [1]:

$$\left.\begin{array}{ccc}\hfill y^{\left(3\right)}\left(x\right)& =& (cos\left(x\right)-1)cos\left(x\right)+y^2-1\hfill \\ \\ \hfill y\left(0\right)& =& 0\hfill \\ \\ \hfill y^{\prime }\left(0\right)& =& 1\hfill \\ \\ \hfill y^{″}\left(0\right)& =& 0\hfill \end{array}\right\},0\le x\le 1,$$

(29)

whose exact solution $y\left(x\right)=sin\left(x\right)$ is known. The function $f(x,y)=(cos\left(x\right)-1)cos\left(x\right)+y^2-1$ is not Lipschitz for all of ${\mathbb{R}}^2$. However, the particular region $D=(0,1)\times (-1,1)\subset {\mathbb{R}}^2$ can be chosen, because with constant $L=2$, function $f(x,y)$ is locally Lipschitz. Therefore, it is possible to apply the proposed method, obtaining sixth-order splines ($m=6$) by taking $h<\sqrt[3]{60}\approx 3.91$. The analytical expressions of the distinct polynomials that compose the spline for $h=0.1$ are listed in

Table 1. Approximation spline for the problem (29).

Subinterval	Approximation Polynomials
$[0,0.1]$	$-0.0000347107x^6+0.00833333x^5-0.166667x^3+x$
$[0.1,0.2]$	$-0.00031184x^6+0.00849961x^5-0.0000415695x^4-0.166661x^3-4.15695\times {10}^{-7}{x}^2+1.x-2.7713\times {10}^{-10}$
$[0.2,0.3]$	$0.0000383134x^6+0.00807943x^5+0.000168523x^4-0.166717x^3+7.9880\times {10}^{-6}{x}^2+0.999999x+2.2133\times {10}^{-8}$
$[0.3,0.4]$	$-0.00190054x^6+0.0115694x^5-0.00244893x^4-0.16567x^3-0.000227583x^2+1.00003x-1.3913\times {10}^{-6}$
$[0.4,0.5]$	$0.00470998x^6-0.00429589x^5+0.0134163x^4-0.174132x^3+0.00231086x^2+0.999621x+0.0000256854$
$[0.5,0.6]$	$-0.0205597x^6+0.0715132x^5-0.081345x^4-0.110957x^3-0.0213795x^2+1.00436x-0.000369154$
$[0.6,0.7]$	$0.0731793x^6-0.265947x^5+0.424846x^4-0.51591x^3+0.160849x^2+0.960625x+0.00400433$
$[0.7,0.8]$	$-0.277197x^6+1.20563x^5-2.15042x^4+1.88767x^3-1.10103x^2+1.31395x-0.0372171$
$[0.8,0.9]$	$1.02995x^6-5.06866x^5+10.3982x^4-11.4975x^3+6.93007x^2-1.256x+0.305443$
$[0.9,1.0]$	$-3.8489x^6+21.2771x^5-48.8798x^4+59.6361x^3-41.0851x^2+16.0295x-2.28737$

. The maximum errors incurred are shown in

Table 2. Approximation errors for the problem (29).

Subinterval	$[0,0.1]$	$[0.1,0.2]$	$[0.2,0.3]$	$[0.3,0.4]$	$[0.4,0.5]$
Max. error	$1.4872\times {10}^{-11}$	$6.3015\times {10}^{-11}$	$6.9036\times {10}^{-10}$	$6.2645\times {10}^{-10}$	$6.6951\times {10}^{-9}$
Subinterval	$[0.5,0.6]$	$[0.6,0.7]$	$[0.7,0.8]$	$[0.8,0.9]$	$[0.9,1.0]$
Max. error	$1.2629\times {10}^{-8}$	$6.8791\times {10}^{-8}$	$2.2188\times {10}^{-7}$	$8.8033\times {10}^{-7}$	$8.7895\times {10}^{-7}$

for each subinterval.

4.2. A Scalar Nonhomogeneous Linear Problem

Here, we consider the following equation, studied in [8]:

$$\left.\begin{array}{ccc}\hfill y^{\left(3\right)}\left(x\right)& =& y+cos\left(x\right)\hfill \\ \\ \hfill y\left(0\right)& =& 0\hfill \\ \\ \hfill y^{\prime }\left(0\right)& =& 0\hfill \\ \\ \hfill y^{″}\left(0\right)& =& 1\hfill \end{array}\right\},0\le x\le 1,$$

(30)

whose exact solution is given by ${\displaystyle y\left(x\right)=\frac{1}{2}\left(e^x-cos\left(x\right)-sin\left(x\right)\right)}$. In this case, function $f(x,y)=y+cos\left(x\right)$ is Lipschitz over all ${\mathbb{R}}^2$, with constant $L=1$. If we choose $h<\sqrt[3]{120}\approx 4.93$, sixth-order splines ($m=6$) will be achieved. The approximation spline for step size $h=0.1$ is given in

Table 3. Approximation spline for the problem (30).

Subinterval	Approximation Polynomials
$[0,0.1]$	$0.00142361x^6+0.166667x^3+0.5x^2$
$[0.1,0.2]$	$2.77797\times {10}^{-10}-1.66678\times {10}^{-8}x+0.5x^2+0.166661x^3+0.0000416695x^4-0.000166678x^5+0.00170141x^6$
$[0.2,0.3]$	$-2.1931\times {10}^{-8}+6.49596\times {10}^{-7}x+0.499992x^2+0.166717x^3-0.000166538x^4+0.000249737x^5+0.0013544x^6$
$[0.3,0.4]$	$1.39603\times {10}^{-6}-0.0000277096x+0.500228x^2+0.165666x^3+0.00245931x^4-0.0032514x^5+0.00329947x^6$
$[0.4,0.5]$	$-0.0000256191+0.000377517x+0.497696x^2+0.174109x^3-0.0133699x^4+0.0125778x^5-0.00329602x^6$
$[0.5,0.6]$	$0.000369394-0.00436264x+0.521397x^2+0.110906x^3+0.0814334x^4-0.0632648x^5+0.0219848x^6$
$[0.6,0.7]$	$-0.0040011+0.0393423x+0.339292x^2+0.515582x^3-0.424411x^4+0.273965x^5-0.0716901x^6$
$[0.7,0.8]$	$0.0372102-0.313898x+1.60086x^2-1.88741x^3+2.15022x^4-1.19726x^5+0.2786x^6$
$[0.8,0.9]$	$-0.305309+2.255x-6.42693x^2+11.4922x^3-10.3932x^4+5.07446x^5-1.02801x^6$
$[0.9,1.0]$	$2.28658-15.0243x+41.571x^2-59.6158x^3+48.8635x^4-21.2619x^5+3.84909x^6$

, and the corresponding maximum errors for each subinterval are shown in

Table 4. Approximation errors for the problem (30).

Subinterval	$[0,0.1]$	$[0.1,0.2]$	$[0.2,0.3]$	$[0.3,0.4]$	$[0.4,0.5]$
Max. error	$1.4881\times {10}^{-11}$	$6.3053\times {10}^{-11}$	$6.9130\times {10}^{-10}$	$6.3575\times {10}^{-10}$	$6.7183\times {10}^{-9}$
Subinterval	$[0.5,0.6]$	$[0.6,0.7]$	$[0.7,0.8]$	$[0.8,0.9]$	$[0.9,1.0]$
Max. error	$1.256\times {10}^{-8}$	$6.8932\times {10}^{-8}$	$2.2147\times {10}^{-7}$	$8.8065\times {10}^{-7}$	$3.2087\times {10}^{-6}$

.

4.3. A Scalar Nonlinear Problem

Our next test example focuses on the problem posed in [8]:

$$\left.\begin{array}{ccc}\hfill y^{\left(3\right)}\left(x\right)& =& -e^{-y}+3e^{-2y}-2e^{-3y}\hfill \\ \\ \hfill y\left(0\right)& =& log\left(2\right)\hfill \\ \\ \hfill y^{\prime }\left(0\right)& =& \frac{1}{2}\hfill \\ \\ \hfill y^{″}\left(0\right)& =& \frac{1}{4}\hfill \end{array}\right\},0\le x\le 1.$$

(31)

The exact solution to this problem is ${\displaystyle y\left(x\right)=log(e^x+1)}$. Considering the region $D=(0,1)\times (0,10)\subset {\mathbb{R}}^2$, function $f(x,y)=-e^{-y}+3e^{-2y}-2e^{-3y}$ has a partial derivative with respect to its second variable y bounded by $L=13$ in D. Thus, function $f(x,y)$ is locally Lipschitz for constant $L=13$. For $h<\sqrt[3]{\frac{120}{13}}\approx 2.09771$, our method is applied, providing this time sixth-order splines ($m=6$).

Table 5. Approximation spline for the problem (31).

Subinterval	Approximation Polynomials
$[0,0.1]$	$log\left(2\right)+0.5x+0.125x^2-0.00520833x^4+0.000346486x^6$
$[0.1,0.2]$	$0.693147+0.5x+0.125x^2+3.51111\times {10}^{-7}x3-0.00521097x^4+0.0000105333x^5+0.00032893x^6$
$[0.2,0.3]$	$0.693147+0.5x+0.125x^2+2.9433\times {10}^{-6}{x}^3-0.00522069x^4+0.0000299747x^5+0.000312729x^6$
$[0.3,0.4]$	$0.693147+0.500001x+0.124989x^2+0.0000484757x^3-0.00533452x^4+0.000181749x^5+0.00022841x^6$
$[0.4,0.5]$	$0.693148+0.499994x+0.125035x^2-0.000102835x^3-0.00505081x^4-0.000101957x^5+0.000346621x^6$
$[0.5,0.6]$	$0.693137+0.500121x+0.124401x^2+0.00158826x^3-0.00758746x^4+0.00192736x^5-0.000329818x^6$
$[0.6,0.7]$	$0.693243+0.499064x+0.128805x^2-0.00819903x^3+0.00464666x^4-0.00622872x^5+0.00193576x^6$
$[0.7,0.8]$	$0.692217+0.50786x+0.0973911x^2+0.0516366x^3-0.059463x^4+0.0304054x^5-0.00678664x^6$
$[0.8,0.9]$	$0.700682+0.444372x+0.29579x^2-0.279029x^3+0.250536x^4-0.124594x^5+0.0255049x^6$
$[0.9,1.0]$	$0.636507+0.872202x-0.892627x^2+1.48159x^3-1.21665x^4+0.527487x^5-0.0952508x^6$

displays the approximation spline for $h=0.1$. In addition, the maximum errors for each subinterval are listed in

Table 6. Approximation errors for the problem (31).

Subinterval	$[0,0.1]$	$[0.1,0.2]$	$[0.2,0.3]$	$[0.3,0.4]$	$[0.4,0.5]$
Max. error	$4.7340\times {10}^{-13}$	$2.5465\times {10}^{-12}$	$3.9703\times {10}^{-11}$	$1.1270\times {10}^{-10}$	$3.9472\times {10}^{-10}$
Subinterval	$[0.5,0.6]$	$[0.6,0.7]$	$[0.7,0.8]$	$[0.8,0.9]$	$[0.9,1.0]$
Max. error	$4.6335\times {10}^{-10}$	$2.6369\times {10}^{-9}$	$3.8893\times {10}^{-9}$	$2.4363\times {10}^{-8}$	$7.5564\times {10}^{-8}$

.

4.4. A Scalar Nonhomogeneous Problem

An interesting scalar problem, being nonhomogeneous and involving convoluted combinations of trigonometric functions, is:    

$$\left.\begin{array}{ccc}\hfill y^{\left(3\right)}\left(x\right)& =& {\displaystyle -\frac{3}{2}ysin\left(2x\right)+sin\left(x\right)\left(1+{sin}^2\left(x\right)\right)cos(cos\left(x\right))}\hfill \\ \\ \hfill y\left(0\right)& =& sin\left(1\right)\hfill \\ \\ \hfill y^{\prime }\left(0\right)& =& 0\hfill \\ \\ \hfill y^{″}\left(0\right)& =& -cos\left(1\right)\hfill \end{array}\right\},0\le x\le 1.$$

(32)

Surprisingly, this complicated problem has an exact solution, for details see [8]. The known solution is ${\displaystyle y\left(x\right)=sin(cos\left(x\right))}$ (there are mistakes both in the initial conditions and in the solution of the problem, ${\displaystyle y\left(x\right)=cos(sin\left(x\right))}$, given in [8], p. 146). In this case, function:

$$f(x,y)=-\frac{3}{2}ysin\left(2x\right)+sin\left(x\right)\left(1+{sin}^2\left(x\right)\right)cos(cos\left(x\right))$$

is Lipschitz with constant $L=3/2$, and the proposed method will yield sixth-order splines ($m=6$) if $h<\sqrt[3]{\frac{240}{3}}\approx 4.30887$. For step size $h=0.1$, this spline is summarized in

Table 7. Approximation spline for problem (32).

Subinterval	Approximation Polynomials
$[0,0.1]$	$sin\left[1\right]-\frac{cos\left(1\right)}{2}{x}^2+\left(\frac{cos\left(1\right)}{24}-\frac{sin\left(1\right)}{8}\right)x^4+0.0279824x^6$
$[0.1,0.2]$	$0.841471+7.96286\times {10}^{-8}x-0.270153x^2+0.0000265429x^3-0.0828703x^4+0.000796286x^5+0.0266552x^6$
$[0.2,0.3]$	$0.841471+3.18879\times {10}^{-6}x-0.270192x^2+0.00028564x^3-0.083842x^4+0.00273951x^5+0.0250359x^6$
$[0.3,0.4]$	$0.841466+0.000109689x-0.27108x^2+0.00423011x^3-0.0937031x^4+0.0158877x^5+0.0177313x^6$
$[0.4,0.5]$	$0.841488-0.000229253x-0.268961x^2-0.00283119x^3-0.0804632x^4+0.0026478x^5+0.0232479x^6$
$[0.5,0.6]$	$0.840672+0.00956559x-0.317935x^2+0.127767x^3-0.27636x^4+0.159365x^5-0.0289912x^6$
$[0.6,0.7]$	$0.847811-0.0618259x-0.0204706x^2-0.533266x^3+0.549931x^4-0.391495x^5+0.124026x^6$
$[0.7,0.8]$	$0.774658+0.565197x-2.25984x^2+3.7322x^3-4.02021x^4+2.22001x^5-0.497762x^6$
$[0.8,0.9]$	$1.36838-3.88771x+11.6555x^2-19.4601x^3+17.7225x^4-8.65135x^5+1.76711x^6$
$[0.9,1.0]$	$-3.15253+26.2517x-72.0651x^2+104.57x^3-85.6362x^4+37.2859x^5-6.73979x^6$

. All corresponding maximum errors are also listed in

Table 8. Approximation errors for problem (32).

Subinterval	$[0,0.1]$	$[0.1,0.2]$	$[0.2,0.3]$	$[0.3,0.4]$	$[0.4,0.5]$
Max. error	$3.4875\times {10}^{-11}$	$1.811\times {10}^{-10}$	$2.989\times {10}^{-9}$	$9.2962\times {10}^{-9}$	$3.2810\times {10}^{-8}$
Subinterval	$[0.5,0.6]$	$[0.6,0.7]$	$[0.7,0.8]$	$[0.8,0.9]$	$[0.9,1.0]$
Max. error	$4.8765\times {10}^{-8}$	$2.1830\times {10}^{-7}$	$2.0686\times {10}^{-7}$	$1.845\times {10}^{-6}$	$5.095\times {10}^{-6}$

.

4.5. A Nonlinear Scalar Problem

The following problem consists of modeling the spreading of a thin oil drop on a horizontal surface. It is the thin-film problem given by Example 4.3 in [16]:

$$\left.\begin{array}{c}{y}^{\left(3\right)}=y^{-2}\\ y\left(0\right)=y^{\prime }\left(0\right)=y^{″}\left(0\right)=1\end{array}\right\},x\in [0,b].$$

(33)

Table 9. Relative errors of the Mechee ($M_1$), Khataybeh ($M_2$), and spline-based ($M_3$) methods in the film flow problem (33).

x	${\mathit{M}}_1$	${\mathit{M}}_2$	${\mathit{M}}_3$
0.2	$8.397402\times {10}^{-7}$	$8.399040\times {10}^{-7}$	$8.397170\times {10}^{-7}$
0.4	$7.596544\times {10}^{-8}$	$7.643561\times {10}^{-8}$	$7.598801\times {10}^{-8}$
0.6	$3.485744\times {10}^{-9}$	$4.315684\times {10}^{-9}$	$3.496411\times {10}^{-9}$
0.8	$4.587537\times {10}^{-11}$	$1.192760\times {10}^{-9}$	$3.625209\times {10}^{-11}$
1	$3.659122\times {10}^{-7}$	$3.674458\times {10}^{-7}$	$3.659145\times {10}^{-7}$

details the relative errors for the method by Mechee et al. ($M_1$) [17], the method by Khataybeh et al. ($M_2$), and our spline-based method of order $m=9$ ($M_3$), for $h=0.01$ and $b=1$. For the numerical comparison, the exact solution was obtained from Duffy and Wilson [1].

These results were produced with the MATLAB code shown in Figure 2, which calculates function y and its derivatives.

4.6. A Linear Matrix Problem

Consider the matrix problem:

$$\left.\begin{array}{ccc}\hfill Y^{\left(3\right)}\left(x\right)& =& AY\left(x\right),A=\left(\begin{array}{ccc}\hfill \frac{817}{68}& \hfill \frac{1393}{68}& \hfill \frac{448}{68}\\ \\ \hfill -\frac{1141}{68}& \hfill -\frac{2837}{68}& \hfill -\frac{896}{68}\\ \\ \hfill \frac{3059}{136}& \hfill \frac{4319}{136}& \hfill \frac{1592}{136}\end{array}\right)\hfill \\ \\ \hfill y\left(0\right)& =& \left(\begin{array}{c}\hfill 2\\ \hfill -2\\ \hfill 12\end{array}\right)\hfill \\ \\ \hfill y^{\prime }\left(0\right)& =& \left(\begin{array}{c}\hfill -12\\ \hfill 28\\ \hfill -33\end{array}\right)\hfill \\ \\ \hfill y^{″}\left(0\right)& =& \left(\begin{array}{c}\hfill 20\\ \hfill -52\\ \hfill 5\end{array}\right)\hfill \end{array}\right\},0\le x\le 1.$$

(34)

As in the other cases above, its solution (see [18], p. 334) is also known:

$${\displaystyle y\left(x\right)=\left(\begin{array}{c}{e}^x-2e^{2x}+3e^{-3x}\\ 3e^x+2e^{2x}-7e^{-3x}\\ -11e^x-5e^{2x}+4e^{-3x}\end{array}\right).}$$

The matrix function $f(x,Y)=AY$ is Lipschitz with constant $L={∥A∥}_2\approx 90.1136$. Therefore, for $h<\sqrt[3]{\frac{210}{90.1136}}\approx 1.32579$, we employ our method and seventh-order splines ($m=7$) are obtained. Maximum errors for $h=0.1$ are given in

Table 10. Approximation errors for the problem (34).

Subinterval	$[0,0.1]$	$[0.1,0.2]$	$[0.2,0.3]$	$[0.3,0.4]$	$[0.4,0.5]$
Max. error	$7.5906\times {10}^{-11}$	$7.2457\times {10}^{-11}$	$4.3608\times {10}^{-10}$	$1.4836\times {10}^{-9}$	$3.7673\times {10}^{-9}$
Subinterval	$[0.5,0.6]$	$[0.6,0.7]$	$[0.7,0.8]$	$[0.8,0.9]$	$[0.9,1.0]$
Max. error	$7.9945\times {10}^{-9}$	$1.5020\times {10}^{-8}$	$2.5835\times {10}^{-8}$	$4.1559\times {10}^{-8}$	$6.3425\times {10}^{-8}$

.

Figure 3 and Figure 4 display errors committed when considering distinct values of h or different spline orders, respectively. As expected, the errors become smaller as either the value of h decreases or as the spline order increases. MATLAB software was employed to provide these graphs.

4.7. A Nonlinear Matrix Problem

In the next example, we focus on the following matrix problem with a nonlinear term:

$$\left.\begin{array}{ccc}\hfill Y^{\left(3\right)}& =& Y^2\hfill \\ \hfill Y\left(0\right)& =& 0_n\hfill \\ \hfill Y^{\prime }\left(0\right)& =& {10}^{-2}{I}_n\hfill \\ \hfill Y^{″}\left(0\right)& =& {10}^{-2}{1}_n\hfill \end{array}\right\},x\in [0,b],$$

(35)

where $0_n$, $I_n$, and $1_n$ are the null, the identity, and the all-ones matrices of dimension n, respectively.

To compute the “exact” solution, we employed the MATLAB Symbolic Math Toolbox (with 256 digits of precision) using the vpa (variable-precision floating-point arithmetic) function. All computations have IEEE double-precision arithmetic with round off $u=2^{-53}\approx 1.11\times {10}^{-16}$. The “exact” solution is obtained whenever two consecutive spline orders (for fixed step size) present a relative error lower than the unit round off for this precision. For $b=10$ and $n=10$, we found that for $h=0.01$ and for the spline orders 10 and 11, the relative error when the “exact” solution is calculated is equal to $1.2019\times {10}^{-17}<u$.

Table 11. Relative errors for the nonlinear matrix problem (35) solved by a spline of order $m=12$.

x	2	4	6	8	10
$h=0.5$	$1.46\times {10}^{-14}$	$4.32\times {10}^{-13}$	$8.68\times {10}^{-12}$	$4.77\times {10}^{-10}$	$4.24\times {10}^{-07}$
$h=0.1$	$6.10\times {10}^{-16}$	$8.35\times {10}^{-16}$	$1.47\times {10}^{-15}$	$3.37\times {10}^{-15}$	$2.48\times {10}^{-13}$
$h=0.05$	$1.05\times {10}^{-15}$	$6.27\times {10}^{-16}$	$1.23\times {10}^{-15}$	$2.27\times {10}^{-15}$	$4.58\times {10}^{-15}$
$h=0.01$	$1.58\times {10}^{-15}$	$2.48\times {10}^{-15}$	$3.16\times {10}^{-15}$	$8.77\times {10}^{-15}$	$2.11\times {10}^{-14}$

shows the relative errors at $x=2$, 4, 6, 8, and 10, with a spline order $m=12$, considering the step sizes $h=0.1$, $h=0.05$, and $h=0.01$. As we can see, in general, the error committed at any point x becomes smaller and smaller for decreasing h values. Nevertheless, for $h=0.01$, the results obtained are worse than those ones for $h=0.05$, due to a larger number of floating-point operations and reaching the limit of the machine accuracy.

Figure 5 plots the relative errors at $x=10$ for distinct step sizes when varying the spline order. The relative errors decrease as we use smaller values of h and as the order of the spline m increases. When $h=0.05$, the error can even be larger if we use a very high-degree spline versus lower degrees. Moreover, for $h=0.01$, using a value of m greater than or equal to nine does not reduce the relative error committed. As explained for Table 11, this is because of the round off errors that appear when computing a large number of floating-point operations.

Additionally, Figure 6 displays the execution time required for different degrees of the spline and for distinct values of h. As expected, a higher value of h and m leads to an increase in the computational cost and, consequently, in the execution time. For the same value of h, it is possible that the execution time will not always rise with the degree of the spline. This is because the number of iterations required by the fixed-point method is not always identical.

Finally, Figure 7 depicts the response time versus the relative error committed, also with different step sizes. Clearly, higher values of h result in larger execution times. On the other hand, for a constant step size, a smaller relative error leads to a higher execution time, since this is achieved by working with higher spline orders.

5. Conclusions

This work outlined a numerical method for solving a specific type of third-order matrix differential equations and presents the corresponding algorithmic recipes. For illustration purposes, several standard test problems were included. All results demonstrated that our procedure provides solutions with very acceptable accuracy. Of course, the degree of accuracy can always be readjusted by decreasing the step size and/or increasing the order of the approximation spline accordingly.