I'm a teacher assistant in an ODE course, and I was thinking about the different methods to solve a linear ODE of order $n$ with constant coefficients, and how long they can get. For example you need to invert a matrix (a very costly operation) and do a bunch of integrals if you want to use variation of parameters, or if you want to use undetermined coefficients, finding the coefficients can be a long process (I don't know the actual computational cost though).
I'd like to tell the students that you can't solve the equations faster than this, or that you can actually do it faster, at least theorically, so my question is: is there an upper/lower bound for the time complexity of analytically solving a linear ODE?
This answer is not about upper/lower bound for the time complexity of analytically solving a linear ODE. It focused on helping students or readers how to solve $n$ order ODE with starting from first order ODE and then to show how to get general solution and show some simple methods can reduce time to get the final result.
Firstly, students should understand how the general solution looks like and how they can reduce time to solve ODEs with some simple methods. I do not claim that it is upper/lower bound time for an ODE question. I recommend that they use these methods to solve any type of constant coefficient linear ODEs.
Firstly ,Let's solve first order constant coefficient linear ODE.
$$y^{'}(x)-r_1 y(x)=f(x)$$
General solution for it:
$$y(x)=e^{r_1x}\int f(x)e^{-r_1x}dx=e^{r_1x}(\int f(x)e^{-r_1x}dx+c)=ce^{r_1x}+e^{r_1x}\int f(x)e^{-r_1x}dx$$ $$y(x)=ce^{r_1x}+e^{r_1x}\int f(x)e^{-r_1x}dx$$ where $c$ is a constant
Note that we will use this solution for second order ODEs
$$y^{''}(x)+a_{1}y^{'}(x)+a_0 y(x)=f(x)$$
Need to find out roots of 2 order characteristic polynomial.
$r^{2}+a_{1}r+a_0=0$
Roots are $r_1,r_2$
$$(r-r_1)(r-r_2)=0$$ $$r^2-(r_1+r_2)r+r_1r_2=0$$
We can rewrite second order ODE as 2 separate first order ODEs:
$$y^{'}(x)-r_1 y(x)=z_1(x)$$ $$z_1^{'}(x)-r_2 z_1(x)=f(x)$$
If we put $z_1(x)$ in second ODE, we can get
$$(y^{'}(x)-r_1 y(x))'-r_2 (y^{'}(x)-r_1 y(x))=f(x)$$ $$y^{''}(x)-(r_1+r_2) y'(x)+r_1r_2 y(x)=f(x)$$
This is the same ODE with $$y^{''}(x)+a_{1}y^{'}(x)+a_0 y(x)=f(x)$$
We separately will solve 2 first order ODEs Firstly,
$$y^{'}(x)-r_1 y(x)=z_1(x)$$
$$y(x)=c_1e^{r_1x}+e^{r_1x}\int z_1(x)e^{-r_1x}dx$$
Secondly, $$z_1^{'}(x)-r_2 z_1(x)=f(x)$$ $$z_1(x)=c_2e^{r_2x}+e^{r_2x}\int f(x)e^{-r_2x}dx$$
$$y(x)=c_1e^{r_1x}+e^{r_1x}\int e^{-r_1x}(c_2e^{r_2x}+e^{r_2x}\int f(x)e^{-r_2x}dx)dx$$
$$y(x)=c_1e^{r_1x}+c_2e^{r_1x}\int e^{(r_2-r_1)x}dx+e^{r_1x}\int e^{(r_2-r_1)x}\int f(x)e^{-r_2x}dxdx$$
Now we will attempt to solve 3 order ODE:
$$y^{'''}(x)+a_{2}y^{''}(x)+a_{1}y^{'}(x)+a_0 y(x)=f(x)$$
$$y^{''}(x)-(r_1+r_2) y'(x)+r_1r_2 y(x)=z_1(x)$$
$$z_1^{'}(x)-r_3 z_1(x)=f(x)$$
$$y(x)=c_1e^{r_1x}+c_2e^{r_1x}\int e^{(r_2-r_1)x}dx+e^{r_1x}\int e^{(r_2-r_1)x}\int z_1(x)e^{-r_2x}dxdx$$
$$z_1(x)=c_3e^{r_3x}+e^{r_3x}\int f(x)e^{-r_3x}dx$$
$$y(x)=c_1e^{r_1x}+c_2e^{r_1x}\int e^{(r_2-r_1)x}dx+e^{r_1x}\int e^{(r_2-r_1)x}\int e^{-r_2x}(c_3e^{r_3x}+e^{r_3x}\int f(x)e^{-r_3x}dx)dxdx$$
$$y(x)=c_1e^{r_1x}+c_2e^{r_1x}\int e^{(r_2-r_1)x}dx+c_3e^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}dxdx+e^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}\int f(x)e^{-r_3x}dxdxdx$$
If we generalize the solution:
$$y^{(n)}+a_{n-1}y^{(n-1)}+a_{n-2}y^{(n-2)}+.....+a_{2}y^{''}+a_{1}y^{'}+a_0 y=f(x)$$
$r^n+a_{n-1}r^{n-1}+a_{n-2}r^{n-2}+.....+a_{2}r^{2}+a_{1}r+a_0=0$
We will get n roots for this characterictic equation ; If characteristic equation ($n>=5$), the roots of this equation cannot be expressed by radicals but Jerrard demonstrated that quintics ($n=5$) can be solved by using ultra radicals (also known as Bring radicals) . 1858 Charles Hermite showed that the Bring radical could be characterized in terms of the Jacobi theta functions and their associated elliptic modular functions. Higher degree polynomial equation roots can be solved by some special functions but I have not seen a general method to solve for all orders yet. Please see my wiki reference link in beyond radicals section . link.
Homogeneous solution : $$y_h(x)=c_1e^{r_1x}+c_2e^{r_1x}\int e^{(r_2-r_1)x}dx+c_3e^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}dxdx+c_4e^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}\int e^{(r_4-r_3)x}dxdxdx+c_5e^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}\int e^{(r_4-r_3)x}\int e^{(r_5-r_4)x}dxdxdxdx+.........+c_ne^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}\int e^{(r_4-r_3)x}\int.....\int e^{(r_{n}-r_{n-1})x}dxdx......dx \tag 1$$
Non-Homogeneous solution (particular solution) : $$y_p(x)=e^{r_1x}\int e^{(r_2-r_1)x}\int e^{(r_3-r_2)x}\int e^{(r_4-r_3)x}\int.....\int e^{(r_{n}-r_{n-1})x}\int f(x)e^{-r_nx}dxdxdx......dx$$
General solution is for n order linear ODE:
$$y(x)=y_h(x)+y_p(x)$$
First method to reduce time:Focus on repeated roots
We can reduce time while placing repeated roots in general solutions above if some roots are repeated
More generally, we can express the root of polynomial as: $$r^n+a_{n-1}r^{n-1}+a_{n-2}r^{n-2}+.....+a_{2}r^{2}+a_{1}r+a_0=0$$ $$(r-r_1)^{m_1}(r-r_2)^{m_2}...(r-r_k)^{m_k}=0$$ where $m_1+m_2+...+m_k=n$
$$y_h(x)=(c_{(1,1)}+c_{(1,2)} x+...+c_{(1,m_1)}x^{m_1-1})e^{r_1x}+(c_{(2,1)}+c_{(2,2)} x+...+c_{(2,m_2)}x^{m_2-1})e^{r_2x}+.....+(c_{(k,1)}+c_{(k,2)} x+...+c_{(k,m_k)}x^{m_k-1})e^{r_kx}$$
where $c_{(i,j)}$ are unknown constants. (Totally $n$ unknown constant must be.)
Non-homogenous solution:$y_p(x)$ We can reorder roots in any combination in solution of $y_p(x)$ and result will be the same because of symmetrical property of roots.
$$y_p(x)=e^{r_1x}\underbrace{\int.\int}_{m_1\text{ times}} e^{(r_2-r_1)x}\underbrace{\int.\int}_{m_2\text{ times}} e^{(r_3-r_2)x}\underbrace{\int.\int}_{m_3\text{ times}} e^{(r_4-r_3)x}.....\underbrace{\int.\int}_{m_{k-1}\text{ times}} e^{(r_{k}-r_{k-1})x}\underbrace{\int.\int}_{m_k\text{ times}} f(x)e^{-r_kx} \underbrace{dx..dx}_{n\text{ times}}$$
General solution is: $$y(x)=y_h(x)+y_p(x)$$ Second method to reduce time:Focus on complex conjugate roots
If the coefficients of characteristic polynomial are real constants and if any root has imaginary part, the characteristic polynomial must have a conjecture complex root too. $$r^n+a_{n-1}r^{n-1}+a_{n-2}r^{n-2}+.....+a_{2}r^{2}+a_{1}r+a_0=0$$
$$(r-r_1)^{m_1}(r-r_2)^{m_2}...(r-r_i)^{m_i}...(r-r_j)^{m_j}...(r-r_k)^{m_k}=0$$ where $m_1+m_2+...+m_i+...+m_j+...+m_k=n$
if ($r_i=\bar{r_j}$) then $m_i=m_j$ $$r_i=a+ib$$ $$r_j=a-ib$$
$m_1+m_2+...+2m_i+...+m_k=n$
$r_i=a+ib$
$$y_h(x)=(c_{(1,1)}+c_{(1,2)} x+...+c_{(1,m_1)}x^{m_1-1})e^{r_1x}+(c_{(2,1)}+c_{(2,2)} x+...+c_{(2,m_2)}x^{m_2-1})e^{r_2x}+.....+(c_{(i,1)}+c_{(i,2)} x+...+c_{(i,m_i)}x^{m_i-1})e^{ax}\cos(bx)+(c_{(j,1)}+c_{(j,2)} x+...+c_{(j,m_i)}x^{m_i-1})e^{ax}\sin(bx)...+(c_{(k,1)}+c_{(k,2)} x+...+c_{(k,m_k)}x^{m_k-1})e^{r_kx}$$