The original results
Running the code you provided and comparing the result with the RK4 code posted in your other question, we get the following responses:
where the blue curve represents the FFT based implementation, and the red curve represents your alternate RK4 implementation. As you have pointed out, the curves are quite different.
Getting the correct response
The most obvious problem is of course the amplitude, and the main sources of the amplitude discrepancies of the code posted in this question are the same as the ones I indicated in my answer to your other question:
- The RK4 implementation performs a numeric integration which correctly scales the summed values by the integration step
dt. This scaling is lacking from the FFT based implementation.
- The impulse response used in the FFT based implementation is consistent with the driving force being scaled by the mass
m, a factor which was missing from the RK4 implementation.
Fixing those two issues results in responses which are a little closer, but still not identical. As you probably found out given the changes in the posted code of your other question, another thing that was lacking was zero padding of the input and of the impulse response, without which you were getting a circular convolution rather than a linear convolution:
f=[cos(t),zeros(1,length(t)-1)]; %force f
h=[1/m/wd*exp(-z*wn*t).*sin(wd*t),zeros(1,length(t)-1)]; %impulse response
Finally the last element to ensure the convolution yields good result is to use a good approximation to the infinite length impulse response. How long is long enough depends on the rate of decay of the impulse response. With the parameters you provided, the impulse response would have died down to 1% of its original values after approximately 11*pi. So extending the time span to tMax=14*pi (to include a full 2*pi cycle after the impulse response has died down), would probably be enough.
Obtaining the steady-state response
The simplest way to obtain the steady-state response is then to discard the initial transient. In this process we discard a integer number of cycles of the reference driving force (this of course requires knowledge of the driving force's fundamental frequency):
T0 = tMax-2*pi;
delay = min(find(t>T0));
sol = sol(delay:end);
plot([0:length(sol)-1]*dt, sol, 'b');
axis([0 2*pi]);
The resulting responses are then:
where again the blue curve represents the FFT based implementation, and the red curve represents your alternate RK4 implementation. Much better!
An alternate method
Computing the response for many cycles waiting for the transient response to die down and extracting the remaining samples corresponding
to the steady state might appear to be a little wasteful, despite the fact that the computation is still fairly fast thanks to the FFT.
So, let's go back a little and look at the problem domain. As you are probably aware,
the mass-spring-damper system is governed by the differential equation:
where f(t) is the driving force in this case.
Note that the general solution to the homogeneous equation has the form:
The key is then to realize that the general solution in the case where c>0 and m>0 vanishes in the steady state (t going to infinity).
The steady-state solution is thus only dependent on the particular solution to the non-homogenous equation.
This particular solution can be found by the method of undetermined coefficients, for a driving force of the form
by correspondingly assuming that the solution has the form
Substituting in the differential equation yields the equations:
thus, the solution can be implemented as:
EF0 = [wn*wn-w*w 2*z*wn*w; -2*z*wn*w wn*wn-w*w]\[1/m; 0];
sol = EF0(1)*cos(w*t)+EF0(2)*sin(w*t);
plot(t, sol);
where w=2*pi in your case.
Generalization
The above approach can be generalized for more arbitrary periodic driving forces by expressing the driving force as a
Fourier Series (assuming the driving force function satisfies the Dirichlet conditions):
The particular solution can correspondingly be assumed to have the form
Solving for the particular solution can be done in a way very similar to the earlier case. This result in the following implementation:
% normalize
y = F/m;
% compute coefficients proportional to the Fourier series coefficients
Yw = fft(y);
% setup the equations to solve the particular solution of the differential equation
% by the method of undetermined coefficients
k = [0:N/2];
w = 2*pi*k/T;
A = wn*wn-w.*w;
B = 2*z*wn*w;
% solve the equation [A B;-B A][real(Xw); imag(Xw)] = [real(Yw); imag(Yw)] equation
% Note that solution can be obtained by writing [A B;-B A] as a scaling + rotation
% of a 2D vector, which we solve using complex number algebra
C = sqrt(A.*A+B.*B);
theta = acos(A./C);
Ywp = exp(j*theta)./C.*Yw([1:N/2+1]);
% build a hermitian-symmetric spectrum
Xw = [Ywp conj(fliplr(Ywp(2:end-1)))];
% bring back to time-domain (function synthesis from Fourier Series coefficients)
x = ifft(Xw);
A final note
I purposely avoided the undamped c=0 case in the above derivation. In this case, the oscillation never die down and the general solution to the homogeneous equation does not have to be the trivial one.
The final "steady-state" in this case may or may not have the same period as the driving force. As a matter of fact, it may not be periodic at all if the period oscillations from the general solution is not related to the period of the driving force through a rational number (ratio of integers).
