3.4 Cubic spline time integral
The time interval \(\left[0,T\right]\) is divided into\(N\) intervals with the length \(k=T/N\) and the time node is\(\ \tau_{n}=\text{nk},\text{nϵ}\left[0,N\right]\). At each time interval \(\left[\tau_{n},\tau_{n+1}\right]\)three spline approximations are found for the given values\(V\left(\tau_{n}\right)\)and the differential terms\(\ V^{{}^{\prime}}\left(\tau_{n}\right)V^{{}^{\prime}}\left(\tau_{n}+\frac{k}{2}\right)V^{{}^{\prime}}\left(\tau_{n}+k\right)\), that is
\begin{equation} V\left(\tau\right)=V\left(\tau_{n}\right)+k\rho_{0}\left(q\right)V^{{}^{\prime}}\left(\tau_{n}\right)+k\rho_{1}\left(q\right)V^{{}^{\prime}}\left(\tau_{n}+\frac{k}{2}\right)+k\rho_{2}\left(q\right)V^{{}^{\prime}}\left(\tau_{n}+k\right),\nonumber \\ \end{equation}
where polynomial \(\rho_{i}\) are
\begin{equation} \text{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ρ}_{0}\left(q\right)=q-\frac{3}{2}q^{2}+\frac{2}{3}q^{3},\rho_{1}\left(q\right)=2q^{2}-\frac{4}{3}q^{3},\rho_{2}\left(q\right)=-\frac{1}{2}q^{2}+\frac{2}{3}q^{3}.\nonumber \\ \end{equation}
Therefore, we obtain
\(V\left(\tau_{n}+k\right)=V\left(\tau_{n}\right)+\frac{k}{6}\left(V^{{}^{\prime}}\left(\tau_{n}\right)+4V^{{}^{\prime}}\left(\tau_{n}+\frac{k}{2}\right)+V^{{}^{\prime}}\left(\tau_{n}+k\right)\right),\)(3.4)
\(V\left(\tau_{n}+\frac{k}{2}\right)=V\left(\tau_{n}\right)+\frac{k}{24}\left(5V^{{}^{\prime}}\left(\tau_{n}\right)+8V^{{}^{\prime}}\left(\tau_{n}+\frac{k}{2}\right)-V^{{}^{\prime}}\left(\tau_{n}+k\right)\right).\)(3.5)
By use of (3.3) and (3.5), we have
\(\left(I-\frac{k}{3}C\right)V\left(\tau_{n}+\frac{k}{2}\right)=\left(I+\frac{5k}{24}C\right)V\left(\tau_{n}\right)-\frac{k}{24}\text{CV}\left(\tau_{n}+k\right)+_{1}\left(\tau_{n}\right)\),
where\(\ {}_{1}\left(\tau_{n}\right)=\frac{k}{24}\left(5b\left(\tau_{n}\right)+8b\left(\tau_{n}+\frac{k}{2}\right)-\theta\left(\tau_{n}+k\right)\right).\)
Removing\(\ V\left(\tau_{n}+\frac{k}{2}\right)\), (3.4) can be written as the following form
\begin{equation} V\left(\tau_{n}+k\right)=V\left(\tau_{n}\right)+_{2}\left(\tau_{n}\right)+\frac{k}{6}C\left(V\left(\tau_{n}\right)+4V\left(\tau_{n}+\frac{k}{2}\right)+V\left(\tau_{n}+k\right)\right),\nonumber \\ \end{equation}
where\({}_{2}\left(\tau_{n}\right)=\frac{k}{6}\left(\theta\left(\tau_{n}\right)+4\theta\left(\tau_{n}+\frac{k}{2}\right)+\theta\left(\tau_{n}+k\right)\right).\)
Time marching algorithm
\begin{equation} \text{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\left(I-\frac{k}{2}C+\frac{k^{2}}{12}C^{2}\right)V\left(\tau_{n}+k\right)=\left(I+\frac{k}{2}C+\frac{k^{2}}{12}C^{2}\right)\ V\left(\tau_{n}\right)+c\left(\tau_{n}\right),\nonumber \\ \end{equation}
where\(\ c\left(\tau_{n}\right)=\frac{2k}{3}C_{1}\left(\tau_{n}\right)+\left(I-\frac{k}{3}C\right)_{2}\left(\tau_{n}\right).\)
Although \({}_{1}\left(\tau_{n}\right)\) and\({}_{2}\left(\tau_{n}\right)\) depend on values at time\(\ \tau_{n}+\frac{k}{2}\) and\({\ \tau}_{n}+k\), these values can be easily calculated because they only depend on terminal conditions.