Question:
Construct a Pade approximation to the exponential, using a quotient of two quadratics. Find the roots of the quadratic in the denominator. Use the two quadratics and the roots of the denominator quadratic to describe a numerical method for solving an initial-value problem with constant coefficients. What order do you expect this method to have? (You do not need to program the method.)
Discussion:
First we assume the pade approximation R(x) to the exponential function $e^x$ using a quotient of two quadratics is \[ e^x \approx R(x) = \frac{P(x)}{Q(x)}\] \[\text{ where } P(x) = a_0+a_1x+a_2x^3 \text{ and } Q(x) = 1+b_1x+b_2x^2. \] Taylor expansion of the exponential yields $e^x = 1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+O(x^5).$ Therefore \[(1+b_1x+b_2x^2) (1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+O(x^5)) = a_0+a_1x+a_2x^3.\] Comparing terms of same orders, we get \[ O(x^0): 1=a_0;\] \[ O(x^1): 1+b_1=a_1;\] \[ O(x^2): b_1+\frac{1}{2}+b2=a_2;\] \[ O(x^3): \frac{1}{6}+\frac{b_1}{2}+b_2=0;\] \[ O(x^4): \frac{1}{24}+\frac{b_1}{6}+\frac{b_2}{2}=0.\] From the above relations the coefficients of the pade approximation can be determined as \[ a_0=1,a_1=\frac{1}{2},a_2=\frac{1}{12},b_1=-\frac{1}{2},b_2=\frac{1}{12}\] and \[ R(x) = \frac{1+\frac{1}{2}x+\frac{1}{12}x^2}{1-\frac{1}{2}x+\frac{1}{12}x^2.}\] This exponential approximation leads to the recurrence \[(1-\frac{1}{2}Ah+\frac{1}{12}(Ah)^2)y_k =( 1+\frac{1}{2}Ah+\frac{1}{12}(Ah)^2)y_{k-1} \] for the solution of $y'(t) = Ay(t).$
Since the roots of the denominator quadratic $ 1-\frac{1}{2}x+\frac{1}{12}x^2 $ are $3\pm \sqrt{3}i$, we know that the matrix $(1-\frac{1}{2}Ah+\frac{1}{12}(Ah)^2)$ is always invertible.
Then given that $X_0 = I,$ the recurrence leads to \[ e^{tA} \approx X_m = [(1-\frac{1}{2}Ah+\frac{1}{12}(Ah)^2)^{-1}( 1+\frac{1}{2}Ah+\frac{1}{12}(Ah)^2)]^m \] where $t = mh.$
Referred to the discussion of the numerical methods for initial value problems with constant coefficients, we can use this exponential approximation to construct the following method. Suppose we have the linear ODE system with constant coefficients \[y'(t)=Ay(t)+b, y(t)=\xi.\] M(t) is approximated by the above recurrence result $X_m$. Then $ y(t) = M(t)[\xi+\int_0^{t}M(s)^{-1}b ds] $ solves the IVP. The integral part in the solution can be solved numerically via different approaches.
Since the pade approximation for the exponential that we used is fourth order, I expect the method to have fourth order accuracy as long as we choose the numerical integration method that is at least fourth order.
Last modified: