[Math] Discretization of differential equation via FFT routine

discrete optimizationfourier analysisMATLABnumerical methodsordinary differential equations

I just have a question related to the following problem:

Find a discrete approximation to the differential equation $u^{\prime \prime} + 2u^{\prime} + 2u = 3\cos(6t)$ using Equation 3.12 for these values of $n$: $n = 16, n = 64$, and $n = 256$. You will need to use MATLAB or similar software to compute the FFT and inverse FFT invovled.

OK. So if we have a differential equation of the type:

$$au^{\prime \prime} + bu^{\prime} + cu = f(t)$$

The Equation referred to here (3.12) states that the discrete solution is given by:

$$\hat{u_j} = \frac{h^{2} \hat{f_j}}{aw^{j} + \beta + \gamma \bar{w}^{j}}$$

Here $h = 2 \pi/n$, $\beta = ch^2 + bh – 2a$, $\gamma = a-bh$, $w = \exp(2 \pi i/n)$, and $f_k = f(2 \pi k/n)$ for $1 \leq k \leq n-1$

I wrote the following algorithm in MATLAB to solve the given problem:

function u = diffu(n)
a = 1;
b = 2;
c = 2;
h = 2*pi/n;
beta = c*h^2 + b*h - 2*a;
gamma = a - b*h;
w = exp(2*pi*i/n);
cw = conj(w);
f = @(t) 2*cos(6*t);
k = 1:(n-1);
functionvalues = f(2*pi*k/n);
ftransform = fft(functionvalues);
utransform = zeros(1,(n-1));
for m = 1:(n-1)
   utransform(m) = ((h^2).*ftransform(m))./(a*(w.^m) + beta + gamma*(cw.^m));
end
u = ifft(utransform);

This algorithm works fine in that it does not give me any errors, and when I plot $u$ as a function of time, I get a nice trigonometric graph which is progressively smoother as I increase $n$. However, my question is as follows: where in this approximation are initial boundary values determined? After all, this type of second-order differential equation could have an infinite number of solutions based on how we vary initial conditions; yet when I plot $u$ as a function of time, a single graph is generated. So where exactly are initial conditions determined here? If anyone can explain this to me, I would be very grateful! I can also give some more information as to how Equation 3.12 is derived, if it is desired.

Best Answer

Better late than never, I guess.

The way you obtained Eq. 3.12 is by

1) discretizing the derivatives using Finite Differences

2) proposing $u_j$ to be expanded as a Discrete Fourier Transform

3) taking the FT of your discretization scheme

4) inserting your DFT approach in your scheme

5) solving for $\hat{u}_j$

But, like in any finite difference discretization, there are extra equations for the boundaries. You have to figure those out. For instance, suppose you wish to solve for $u_1$. How does the FD equation in step 1) look for this particular time step? There will be a $u_0$ somewhere in your discretization, so you insert the boundary value there, and then take the FT of the equation, finally solving for $u_1$. This should be done for each equation where a boundary value appears.

Hope this helps! Good luck!

Related Question