I need to solve numerically a PDE of the form $$ u_t(x,t)=u_{xx}(x,t)+u_x(x,t)^2-a(x)u_x(x,t)-a_x(x) $$ with initial condition $u(x,0)=u_0(x)$. I can assume that both $u(\cdot,t)$ and $a(\cdot)$ are periodic. I did an implementation of the method of lines and so far it works OK, but I was wondering whether I could improve on it.

So I converted the PDE into a coupled system of ODEs by expanding $u(x,t)$ into a Fourier series to obtain the corresponding system on the coefficients. This is, by expressing $u(x,t)=\sum_{n\in\mathbb Z} c_n(t) e^{inx}$ and $a(x)=\sum_{n\in\mathbb Z} a_n e^{inx}$, I have: \begin{eqnarray} u_t(x,t)&=&\sum_{n\in\mathbb Z} c_n'(t) e^{inx},\\ u_x(x,t)&=&\sum_{n\in\mathbb Z} c_n(t)in e^{inx},\\ u_{xx}(x,t)&=&-\sum_{n\in\mathbb Z} c_n(t)n^2 e^{inx},\\ a_x(x)&=&\sum_{n\in\mathbb Z} a_nin e^{inx},\\ a(x)u_{x}(x,t)&=&\sum_{n\in\mathbb Z} (a*u_x)_n e^{inx}=\sum_{n\in\mathbb Z} \big(\sum_{k\in\mathbb Z} a_k c_{n-k}(t)i(n-k)\big) e^{inx},\\ u_{x}(x,t)^2&=&\sum_{n\in\mathbb Z} (u_x*u_x)_n e^{inx}=-\sum_{n\in\mathbb Z} \big(\sum_{k\in\mathbb Z} c_k(t) c_{n-k}(t)k(n-k)\big) e^{inx}, \end{eqnarray} where $(a*b)_n=\sum_{k\in\mathbb Z} a_k b_{n-k}$. Then by matching the Fourier coefficients, the original PDE turns into the system of ODEs $$ c_n'(t)=-c_n(t)n^2-\sum_{k\in\mathbb Z} c_k(t) c_{n-k}(t)k(n-k)-\sum_{k\in\mathbb Z} a_k c_{n-k}(t)i(n-k)-a_nin,\quad n\in\mathbb Z. $$ What I do next is to solve the truncated system of ODEs for $n\in\{-N,\ldots,N\}$ with an ODE solver. This implies truncating the convolutions up to the coefficients with indexes in $\{-N,\ldots,N\}$. The solver receives the real coefficients $(a_0,a_1,\ldots,a_n,b_1,\ldots,b_n)$ corresponding to $c_n$ using equation (29) from here.

I am not getting the results that I might expect, and in particular when I increase $N$ the solution seems to shrink towards the initial condition. Here are some pictures illustrating this phenomena for simple choices of $a(\cdot)$ and $u_0(\cdot)$. (The method of lines gives the right solution.) Method of lines Fourier1 Fourier2 Fourier10

I wonder whether what I am experiencing is just a bug in my code or something is inherently wrong in the math. So my question is: is this a valid numerically scheme? I.e., is it guaranteed that when $N\to\infty$ the solution will converge to the solution of the PDE?

Hints on what can be causing the shrinking or related alternative approaches are very welcomed.

  • $\begingroup$ Very interesting question, would give +$\sin(\pi/2)$ if I could :) $\endgroup$ – mathreadler Oct 5 '15 at 16:45
  • 1
    $\begingroup$ This technique is usually referred to as a spectral method, and it is really valid. I suspect a bug in your convolution code. Try testing it by solving some simpler equations like $u_t = u_{xx}$, $u_t = a(x)u_x$, etc. $\endgroup$ – uranix Oct 6 '15 at 21:14
  • $\begingroup$ Thanks for pointing this out @uranix. $\endgroup$ – epsilone Oct 6 '15 at 21:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.