We split the matrix into five pieces, M=Me+Mo+Mt1+Mt2+Mt3, where
In the same manner as the tridiagonal matrix, we can approximate the time-stepping operator as
Using the series definition of the exponential function, we see that
Exponentiating Mt1, Mt2, and Mt3 amounts to exponentiating F. The terms , and are also block diagonal. Since the eigenvalues of are 1 and with imaginary c, it follows that .As before, to prove the unconditional stability of the algorithm, we need only show that . This follows immediately since
each of which is equal to 1 according to the preceding argument.
In Figure , we compared the impulse responses of first-order approximation with the second-order approximation in the Taylor-series expansion of the square-root operator. We also superposed semicircles which is the theoretical solution of the extrapolation operator on Fig and the higher order approximation shows better fitting to semicircles than the first-first order approximation.
With the same manner as showed in this section, we can get more accurate operator by taking more terms in a Taylor series expansion of the square-root operator. However, it will produce a matrix with increasing width of the band and thus will cause to an increasing in computation cost.