0
$\begingroup$

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)$?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.