Discrete State Space Representation of a second order system A second order disturbance-based model can be represented as follows:
$$\ddot{y}=f(y,\dot{y},w,t)+bu(t)$$
where $f(y,\dot{y},w,t)$ is generalized disturbance. This can be represented in continuous time state as:
$$\dot{x}=Ax(t)+Bu(t)+E\dot{f}(t)$$
$$y(t)=Cx(t)$$
where, $A=
    \begin{bmatrix}
    0 & 1 & 0 \\
    0 & 0 & 1 \\
    0 & 0 & 0 \\
    \end{bmatrix}$, $B=\begin{bmatrix}0\\b\\0\end{bmatrix}$, $C=\begin{bmatrix}1 & 0 & 0\end{bmatrix}$, $E=\begin{bmatrix}0\\0\\1\end{bmatrix}$ and $x=\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}y\\\dot{y}\\f\end{bmatrix}$.
When the above state space equation is discretized, we can get the following:
$$x(k+1)=\Phi x(k)+\Gamma u(k)+E_df(k+1)$$
$$y(k)=Cx(k)$$
My question is,


*

*Is the term $f(k+1)$ for discretized $\dot{f}(t)$ correct or should it be $f(k)$?

*How do you determine the term $E_d$ associated with $\dot{f}(t)$ using zero order hold?

 A: In continuous time, you assume that $\dot{f}$ exists and can be written as a function $\xi(t):=\dot{f}(t)$. Then you can write your model as $$\dot{x}=Ax + Bu + E\xi$$ making the difference between the states and external inputs more explicit.
Following the same idea, in discrete time you have to assume the existence of a function $\Delta(k)$ such that $f(k+1)=f(k)+T_s\Delta(k)$, where $T_s$ is the sampling time. Then your model becomes $$x(k+1) = \Phi x(k) + \Gamma u(k) + E_d \Delta(k),$$ where the last line in $\Phi$ is $\begin{bmatrix}0 & 0 & 1\end{bmatrix}$ and $E_d = \begin{bmatrix}0 & 0 & T_s\end{bmatrix}^\top$.
Obviously, a good guess for $\Delta(k)$ is $$\Delta(k) = \frac{1}{T_s}\int_{t_k}^{t_{k+1}}\dot{f}(s)ds.$$
A: The analytical solution of the differential equation is as follows:
$$ x(t) = e^{At} x(t_0) + \int_{t_0}^t e^{A(t-\tau)} B u(\tau) d\tau + \int_{t_0}^t e^{A(t-\tau)} E \dot{f}(\tau) d\tau $$
To discretize this system we have to find the solution at $(k+1) T$ (assuming $t_0=0$). So,
$$\begin{align}
x((k+1) T) &= e^{A(k+1) T} x(0) + \int_{0}^{(k+1) T} e^{A((k+1) T-\tau)} B u(\tau) d\tau + \int_{0}^{(k+1) T} e^{A((k+1) T-\tau)} E \dot{f}(\tau) d\tau \\
&= e^{A(k+1) T} x(0) + \int_{0}^{kT} e^{A((k+1) T-\tau)} B u(\tau) d\tau + \int_{kT}^{(k+1) T} e^{A((k+1) T-\tau)} B u(\tau) d\tau + \\
&\phantom{=} \int_{0}^{kT} e^{A((k+1) T-\tau)} E \dot{f}(\tau) d\tau + \int_{kT}^{(k+1) T} e^{A((k+1) T-\tau)} E \dot{f}(\tau) d\tau \\
&= e^{AT} x(kT) + \int_{kT}^{(k+1) T} e^{A((k+1) T-\tau)} B u(\tau) d\tau + \int_{kT}^{(k+1) T} e^{A((k+1) T-\tau)} E \dot{f}(\tau) d\tau
\end{align}$$
Define $\eta = (k+1) T-\tau$. Then,
$$x((k+1) T) = e^{AT} x(kT) + \int_{0}^{T} e^{A\eta} B u((k+1) T-\eta) d\eta + \int_{0}^{T} e^{A\eta} E \dot{f}((k+1) T-\eta) d\eta $$
Assuming zero-order hold on $u(t)$ we can write $u((k+1) T-\eta) = u(kT)$ to obtain the well-known discretization result
$$x((k+1) T) = e^{AT} x(kT) + \left( \int_{0}^{T} e^{A\eta} d\eta \right) B u(kT)  + \int_{0}^{T} e^{A\eta} E \dot{f}((k+1) T-\eta) d\eta $$
Obviously, if you have further information on $\dot{f}$, you can go further. But taking the average of $\dot{f}$ gives you only an approximation, which gets worse when $T$ gets bigger.

Edit: I've misread the question. The above calculations is true only when $f$ is independent from the states. Otherwise one needs to linearize the system around an equilibrium point.
A: Your dynamics would be equivalent to
$$
\dot{x}(t) = A\,x(t) + B^*\,u^*(t), \tag{1}
$$
with
$$
B^* = 
\begin{bmatrix}
B & E
\end{bmatrix}, \quad
u^*(t) = 
\begin{bmatrix}
u(t) \\ \dot{f}(t)
\end{bmatrix}. \tag{2}
$$
This more standard expression for a continuous linear state space model has a known solution for discretization using zero order hold (so $u^*(t)=u^*[k]\,\forall\ k\,T\leq t<(k+1)T$ with $T$ the time step size), namely
$$
\begin{bmatrix}
A_d & B^*_d \\ 0 & I
\end{bmatrix} = 
e^{\begin{bmatrix}
A & B^* \\ 0 & 0
\end{bmatrix}T}, \tag{3}
$$
such that
$$
x[k+1] = A_d\,x[k] + B^*_d\,u^*[k]. \tag{4}
$$
Using this for your system yields
$$
A_d = 
\begin{bmatrix}
1 & T & \frac{1}{2}T^2 \\
0 & 1 & T \\
0 & 0 & 1
\end{bmatrix}, \quad
B^*_d = 
\begin{bmatrix}
\frac{b}{2}T^2 & \frac{1}{6}T^3 \\
b\,T & \frac{1}{2}T^2 \\
0 & T
\end{bmatrix}. \tag{5}
$$
If you want to express it using $f[k+1]$ instead of $\dot{f}[k]$ you can use that from zero order hold one gets
$$
f[k+1] = f[k] + T\,\dot{f}[k] \leftrightarrow \dot{f}[k] = \frac{f[k+1] - f[k]}{T}. \tag{6}
$$
Substituting $(6)$ into $(4)$ yields
$$
x[k+1] = \hat{A}\,x[k] + \hat{B}\,v[k], \tag{7}
$$
with
$$
v[k] = 
\begin{bmatrix}
u[k] \\ f[k+1]
\end{bmatrix}, \quad
\hat{A} = 
\begin{bmatrix}
1 & T & \frac{1}{3}T^2 \\
0 & 1 & \frac{1}{2}T \\
0 & 0 & 0
\end{bmatrix}, \quad
\hat{B} = 
\begin{bmatrix}
\frac{b}{2}T^2 & \frac{1}{6}T^2 \\
b\,T & \frac{1}{2}T \\
0 & 1
\end{bmatrix}. \tag{8}
$$
