We are considering the spring model
\begin{array}{l}
\frac{d^2 z}{d t^2} + C \frac{d z}{d t} + K z = 0, \\[10pt]
z(0) = 2 \quad , \quad \frac{dz}{dt}(0) = -C
\end{array}
\)
with parameters \(\Theta = [K,C]\) (here \(K\) denotes stiffness due to Hooke’s Law and \(C\) denotes a damping action; the mass \(m\) of the object is seen to be one from the presentation of the differential equation and initial conditions).
Our goal is to use the complex-step approximation to compute the sensitivities \(\frac{\partial z(t)}{\partial K}\) and \(\frac{\partial z(t)}{\partial C}\) for \(K=20.5\) and \(C = 1.5\). We also employ finite-difference methods, and compare them to the analytic values:
\begin{array}{l}
\frac{\partial z}{\partial K} = e^{-Ct/2} \frac{-2t}{\sqrt{4K – C^2}} \sin\left(\sqrt{K-C^2/4} \cdot t \right) , \\[15pt]
\frac{\partial z}{\partial C} = e^{-Ct/2} \left[ \frac{Ct}{\sqrt{4K-C^2}} \sin \left(\sqrt{K-C^2/4} \cdot t \right) – t \cos\left(\sqrt{K-C^2/4} \cdot t \right) \right]
\end{array}
\)
In general, complex-step analysis is performed by computing \(f'(t) \approx \frac{\text{Im}\left[f(t+ih)\right]}{h}\) for small \(h\). This is similar to the finite-difference approximation \(f'(t) \approx \frac{f(t+h)-f(t)}{h}\), though when dealing with analytic functions, complex-step is often preferred. We can use these methods to approximate our sensitivities. In doing so, we allow time to vary from 0 to 5 and choose a step-size of \(h=10^{-10}\).
For the first parameter \(K\), our approximations match the analytic solution nicely. For the second parameter \(C\), our numerical approximations agree, but do not do a great job approximating our solution. One potential explanation is that \(\frac{\partial z}{\partial C}\) has an extra trigonometric function, so may be harder to approximate.
%%%1a. Specify Parameters And Initial Conditions%%%
k = 20.5; %Stifness%
c = 1.5; %Damping%
z0 = 2; %IC 1%
dz0 = -c; %IC 2%
t = linspace(0,5,1000); %generate 1k/5=200 pts times every s%
%%%1b. Specify Step Sizes%%%
%Complex Step: f' \approx \frac{Im(f(x+ih))}{h}%
%Finite Difference: f' \approx \frac{f(x+h)}{h}
h_cs = 10^(-10); %Complex step step size%
h_fd = 10^(-10); %Finite difference step size%
%%%1c. Solve ODE%%%
function z = z_solution(t, k, c, z0, dz0)
omega = sqrt(k - (c^2)/4);
A = z0;
B = (dz0 + (c/2) * z0) / omega;
z = exp(-c*t/2) .* (A * cos(omega*t) + B * sin(omega*t));
end
spring_model = @(t, k, c) z_solution(t, k, c, z0, dz0);
%%%1d. Complex Step Approximation%%%
dzdk_cs = imag(spring_model(t, k + h_cs*1i, c)) / h_cs;
dzdc_cs = imag(spring_model(t, k, c + h_cs*1i)) / h_cs;
%%%1e. Finite Difference Approximation%%%
dzdk_fd = (spring_model(t, k + h_fd, c) - spring_model(t, k, c)) / (h_fd);
dzdc_fd = (spring_model(t, k, c + h_fd) - spring_model(t, k, c)) / (h_fd);
%%%1f. Analytic Solution (Given)%%%
dzdk_true = @(t, k, c) exp(-c*t/2) .* ...
(-2*t ./ sqrt(4*k - c^2)) .*...
sin(sqrt(k - c^2/4) * t);
dzdt_true = @(t, k, c) exp(-c*t/2) .* ...
(t .* cos(sqrt(k - c^2/4) * t) - c*t ./ sqrt(4*k - c^2)) .* ...
sin(sqrt(k - c^2/4) * t);
%%%1g. Plot Results%%%
figure;
plot(t, dzdk_cs, 'r-', t, dzdk_fd, 'b--', t, dzdk_true(t, k, c), 'g-.');
legend('Complex Step', 'Finite Difference', 'Analytic');
xlabel('Time'); ylabel('d/dk z(t)');
title('Sensitivity Of z(t) To k');
figure;
plot(t, dzdc_cs, 'r-', t, dzdc_fd, 'b--', t, dzdk_true(t, k, c), 'green');
legend('Complex Step', 'Finite Difference', 'Analytic');
xlabel('Time'); ylabel('d/dc z(t)');
title('Sensitivity Of z(t) To c');