I am not very familiar with control theory but have a decent bit of experience with classical numerical integration.
I am looking at a the equation
$$\dot{x}(t) = Ax(t) + Bu(t) \hspace{1cm} x(0) =0 $$ (coming from a state-space model) and the way that I would handle this is by multiplying by the matrix exponential and integrating to get:
$$x(t) = \int_0^t e^{A(t - \tau)}Bu(\tau) \,d\tau$$
and then handling the convolution integral on the right hand side using some standard quadrature. And if I wanted to convert this continuous system to a discrete one I would use something like:
$$x(t + h) = x(t) + \int_t^{t + h} e^{A(t + h - \tau)} Bu(\tau) \,d\tau$$
with a quadrature rule for the integral on the RHS (which is of course just the quadrature rule of the full interval applied to a uniform partition.)
However, the authors of a paper I'm looking at use the bilinear transform from control theory.
My understanding of the bilinear transform is that it computes the transfer function of a discrete system coming from the trapezoid rule and the transfer function of the continuous system coming from the continuous integral and equates them to give a transformation to apply in "Laplace transform space" to go from continuous to discrete time. Is this correct?
Now, If I take the bilinear transform as a given and apply a Laplace transform to the above system, I understand that after a bit of algebra I will arrive at the discrete Laplace transform (or z-transform?) of the discrete system: $$x(t + h) = (I - h/2 A)^{-1}(I + h/2 A)x(t) + h(I - h/2 A)^{-1}Bu(t)$$
which I am not sure how to derive from the original system using a traditional numerical method (I tried just applying trapezoid rule + pade approximation but don't see immediately how to get the $$h(I - h/2 A)^{-1}B$$ portion).
My questions are:
How does this compare to other numerical methods? I can compute the order numerically so I'm not so concerned about that, but why would one take the approach of using the bilinear transfer function over the type of numerical integral/differential equation solve that I outlined initially?
The main reasons why we use (forward/backward) Euler or bilinear transformations is that they allow us to obtain an explicit form for the discretized model associated with a continuous-time system.
In control, the input is usually not known beforehand and must be computed based on the model of the system and the different control objectives. In this regard, it is not really possible to use numerical integration methods, as you first mentioned.
It could be possible to assume that the input is constant over each component of the partitions and obtain a model this way. In such a case, though, we have an explicit solution for a system whose input is piecewise constant and is given by
$$\displaystyle x(t+h)=e^{Ah}x(t)+\left(\int_0^he^{A(h-s)}Bds\right) u_t$$ where $u(t+s)=u_t$ for all $s\in[0,h)$.
Regarding the bilinear transform, it is enough to note that
$$x(t+h)-x(t)=\int_0^h(Ax(t+s)+Bu(t+s))ds.$$
The trapezoidal rule says that
$$\int_0^hf(s)ds\approx \dfrac{h}{2}(f(0)+f(h)),$$ which yields
$$\int_0^hAx(t+s)+Bu(t+s))ds\approx \dfrac{h}{2}(Ax(t)+Bu(t)+Ax(t+h)+Bu(t+h)).$$
Now define the shift operator $\sigma$ as $\sigma u(t)=u(t+h)$, which yields $$\dfrac{h}{2}(Ax(t)+Bu(t)+Ax(t+h)+Bu(t+h))\approx \dfrac{h}{2}(\sigma+1)(Ax(t)+Bu(t))$$
Therefore,
$$(\sigma -1)x(t)\approx\dfrac{h}{2}(\sigma+1)(Ax(t)+Bu(t))$$
and thus
$$\dfrac{h}{2}\dfrac{\sigma-1}{\sigma+1}x(t)\approx Ax(t)+Bu(t).$$
This shows that the derivative operator (i.e. $s$) can be approximated by the operator $\dfrac{h}{2}\dfrac{\sigma-1}{\sigma+1}$