Short version: I am trying to implement a time response for analog transfer functions, based on Heun's method, to be implemented in C++, but I can't.
Long version (this may get long, so if you don't have the patience, skip this question):
Consider this test transfer function (TF):
$$H(s)=\frac{Y(s)}{X(s)}=\frac{a_2s^2+a_1s+a_0}{b_2s^2+b_1 s+b_0}=\frac{2.3s^2+0.1s+4.5}{1.3s^2+0.7s+1.1}$$
where the input, X(s), can be any signal, Dirac, Heaviside, sine, noise, etc, but sampled, because I need this to be implemented in C++; thus, it's a vector. In order to make a time reponse formula out of it, and to avoid derivatives, I had integrated the TF to give:
$$y(t)=a_2x(t)+\omega\int{\left[a_1x(t)-b_1 x(t)+\omega\int{\left(a_0 x(t)-b_0 y(t)\right)\text{d}t}\right]\text{d}t}$$
This was used to implement it with the basic forward Euler (T=sampling period, and from now on consider all $a_n$ and $b_n$ terms as scaled by $1/b_2$):
$$\begin{align}y''_{n}&=T\left(a_0 x_{n-2} - b_0 y_{n-2}\right) + y''_{n-1}\\ y'_{n}&=T\left(a_1 x_{n-1} - b_1 y_{n-1} + y''_{n}\right) + y'_{n-1}\\ y_n&=a_2 x_n + y'_{n}\end{align}$$
where the initial two samples were given by $y_0=a_2$ and $y_1=T(a_1x_0-b_1y_0)$. While it works, it suffers from the typical cumulative error of forward Euler, meaning I need a sampling time of at least 0.01 (relative to cutoff frequency) in order to have a fairly good accuracy. This also meant the resulting vector (in C++) was large for plotting.
After some searching, I found out about other methods, among them Heun's (HM), so I tried to convert what I had. The problems are as follows:
Since all I read about the method implied it being used as a means to solve ODEs, I gave up the integration in favour of the direct ODE that results from the TF. I tried solving it in terms of a generic input x(t), which, at its worst, is a Dirac function. This means that I would have this:
$$b_2y''(t)+b_1y'(t)+b_0y(t)=a_2\delta''(t)+a_1\delta'(t)+a_0\delta(t)$$
so two derivatives of Dirac. User @LutzL made it clear that I need that the right hand side be at least as smooth a function as the left side, for which he was kind enough to provide a solution, because I realized I need to calculate the initial conditions in order to go further.
But this means a continuous function, which I don't have, since the input is a vector. What's more, all the descriptions so far give solutions for continuous time ODEs, or some function to be approximated, which I don't have. But if I did, this is what I would have come up with:
$$z=y'\Rightarrow$$ $$b_2z'=a_2x''+a_1x'+a_0x-b_1z-b_0y$$
But what do I do with $x''$ and $x'$? Maybe transform them as (i.e.) $x'=w$, in order to keep the ODE as a 1st order, but then what do I do with $w$? Do I use it in the formula of $z=y'-x'$? Do I leave it in the formula for $z'$? If I follow the linked answer from @LutzL, I get the analog, continuous time function, which I don't have because the input is some unknown series of discrete samples.
Also, if I am to try to determine whether the input is steep, and, thus, risking discontinuity, how would I do that with discrete values? Because I approximate the Dirac (or Kronecker delta) by having the first sample as 1/T (so the area is 1). What's stopping me from having a cosine input? If it's low frequency, it will seem as a Heaviside function; if it's conveniently sampled at twice its frequency, I have a train of alternating signs Kronecker deltas; and even if I have a nice sine input, what if the amplitude of it is high enough that sin(T)=1/T? Will that count as discontinuity?
And yet, there has to be a way this is done since they work just fine in SPICE solvers, for example (I couldn't find examples of these). So I thought I'd give up ODEs in favour of, again, integration, which seemed to work before, but needed only adaptation to HM. I won't post everything I have tried, there was a lot of hammering in the process (lacking math skills provides you with a strong arm), but here's the last try, in wxMaxima code. I kept the first two samples as in forward Euler, since they seem quite close (<1%) to what LTspice shows (which I know and trust since I've been working with it for many years):
y[1] : A2*u[1],
z : A1*u[1] - B1*y[1],
ky1 : z*t,
kz1 : z + ky1,
ky2 : ky1*(1 + t),
kz2 : kz1*(1 + t),
y[2] : A2*u[2] + (ky1 + ky2)/2,
z : z + (kz1 + kz2)/2,
u is the input, ky and kz are the states and the estimated states. t is the sampling period. At this point, another confusion: the formula shows derivatives, that is it should be:
$$y'=a_2x'+a_1x-b_1y$$
but isn't $x'$ simply the current sample (which is known, true)? At any rate, other examples seem to treat it this way, so that's what I'll do, too. Next there's a loop:
ky1 : z,
kz1 : (A1*u[i-1] - B1*y[i-1] + (A0*u[i-2] - B0*y[i-2])*t)*t,
ky2 : ky1*(1 + t),
kz2 : (A1*u[i] - B1*y[i-1] + (A0*u[i-1] - B0*y[i-1])*t + kz1)*t,
y[i] : A2*u[i] + (ky1 + ky2)/2,
z : z + (kz1 + kz2)/2
The biggest problem is that I don't know how to transform the functions (i.e. $f_1(t,z,y)$) shown in all the examples into what I have, samples. Whatever came out was by ogling to the forward Euler. Whatever I tried, the results are either diverging (as if the poles on the right hand side), or come as ultra damped. But if I try with a sine input, it seems to converge to ..some sort of response, wrong, but it converges. Is this because the right side isn't smooth if it's a step, or impulse? At this point I simply don't understand how to do it and am confused a lot, and I would rather have someone provide the understanding for it, as opposed as simply replicating some steps like a robot.
How do I make the Heun's method work in code that works with a vector of samples as input? What do I use as the functions that are used in the examples for ODE of continuous time examples? But, mostly, why? It doesn't matter whether it turns out that the equations are based on the derivative, or on the integration method, as long as I understand it and it works.
You get a sampled input signal $x$, so that you either need to interpolate it to apply general Runge-Kutta methods that use values inside the interval, or apply multi-step methods that only use samples at the given sample points. Heun's method can be interpreted as predictor-corrector implementation of the implicit trapezoidal method $$ \vec u_{n+1}-\vec u_n=\frac{h}2(f(t_n,\vec u_n)+f(t_{n+1},\vec u_{n+1})), $$ so that it fits into both frameworks.
The next point is the derivatives on the right side. You could use divided differences, but that will amplify the error. The other variant is to compensate for the derivatives of $x$ in $$ b_2y''+b_1y'+b_0y=a_2x''+a_1x'+a_0x $$ by explicitly incorporating terms based on anti-derivatives of $x$ in a parametrization of the solution, as integration is usually more robust than differentiation. $$ y = u + c_2x+c_1X+c_0\Xi,~~\text{ where }~~ \Xi'=X,~X'=x $$ from where you get $$ y'=u' + c_2x'+c_1x+c_0X\\ y''=u'' + c_2x''+c_1x'+c_0x $$ so that in the differential equation $$ b_2u''+b_1u'+b_0u=a_2x''+a_1x'+a_0x-b_2c_2x''-(b_2c_1+b_1c_2)x'-(b_2c_0+b_1c_1+b_0c_2)x-(b_1c_0+b_0c_1)X-b_0c_0\Xi $$ the derivative terms are removed by setting $b_2c_2=a_2$, $b_2c_1+b_1c_2=a_1$, $b_2c_0+b_1c_1+b_0c_2=a_0$ and thus obtaining the first-order system \begin{align} X'&=x\\ \Xi'&=X\\ u'&=v\\ v'&=-\frac1{b_2}\left(b_1v+b_0u+(b_1c_0+b_0c_1)X+b_0c_0\Xi\right) \end{align} Now apply the Heun method or the implicit trapezoidal method (as the system is linear, this can be solved directly) and then reconstruct $y$ from $u$.