Using scipy fft and ifft to solve ordinary differential equation numerically

Question

I have an `ordinary differential equation` in time domain as follows: ``` C*du/dt = -g*u + I ``` where `I = A*t/tau*exp^(1-t/tau)` in the freq domain: ``` u(w) = I(w)/(g*(1+C/g*j*w)) ``` `j` being the complex number `sqrt(-1)` hence i can get `u(t)` by going into the freq domain using fast Fourier transform (`fft`) and then back using `ifft`. the codes: ``` t = np.linspace(0.,499.9,5000) I = q*t*np.exp(1-t/Tau_ca)/Tau_ca u1 = np.fft.ifft(np.fft.fft(I)/(G_l*(1.+1.j*(C_m/G_l)*np.fft.fftfreq(t.shape[-1])))) ``` However, when I compare this `u(t)` obtained with other methods such as numerical integration of the differential equation or its analytical form, it is not correct. I have tried and been unsuccessful at figuring out where my mistakes are. please enlighten.

Accepted Answer

The derivative of a sinusoid, or complex exponential, is directly proportional to its frequency, and phase shifted by `π/2`. For a complex exponential the phase shift is equivalent to multiplying by `j`. For example, `d/dt exp(j*Ω*t)` `== j*Ω * exp(j*Ω*t)` `== Ω * exp(j*π/2) * exp(j*Ω*t)` `== Ω * exp(j*(Ω*t + π/2))`. So if you have the Fourier transform pair `u(t) <=> U(Ω)`, then `du/dt <=> jΩ * U(Ω)`. Integration is pretty much the inverse operation, but it may add an impulse if there's a DC component: `∫udt <=> U(Ω) / jΩ + π*U(0)*δ(Ω)` To approximate the derivative using a sampled sequence, you have to scale the discrete-time angular frequency ω (which ranges from 0 to 2π, or -π to π) by the sample rate, `fs`. For example, say the discrete-time frequency is π/2 radians/sample, such as the sequence `[1, 0, -1, 0, ...]`. In the original signal this corresponds to `fs/4`. The derivative is `d/dt cos(2*π*fs/4*t) == d/dt cos(π/2*fs*t)` `== -π/2*fs*sin(π/2*fs*t)` `== π/2*fs*cos(π/2*fs*t + π/2)`. You have to sample at an `fs` that's more than twice the bandwidth of the signal. The sampled component at exactly `fs/2` is unreliable. Specifically, with only 2 samples per cycle the amplitude of the `fs/2` component alternates the sign of the first sample in the sequence. So for a real-valued signal the `fs/2` DFT component is real valued, with a phase of either 0 or π radians, i.e. `a*cos(π*fs*t + {0, π})`. Given the latter, the derivative of the `fs/2` component is `-a*π*fs*cos(π*fs*t + {π/2, 3*π/2})`, which is 0 for the sample times `t == n/fs`. (Regarding the latter, the standard trigonometric interpolation of the DFT uses a cosine, and in that case the derivative will be zero at the sample points. The isn't necessarily true if you sample the signal and its derivative simultaneously. Sampling loses the phase of the `fs/2` component relative to the signal, but not relative to its derivative. Depending on the time you start sampling, both the `fs/2` component and its derivative may be non-zero at the sample points. If by luck one of them is 0 at the sample times, the other will be at its peak, since they're `π/2` radians apart.) Given that the `fs/2` DFT component is always real valued for a real-valued signal, when you multiply it by `j` in the course of computing the derivative or integral, this introduces an imaginary-valued component in the result. There's a simple workaround. If `N` is even, just zero out the `fs/2` component at index `N/2`. Another problem is division by zero when dividing by `jΩ` for integration. This can be solved by adding a small value to index 0 of the `Ω` vector (e.g. `finfo(float64).tiny` is the smallest double precision float). Given `Ω = fs * ω`, the system shown in the question has the following form in the frequency-domain: ``` H(Ω) = 1 / (g + j*Ω*C) U(Ω) = I(Ω) * H(Ω) ``` It's a single-pole low-pass filter. The solution you derived has 2 problems. - You aren't scaling the frequency variable `w` by `fs`. - You use `fftfreq`, which uses the range -0.5 to 0.5. You need -π to π. Actually you only need 0 to π because `i(t)` is real-valued. In this case you can use `rfft` and `irfft` for real-valued signals, which skips computing the negative frequencies. All that said, you may still be disappointed with the result because the DFT uses the periodic extension of your signal. Examples First, here's a simple example of a 1 Hz sinusoid (plotted in red) sampled at 1024 samples/s for 2 seconds, and its derivative computed via the DFT (plotted in blue): ``` from pylab import * fs = 1024 t = arange(2*fs, dtype=float64) / fs N = len(t) // 2 + 1 # non-negative frequencies w = 2 * pi / len(t) * arange(N) Omega = w * fs x0 = cos(2*pi*t) # w0 == 2*pi/fs X0 = rfft(x0); # Zero the fs/2 component. It's zero here anyway. X0[-1] = 0 dx0 = irfft(X0 * 1j*Omega) plot(t, x0, 'r', t, dx0, 'b') show() ``` This is an easy case -- a periodic signal with finite bandwidth. Just make sure to sample an integer number of periods at a high enough rate to avoid aliasing. The next example is a triangle wave, with a slope of 1 and -1, and a discontinuity in the derivative at the center and edges. Ideally, the derivative should be a square wave, but computing that perfectly would require infinite bandwidth. Instead there's Gibbs ringing around the discontinuity: ``` t2 = t[:fs] m = len(t) // (2*len(t2)) x1 = hstack([t2, 1.0 - t2] * m) X1 = rfft(x1) X1[-1] = 0 dx1 = irfft(X1 * 1j*Omega) plot(t, x1, 'r', t, dx1, 'b') show() ``` The DFT's implicit periodic extension is problematic if you're solving a non-periodic system. Here's a solution to the system in question using both `odeint` and the DFT (`tau` is set to 0.5s; `g` and `C` are set for a 1 Hz corner frequency): ``` from scipy.integrate import odeint a = 1.0; tau = 0.5 g = 1.0; C = 1.0 / (2 * pi) i = lambda t: a * t / tau * exp(1 - t / tau) f = lambda u, t: (-g*u + i(t)) / C H = 1 / (g + 1j*Omega*C) # system function I = rfft(i(t)) I[-1]…

Using scipy fft and ifft to solve ordinary differential equation numerically

Solution

Problem