Are there numerical methods to solve a differential equation which are actually faster than numerically computing its analytical solution?

Question

In the topic of numerical solutions of ODE and PDE, usually it's said that many times it's impractical to try to find an analytical solution due to the complexity of the boundary conditions, or even outright impossible due to the nature of the equation, for example: \begin{equation} y'=e^{-x^{2}} \end{equation} doesn't have an analytical solution we can write using other known functions. It is implied, though, that an analytical solution it's always preferable whenever possible to speed up calculations.

Given that known solutions of differential equations can turn out to be a really complicated mess of nested trigonometric and exponential functions (as an example, this is the solution to the heat equation in 2 spatial coordinates:

$$\theta(x,y)=\frac{2}{\pi}\sum_{n=1}^\infty\frac{(-1)^{n+1}+1}{n}\sin\frac{n\pi x}{L}\frac{\sinh(n\pi y/L)}{\sinh(n\pi W/L)}$$ ) and given that, for practical purposes, to output a value, the known functions like sines, cosine and exponentials which make up the solution eventually have to be computed numerically, can there be an equation for which the analytical solution is known, but so complex that solving the DE with a common numerical method like Euler or Runge-Kutta turns out to be actually faster than computing the analytical solution using Taylor series?

[Edit: It was pointed out that $ y'=e^{-x^{2}} $ does indeed have an analytical solution, just not a closed form in terms of elementary functions]

Answer 1

I've found a complete answer to my question in this paper:

Ardourel, Vincent; Jebeile, Julie, On the presumed superiority of analytical solutions over numerical methods, Eur. J. Philos. Sci. 7, No. 2, 201-220 (2017). ZBL1384.00034.

There are many differential equations for which an analytical solution is known, but the Taylor series for computing the function converges so slowly that the solution is of no practical use.

Examples given by the authors are analytical solutions to the N-body problem proposed by Sundman (1907, 1909) and Wang (1991):

There are analytical solutions to the N-body problem (Sundman 1907, 1909; Wang 1991), which are perfectly general and have been viewed as theoreti- cal successes. However, these solutions remain useless for physicists: neither Sundman’s solution nor Wang’s solution can be used to compute the trajecto- ries of an N-body system. This failure is intrinsically due to the form of these solutions, which are infinite convergent series in powers of t1/3, and thus have no practical use to make quantitative predictions. As Wang himself writes:

Although the conclusion given here provides a way to integrate the N- body problem, one does not obtain a useful solution in series expansion. The reason for this is because the speed of convergence of the resulting solution is terribly slow. One has to sum, for example, an incredible number of terms, even for an approximate solution of first order in q, p, t. (Wang 1991, p. 87)8

[...] As Florin Diacu asserts about Sundman’s solution and Wang’s solution:

One would have to sum up millions of terms to determine the motion of the particles for insignificantly short intervals of time. The round-off errors make these series unusable in numerical work. From the theoret- ical point of view, these solutions add nothing to what was previously known about the n-body problem. (1996, p. 70).

Another example given by the authors is the Airy function, for which we have an analytical solution whose power series converges analytically, even though it doesn't converge in a floating point calculation:

Infinite series are not always convenient for numerical application. The reason, as Fillion and Bangu claim, is that one cannot sum an infinite number of terms. In such a case, one can only sum the terms of a truncated version of the series. But such finite summation may lead to an inaccurate numerical result. Fillion and Bangu (2015) illustrate this with the Taylor series of the Airy function (a special function). They write “[n]umerically, even if the series converges for all x, it might be of little practical use, since the theoretical uniform conver- gence might not translate to success in numerical contexts”(p. 6). Ai(−t) is the analytical solution to d2x/dt2 = −tx. The equation describes an harmonic oscillator whose stiffness increases over time. Fillion and Corless (2014) insist that, although the infinite Taylor series of Airy function “converges uniformly, the floating-point computation diverges” (p. 1461). This means that a numer- ical computation of the Taylor series of the function – e.g. Ai(−12.82) – with floating-point arithmetic does not necessary tend to the exact value while the numerical computation includes more and more terms of the series. Such a numerical computation thus seems to be useless for the purpose of making quantitative predictions.

Answer 2

On the top of the comments, I'd like to add an example where your way may fail even for simple autonomous linear odes. Consider the following $$y'=Ay, y(0)= \mathbf{1}$$ with $A = \begin{bmatrix} -49 & 24 \\ -64 &31 \end{bmatrix}$

The exact solution is of course $y = e^A y_0$, and you may think to approximate $e^A$ with its definition $$e^A = I +A +\frac{A^2}{2} + \frac{A^3}{3!} + ... \frac{A^k}{k!}$$ and truncate the series at some proper $k$.

This is a fairly legitimate way from the theoretical viewpoint: however, what you're going to get is totally wrong. See the discussion about this at section 3 "Series methods" (the matrix $A$ is the same I wrote there) here on the seminal paper Nineteen dubious ways to compute the exponential of a matrix by Cleve Moler and Charles Van Loan