# How to obtain the transfer function between boundary condition and point of wave equation?

I am considering the wave equation with position varying material properties

$$m(x) \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x}\left(k(x) \frac{\partial u}{\partial x}\right), \tag{1}$$

with $$u(t,x)\in\mathbb{R}$$ for $$t\geq 0$$, $$x\in[-L,L]$$, $$k(x)$$ the known local stiffness, $$m(x)$$ the known local density, initial conditions

\begin{align} u(0,x) &= g(x), \tag{2a}\label{ic1} \\ \frac{\partial u}{\partial t}(0,x) &= h(x), \tag{2b}\label{ic2} \end{align}

and the boundary conditions

\begin{align} u(t,-L) &= f(t), \tag{3a}\label{bc1} \\ \frac{\partial u}{\partial t}(t,L) &= -\sqrt{\frac{k(L)}{m(L)}} \frac{\partial u}{\partial x}(t,L), \tag{3b}\label{bc2} \end{align}

with $$f(t)$$ a known function of time (potentially free to be chosen by the user), such that $$\eqref{bc1}$$ is a Dirichlet boundary condition at $$x=-L$$, and $$\eqref{bc2}$$ at $$x=L$$ a non-reflecting boundary condition (the wave only moves to the right at the boundary).

I am interested in calculating the transfer function (or the frequency response function) from $$f(t)$$ to $$u(t,L)$$, thus

$$G(s) = \frac{\mathcal{L}\{u(t,L)\}(s)}{\mathcal{L}\{f(t)\}(s)},$$

with $$\mathcal{L}\{f(t)\}(s)$$ denoting the Laplace transform of $$f(t)$$. I am interested in this in order to learn about what impact spatially varying material properties has on the transmission of given frequencies.

How could one go about solving for this transfer function numerically?

In order to calculate the transfer function one can assume that $$f(t) = e^{s\,t}$$, with $$s$$ representing the Laplace variable. Using separation of variables $$u(t,x) = T(t)\,X(x)$$ in $$(1)$$ yields

$$\frac{\ddot{T}(t)}{T(t)} = \frac{\dot{k}(x)\,\dot{X}(x) + k(x)\,\ddot{X}(x)}{m(x)\,X(x)} = c. \tag{4}\label{sov_eq1}$$

The separation of variables together with the first boundary condition $$(3a)$$ can be satisfied using $$T(t) = e^{s\,t}$$ and $$X(-L) = 1$$. This also implies $$c = s^2$$, which can be used to write the right hand side of $$\eqref{sov_eq1}$$ as

$$\ddot{X}(x) = \frac{s^2\,m(x)\,X(x) - \dot{k}(x)\,\dot{X}(x)}{k(x)}. \tag{5}\label{ode_x}$$

By using the separation of variables with $$T(t) = e^{s\,t}$$, the second boundary condition $$(3b)$$ can also be written as

$$s\,X(L) = -\sqrt{\frac{k(L)}{m(L)}} \dot{X}(L). \tag{6}\label{ode_x_bc}$$

The second order ordinary differential equation from $$\eqref{ode_x}$$ together with the constraints $$X(-L) = 1$$ and $$\eqref{ode_x_bc}$$ is equivalent to the following linear "time" varying autonomous system

\left\{ \begin{align} \dot{z}(x) &= \begin{bmatrix} 0 & 1 \\ \frac{s^2\,m(x)}{k(x)} & \frac{-\dot{k}(x)}{k(x)} \end{bmatrix} z(x), \\ \text{s.t.} &\ \begin{bmatrix}1 & 0\end{bmatrix} z(-L) = 1, \quad \begin{bmatrix}s & \sqrt{\frac{k(L)}{m(L)}}\end{bmatrix} z(L) = 0, \end{align}\right. \tag{7}\label{LTV}

with $$z(x) = \begin{bmatrix}X(x) & \dot{X}(x)\end{bmatrix}^\top$$. It can be noted that $$z(L)$$ is indirectly also a function of $$s$$ and the transfer function can be calculated with

$$G(s) = \frac{\mathcal{L}\{u(t,L)\}(s)}{\mathcal{L}\{f(t)\}(s)} = \{\begin{bmatrix}1 & 0\end{bmatrix} z(L)\}(s). \tag{8}\label{tf}$$

Linear time varying systems in general don't have a closed form analytical solution. Therefore, one has to resort to a numerical solution in order to solve $$\eqref{LTV}$$. However, that system is linear, which means that the final state $$z(L)$$ is linear in the initial state $$z(-L)$$ according to the state transition matrix $$\Phi(L,-L)$$, such that $$z(L) = \Phi(L,-L)\,z(-L)$$. The two columns of the state transition matrix can be calculated numerically by using $$\begin{bmatrix}1 & 0\end{bmatrix}^\top$$ and $$\begin{bmatrix}0 & 1\end{bmatrix}^\top$$ as initial condition for $$z(-L)$$ respectively. Once an approximation of $$\Phi(L,-L)$$ is obtained, the initial condition condition can be calculated by solving the following system of linear equations for $$z(-L)$$

$$\begin{bmatrix} \begin{bmatrix}1 & 0\end{bmatrix} \\ \begin{bmatrix}s & \sqrt{\frac{k(L)}{m(L)}}\end{bmatrix} \Phi(L,-L) \end{bmatrix} z(-L) = \begin{bmatrix}1 \\ 0\end{bmatrix}, \tag{9}\label{LTV_sol}$$

and thus the value for the transfer function from $$\eqref{tf}$$ for a given $$s$$ can be obtained by solving $$\eqref{LTV_sol}$$ for $$z(-L)$$ and use $$z(L) = \Phi(L,-L)\,z(-L)$$.

When there is a lot of attenuation (so when $$|G(s)| \ll 1$$) the obtained value for $$z(L)$$ would be very small, in which case the obtained numerical solution might be badly conditioned. Instead, in such case it might be more accurate to solve $$\eqref{LTV}$$ in the reverse $$x$$ direction. For example by using $$\hat{z}(x) = z(-x)$$, such that the dynamics from $$\eqref{LTV}$$ can be written as

$$\dot{\hat{z}}(x) = \begin{bmatrix} 0 & -1 \\ \frac{-s^2\,m(-x)}{k(-x)} & \frac{\dot{k}(-x)}{k(-x)} \end{bmatrix} \hat{z}(x), \tag{10}$$

with $$\hat{z}(-L) = z(L)$$, thus the calculated state transition matrix would be $$\Phi(-L,L)$$ and system of linear equations from $$\eqref{LTV_sol}$$ would become

$$\begin{bmatrix} \begin{bmatrix}1 & 0\end{bmatrix} \Phi(-L,L) \\ \begin{bmatrix}s & \sqrt{\frac{k(L)}{m(L)}}\end{bmatrix} \end{bmatrix} z(L) = \begin{bmatrix}1 \\ 0\end{bmatrix}. \tag{11}$$

This can be implemented in Julia using the following code:

using DifferentialEquations, ControlSystems, ControlSystemIdentification, Plots

function main(N)
Omega = exp.(range(log(1e-2), log(1e1), N))
G = [eval_TF(1im * omega) for omega = Omega]
bodeplot(FRD(Omega, G))
end

function eval_TF(s)
L = 1.0
tspan = (-L, L)
X0 = [1.0 0.0; 0.0 1.0]
p = s^2
b = [1.0 0.0]
MdivK, KxdivK = material_parameters(L)

prob1_back = ODEProblem(dx_backward!, X0[:,1], tspan, p)
prob2_back = ODEProblem(dx_backward!, X0[:,2], tspan, p)
sol1 = solve(prob1_back, Tsit5(); reltol = 1e-3, abstol = 1e-4, saveat = tspan[2])
sol2 = solve(prob2_back, Tsit5(); reltol = 1e-3, abstol = 1e-4, saveat = tspan[2])
Phi = [sol1.u[end] sol2.u[end]]
A = [b * Phi; [s*sqrt(MdivK) 1]]
if abs.(b * (A \ b'))[1] < 1
sol1 = solve(prob1_back, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
sol2 = solve(prob2_back, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
Phi = [sol1.u[end] sol2.u[end]]
A = [b * Phi; [s*sqrt(MdivK) 1]]
return (b * (A \ b'))[1]
else
prob1_forw = ODEProblem(dx_forward!, X0[:,1], tspan, p)
prob2_forw = ODEProblem(dx_forward!, X0[:,2], tspan, p)
sol1 = solve(prob1_forw, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
sol2 = solve(prob2_forw, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
Phi = [sol1.u[end] sol2.u[end]]
A = [b; [s*sqrt(MdivK) 1]*Phi]
return (b * Phi * (A \ b'))[1]
end
end

function material_parameters(x)
L = 1
B = 5 * pi / L
MdivK = 1 # Constant wave speed
KxdivK = -90 * sin(B * x)^9 # Varying impedance
return MdivK, KxdivK
end

function dx_forward!(dx,x,p,t)
MdivK, KxdivK = material_parameters(t)
dx[1] = x[2]
dx[2] = p[1] * MdivK * x[1] - KxdivK * x[2]
end

function dx_backward!(dx,x,p,t)
MdivK, KxdivK = material_parameters(-t)
dx[1] = -x[2]
dx[2] = -p[1] * MdivK * x[1] + KxdivK * x[2]
end


Which for example can generate the following Bode plot: