Heat equation with variable coefficients

[Telemachus]

$\displaystyle \frac{\partial u(x,t)}{\partial t}-p(x) \frac{\partial^2 u(x)}{\partial x^2}-\frac{\partial p(x)}{\partial x}\frac{\partial u(x)}{\partial x}=f(x,t)$.

With periodic boundary conditions: $u(x,t)=u(x+2\pi,t)$. The right hand side is also a periodic function of space, and it can be saparated: $f(x,t)=f(x)T(t)$. I wanted to do this by means of Fourier series, so I've expanded everything in Fourier series:

$\displaystyle u(x)=\sum_{k=-\infty}^{\infty} \hat u_k e^{ikx},\frac{\partial u(x)}{\partial x}=\sum_{k=-\infty}^{\infty} ik \hat u_k e^{ikx},\frac{\partial^2 u(x)}{\partial x^2}=-\sum_{k=-\infty}^{\infty} k^2 \hat u_k e^{ikx},f(x)=\sum_{k=-\infty}^{\infty} \hat f_k e^{ikx},p(x)=\sum_{k=-\infty}^{\infty} \hat p_k e^{ikx},\frac{\partial p(x)}{\partial x}=\sum_{k=-\infty}^{\infty} ik \hat p_k e^{ikx}$.

So when I replace all this in the differential equation, with the aim to get a differential equation for the $\hat u_k$ in Fourier space, I get some products of the form:

$\displaystyle \left( \sum_{k=-\infty}^{\infty} \hat p_k e^{ikx} \right ) \left( -\sum_{k=-\infty}^{\infty} k^2 \hat u_k e^{ikx} \right)$ and
$\displaystyle \left( \sum_{k=-\infty}^{\infty} ik \hat p_k e^{ikx} \right) \left( \sum_{k=-\infty}^{\infty} ik \hat u_k e^{ikx} \right)$,

So, my idea was (which worked for the constant coefficients situation) then project with $\int_0^{2\pi}e^{-ik'x}dx$ in the equation, and use the orthogonality relation $\int_0^{2\pi} e^{ikx}e^{-ik'x}dx=2\pi \delta_{k,k'}$ in order to get an equation for the $\hat u_k$. However, in this situation it is not clear to me what should I get from:
$\displaystyle \int_0^{2\pi}e^{-ik'x}dx \left( \sum_{k=-\infty}^{\infty} \hat p_k e^{ikx} \right ) \left( -\sum_{k=-\infty}^{\infty} k^2 \hat u_k e^{ikx} \right)$

because of the product of sums. So, I'm not sure if this is a trivial question and I am not seeing it, or if I should look for another approach.

Any comments and ideas are welcome, thanks in advance.

 

[Orodruin]

What is your function $p(x)$?

Due to the appearance of $p(x)$, you will obtain coupled ODEs rather than decoupled ones, i.e., you are working in a basis where the basis functions are not eigenfunctions of your differential operator. The spatial operator in your problem seems to be a periodic Sturm-Liouville operator $-\partial_x p(x) \partial_x$. If you want to get rid off the cross terms, you should expand the solution in terms of the eigenfunctions of that operator, not in terms of $e^{ikx}$, which are the eigenfunctions of $-\partial_x^2$ - a different Sturm-Liouville operator.

 

[Telemachus]

I haven't defined p(x), but the idea was to use a periodic function in order to have everything smooth and periodic in the equation. However, it seems like I have chosen the wrong strategy. I mean, I want to use Fourier series (not another function space, because I am thinking in some general p(x), and I don't want to generate this space at every time I change p(x)).

 

[Orodruin]

It is still the same function space ($L^2$ integrable periodic functions with weight function one). The difference is the basis you use. If you are not willing to change the basis when you change the operator, then I am afraid you will be stuck with the coupled differential equations.

Wanting to use the same basis functions for different differential operators to decouple the differential equations is akin to wanting all $n\times n$ matrices to have the same eigenvectors. If you take a new matrix, you will have to compute its eigenvectors and not use the eigenvectors of a different matrix. If you take a new differential operator, you will have to use its eigenfunctions to decouple the differential equation, no use the eigenfunctions of a different operator.

Edit: Note that your sums contain expressions on the form
$$
\int_0^{2\pi}e^{-ik'x}dx \left( \sum_{k=-\infty}^{\infty} \hat p_k e^{ikx} \right ) \left( -\sum_{k=-\infty}^{\infty} \ell^2 \hat u_k e^{ikx} \right)
=
\int_0^{2\pi} e^{-ik'x}dx \sum_{k=-\infty}^\infty \sum_{\ell = -\infty}^\infty \hat p_k \ell^2 \hat u_\ell e^{i(k+\ell)x}
= \sum_{k=-\infty}^\infty \sum_{\ell = -\infty}^\infty \hat p_k (k-k')^2 \hat u_\ell 2\pi \delta_{k',k+\ell}
= \sum_{k=-\infty}^\infty \hat p_k (k-k')^2 \hat u_{k'-k}.
$$
Thus, for all $k$, the differential equation is a priori coupled with all other differential equations (assuming that the corresponding coefficients are non-zero - for $p = 1$ only the $k = 0$ term is non-zero and you do decouple the differential equation).

 

[Telemachus]

But there is no other approach I can use, like a collocation method? I mean, for example, I could start by discretization of the time operator, just using some finite difference scheme in time, and arrive to a linear system of equations to solve it numerically.

This is what I mean, I could use a Crank-Nicolson scheme:
$\displaystyle \frac{u^{n+1}-u^{n}}{\Delta t}-p(x) \frac{\partial^2 }{\partial x^2} \left(\frac{u^{n+1}-u^{n}}{2} \right) -\frac{\partial p(x)}{\partial x} \frac{\partial }{\partial x} \left( \frac{u^{n+1}-u^{n}}{2} \right) =\frac{f^{n+1}+f^{n}}{2}$,

So:
$\displaystyle 2u^{n+1}-\Delta t p(x) \frac{\partial^2 u^{n+1}}{\partial x^2} - \Delta t \frac{\partial p(x)}{\partial x} \frac{\partial u^{n+1}}{\partial x} = \Delta t (f^{n+1}+f^{n}) + 2u^n - \Delta t p(x) \frac{\partial^2 }{\partial x^2} u^n -\Delta t \frac{\partial p(x)}{\partial x}u^n$

And then purpose a solution of the type: $u=\sum \hat u_k^n e^{ikx}$ at a set of grid points $x_j=(j-1)h,\, j=1,2,...,N$

BTW, I've seen this approach in this paper: https://arxiv.org/abs/1209.0751 (eq. 22).

However, I'm not sure on how to implement it. I'll keep studying and I'll see if I can construct the linear system on my own. I thought it was a good idea to talk about it in this forum, because every time I come I always get the ideas more clear.

Alternatively, is it possible to solve the coupled system with an iterative scheme?

 

[Chestermiller]

Shouldn't those be + signs between the u's in the spatial terms of the first equation?

 

[Orodruin]

Shouldn't those be + signs between the u's in the spatial terms of the first equation?

Why? Without the minus signs, the spatial part is not a Sturm-Liouville operator (it is the negative of a SL operator) and would have negative eigenvalues. If there were plusses instead it would no longer be a generalisation of the heat equation and the solutions (with a constant source) would grow exponentially rather than settle around the stationary solution.

Let $p(x) = a$ and you recover the heat equation
$$
u_t - a u_{xx} = f(x,t)
$$
with $f(x,t)$ as the source term.

 

[Chestermiller]

Why? Without the minus signs, the spatial part is not a Sturm-Liouville operator (it is the negative of a SL operator) and would have negative eigenvalues. If there were plusses instead it would no longer be a generalisation of the heat equation and the solutions (with a constant source) would grow exponentially rather than settle around the stationary solution.

Let $p(x) = a$ and you recover the heat equation
$$
u_t - a u_{xx} = f(x,t)
$$
with $f(x,t)$ as the source term.

The Crank-Nicholson scheme uses values (arithmetically) averaged over the time interval to express the transport terms. This makes it second order accuracy with respect to time.

 

[Orodruin]

What does the Crank-Nicolson method have to do with whether or not the given equation is the proper generalisation of the heat equation or not?

Without the minuses, the PDE in the OP is simply not a generalisation of the heat equation.

 

[Chestermiller]

What does the Crank-Nicolson method have to do with whether or not the given equation is the proper generalisation of the heat equation or not?

Without the minuses, the PDE in the OP is simply not a generalisation of the heat equation.

We're talking about the following equation, correct?
$$\displaystyle \frac{u^{n+1}-u^{n}}{\Delta t}-p(x) \frac{\partial^2 }{\partial x^2} \left(\frac{u^{n+1}-u^{n}}{2} \right) -\frac{\partial p(x)}{\partial x} \frac{\partial }{\partial x} \left( \frac{u^{n+1}-u^{n}}{2} \right) =\frac{f^{n+1}+f^{n}}{2}$$

 

[Orodruin]

No, I thought you were referring to the original PDE in the OP. For some reason I missed the alert from post #5 and did not know it existed. The only sense I could make of it was in reference to post #1, which seemed strange and is what I reacted to. I agree that the equation you quote in #10 should have the same sign for $u^{n+1}$ and $u^n$ (in the spatial part).

 

[Orodruin]

But there is no other approach I can use, like a collocation method? I mean, for example, I could start by discretization of the time operator, just using some finite difference scheme in time, and arrive to a linear system of equations to solve it numerically.

This is what I mean, I could use a Crank-Nicolson scheme:
$\displaystyle \frac{u^{n+1}-u^{n}}{\Delta t}-p(x) \frac{\partial^2 }{\partial x^2} \left(\frac{u^{n+1}-u^{n}}{2} \right) -\frac{\partial p(x)}{\partial x} \frac{\partial }{\partial x} \left( \frac{u^{n+1}-u^{n}}{2} \right) =\frac{f^{n+1}+f^{n}}{2}$,

So:
$\displaystyle 2u^{n+1}-\Delta t p(x) \frac{\partial^2 u^{n+1}}{\partial x^2} - \Delta t \frac{\partial p(x)}{\partial x} \frac{\partial u^{n+1}}{\partial x} = \Delta t (f^{n+1}+f^{n}) + 2u^n - \Delta t p(x) \frac{\partial^2 }{\partial x^2} u^n -\Delta t \frac{\partial p(x)}{\partial x}u^n$

And then purpose a solution of the type: $u=\sum \hat u_k^n e^{ikx}$ at a set of grid points $x_j=(j-1)h,\, j=1,2,...,N$

BTW, I've seen this approach in this paper: https://arxiv.org/abs/1209.0751 (eq. 22).

However, I'm not sure on how to implement it. I'll keep studying and I'll see if I can construct the linear system on my own. I thought it was a good idea to talk about it in this forum, because every time I come I always get the ideas more clear.

Alternatively, is it possible to solve the coupled system with an iterative scheme?

If you anyway want a numerical method, I would suggest discretising your PDE in $x$ rather than $t$ with the spatial derivatives replaced by finite differences. The result would be a finite set of linear partial differential equations in time where the SL operator is replaced by a matrix acting on a vector containing the values of $u$ at the discretisation points. This differential equation can be solved exactly by diagonalising a finite $n \times n$ matrix. It essentially corresponds to truncating your Fourier series at a finite $k$.

 

[Telemachus]

Shouldn't those be + signs between the u's in the spatial terms of the first equation?

Yes, I'm sorry. It was too late in my country and I was tired when I wrote that. You are right. It should read:

$$\displaystyle \frac{u^{n+1}-u^{n}}{\Delta t}-p(x) \frac{\partial^2 }{\partial x^2} \left(\frac{u^{n+1}+u^{n}}{2} \right) -\frac{\partial p(x)}{\partial x} \frac{\partial }{\partial x} \left( \frac{u^{n+1}+u^{n}}{2} \right) =\frac{f^{n+1}+f^{n}}{2}$$

 

[Telemachus]

If you anyway want a numerical method, I would suggest discretising your PDE in $x$ rather than $t$ with the spatial derivatives replaced by finite differences. The result would be a finite set of linear partial differential equations in time where the SL operator is replaced by a matrix acting on a vector containing the values of $u$ at the discretisation points. This differential equation can be solved exactly by diagonalising a finite $n \times n$ matrix. It essentially corresponds to truncating your Fourier series at a finite $k$.

Yes, but I wanted to do exactly what is done in that paper, because it is computationally better when you go to higher dimensions (with the ADI scheme you gain a lot in computational terms, and with that method it gives good accuracy and stability). My idea was to solve the problem I posed at the beginning of the thread to have the tools to apply it in higher dimensions as it is done in that paper (diagonalization of a matrix is too expensive in higher dimensions, I thought of that too). But what is done in that paper seems quiet complicated, so I don't know if I'll go further. I don't exactly understand what it is done in equation (22), when they talk about the Fourier collocation. The equation is there, the system is there, but I don't know what the matrix should look like.

 

[hilbert2]

For that kind of an equation, the static solutions for which $\frac{\partial u}{\partial t} = 0$ at all points $(x,t)$ are not the usual linearly changing solutions $u(x,t) = A + Bx$ of the constant coefficient heat equation. The static solutions are easier to find (by setting the time derivative to value zero) than the solutions of the whole PDE and are a good test for any numerical scheme.

 

[Telemachus]

Ok, I think that I know what I have to do now. Thank you all.