# More efficient way of solving for transfer function between boundaries 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 $$k(x)$$, $$m(x)$$ the known material properties and the boundary conditions

$$u(t,-L) = f(t), \quad \frac{\partial u}{\partial t}(t,L) = -\sqrt{\frac{k(L)}{m(L)}} \frac{\partial u}{\partial x}(t,L). \tag{2}$$

I am interested in calculating the transfer function from $$f(t)$$ to $$u(t,L)$$. In order to do this 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)$$ to $$(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{3}$$

The separation of variables together with the first boundary condition from $$(2)$$ 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 $$(3)$$ as

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

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

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

The second order ordinary differential equation from $$(4)$$ together with the constraints $$X(-L) = 1$$ and $$(5)$$ 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{6}

It can be noted that $$z(L)$$ is indirectly also a function of $$s$$ and the transfer function I am after 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{7}$$

Linear time varying systems in general don't have a closed form analytical solution. Therefore, I think I have to resort to a numerical solution in order to solve $$(6)$$. 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{8}$$

and thus the value for the transfer function from $$(7)$$ for a given $$s$$ can be obtained by solving $$(8)$$ 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)$$ should be very small, in which case the obtained numerical solution might be very inaccurate. Instead, in such case it might be more accurate to solve $$(6)$$ in the reverse $$x$$ direction. For example by using $$\hat{z}(x) = z(-x)$$, such that the dynamics from $$(6)$$ 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{9}$$

with $$\hat{z}(-L) = z(L)$$, thus the calculated state transition matrix would be $$\Phi(-L,L)$$ and system of linear equations from $$(8)$$ 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{10}$$

I implemented this in Matlab using the following code:

L = 1;
B = 3 * pi / L;
r = 1.001;
K = @(x) 1 ./ (r + exp(-9 * cos(B * x))) + 1 - 1 / r;
M = @(x) 1 ./ (r + exp(-9 * cos(B * x))) + 1 - 1 / r;
Kx = @(x) -a * B * sin(B * x) .* exp(-9 * cos(B * x)) .* (K(x) - 1 + 1 / r).^2;
N = 1e3;
S = logspace(-2, 2, N);
G = zeros(1, N);
opts_low = odeset('RelTol',1e-2,'AbsTol',1e-3);
opts_high = odeset('RelTol',1e-6,'AbsTol',1e-8);
b = [1; 0];
for n = 1 : N
s = 1i * S(n);
f_back = @(x, z) [0 -1; -s^2*M(-x)/K(-x) Kx(-x)/K(-x)] * z;
[~, z1] = ode45(f_back, [-L L], [1; 0], opts_low);
[~, z2] = ode45(f_back, [-L L], [0; 1], opts_low);
Phi = [z1(end,:)' z2(end,:)'];
A = [b'*Phi; [s*sqrt(M(L)/K(L)) 1]];
if abs(b' * (A \ b)) < 1
[~, z1] = ode45(f_back, [-L L], [1; 0], opts_high);
[~, z2] = ode45(f_back, [-L L], [0; 1], opts_high);
Phi = [z1(end,:)' z2(end,:)'];
A = [b'*Phi; [s*sqrt(M(L)/K(L)) 1]];
G(n) = b' * (A \ b);
else
f_frwd = @(x, z) [0 1; s^2*M(x)/K(x) -Kx(x)/K(x)] * z;
[~, z1] = ode45(f_frwd, [-L L], [1; 0], opts_high);
[~, z2] = ode45(f_frwd, [-L L], [0; 1], opts_high);
Phi = [z1(end,:)' z2(end,:)'];
A = [b'; [s*sqrt(M(L)/K(L)) 1]*Phi];
G(n) = b' * Phi * (A \ b);
end
end
figure
bode(frd(G, S))


However, this approach seems a bit cumbersome to me. So I wonder if maybe there are other (more efficient) ways of calculating the transfer function $$G(s)$$?