Numerisk løsning av diff.ligning ved Fast Fourier Transform

[krje1980]

Hei.

Har et lite spørsmål relatert til følgende oppgave:

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. Så hvis vi har diff-ligninger av typen:

$a{u}^{\prime \prime }+b{u}^{\prime }+cu=f(t)$

så sier ligningen det refereres til her (3.12) at den diskré løsningen er gitt ved:

$\widehat{u_j}=\frac{h^2\widehat{f_j}}{a{w}^j+\beta +\gamma {\overline{w}}^j}$

Her er $h=2\pi /n$, $\beta =c{h}^2+bh-2a$, $\gamma =a-bh$, $w=\exp (2\pi i/n)$, og $f_k=f(2\pi k/n)$ for $1\le k\le n-1$

Jeg klarer å fint å skrive en MATLAB algoritme som løser dette for et gitt input av $n$. Først defineres alle variablene, før jeg i en vektor regner ut alle verdiene for $f$. Deretter bruker jeg MATLAB sin fft-rutine til å finne den diskré Fourier transformen til $f$. Og til slutt setter jeg opp alt slik som i likning 3.12, og tar den inverse transformasjonen til uttrykket.

Det jeg imidlertid lurer på er følgende: Når jeg skriver dette som en algoritme og senere plotter løsningen $u$ som en funksjon av tid, så får jeg en fin trigonometrisk kurve. Men med slike andreordens diff.ligninger vil det jo kunne være et uendelig antall mulige løsningskurver basert på initialverdiene. Så hvor kommer initialverdiene inn i denne oppgaven? Som sagt så får jeg en fin løsningskurve når jeg plotter dette, men kan ikke se noen sted at vi tar hensyn til initialverdier. Hva er det da som gjør at vi får en unik løsning?

Dersom noen med peiling på dette kan svare meg så vil jeg være ekstremt takknemlig! Jeg kan godt supplementere med mer informasjon angående hvordan formelen over et utredet, samt MATLAB koden min, dersom det er ønskelig. Men hvis noen kan svare meg på dette basert på det jeg har skrevet over, så hadde det vært glimrende!

[Gustav]

kan du skissere hvordan du har kommet frem til formelen for u_j ?

[krje1980]

Hei.

Fra læreboken:

As anoter application of the DFT, we consider the differential equation $a{u}^{\prime \prime }+b{u}^{\prime }+cu=f(t)$, and take $f$ to be a continuous, $2\pi$-periodic function of $t$. There is a well-known analytic method for finding the unique periodic solution of this equation (cf. [3. §3.7.2]), provided that $f$ is known for all $t$. On the other hand, if we only know $f$ at the points $t_j=2\pi j/n$, for some integer $n\ge 1$, this method is no longer applicable.

Instead of trying to work directly with the differential equation itself, we work with a discretized version of it. There are many ways of discretizing. With $h=2\pi /n$, the one that we use here involves the following approximations:

$u^{\prime }(t)\approx \frac{u(t)-u(t-h)}{h}$

$u^{\prime \prime }(t)\approx \frac{u(t+h)+u(t-h)-2u(t)}{h^2}$

We let $t_k=2\pi k/n$ and $u_k=u(t_k)$ for $0\le k\le n$. Since $u$ is periodic, $u_n=u_0$.

At $t=t_k$, the preceding difference formulas become

$u^{\prime }(t_k)\approx \frac{u_k-u_{k-1}}{h}$

$u^{\prime \prime }(t_k)\approx \frac{u_{k+1}+u_{k-1}-2u_k}{h^2}$

for $1\le k\le n-1$. Substituting these values into the differential equation and collecting terms yields the following difference equation:

$a{u}_{k+1}+\beta u_k+\gamma u_{k-1}=h^2{f}_k$ for $1\le k\le n-1$

where $f_k=f(2\pi k/n),\beta =c{h}^2+bh-2a$, and $\gamma =a-bh$.

Suppose that $u\in S_n$ is a solution to the difference equation above. Let $[tex]$\hat{u} = \mathcal{F}[/tex] and $\hat{f}=\mathcal{F}[f]$ and $w=\exp (2\pi i/n)$. By taking the DFT of both since and using the first part of Theorem 3.4 we can show that

$\widehat{u_j}=\frac{h^2\widehat{f_j}}{a{w}^j+\beta +\gamma {\overline{w}}^j}$

provided that $a{w}^j+\beta +\gamma {\overline{w}}^j$ is never $0$. This gives us $\hat{u}$. We get $u$ by applying the inverse DFT to $\hat{u}$

[Gustav]

Ok, du har jo her periodiske grensebetingelser på u: $u_0=u_n$. Dette gir da en unik løsning på diff.ligningen slik jeg tolker det.

[krje1980]

Ah OK. Sevlsagt. Da er det mer logisk. Men hvor er det vi her definerer hva verdien skal være for $u_0$ og $u_n$?