Solving a PDE with tricky boundary conditions

Question

I'm reading this paper by Alexander Bloemendal currently and focusing on solving the following PDE numerically: \begin{align*} \frac {\partial F}{\partial x}+\frac {2}{\beta}\frac {\partial^{2} F}{\partial \omega^{2}}+\left(x-\omega^2\right)\frac {\partial F}{\partial \omega}=0\quad \text{for }(x,\omega)\in \mathbb {R}^{2}, \end{align*} \begin{align*} &F(x,\omega)\to 1\quad \text{as }x,\omega\to \infty\text{ together,}\\ &F(x,\omega)\to 0\quad \text{as }\omega\to -\infty\text{ with }x\text{ bounded above.} \end{align*} This is a fairly standard diffusion-advection equation with space variable $\omega$ and time variable $-x$. The main apparent difficulty with the formulation of the boundary value problem is that the boundary conditions and the desired slice of the solution are all at infinity. Therefore, a change of variable is performed, $\omega=-\cot \theta$, and we get \begin{align*} \frac {\partial F}{\partial x}+\left(\frac {2}{\beta}\sin^{4} \theta\right)\frac {\partial^{2} F}{\partial \theta^{2}}+\left(\left(x+\frac {2}{\beta}\sin 2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right)\frac {\partial F}{\partial \theta}=0, \end{align*} \begin{align*} &F(x,\theta)\to 1\quad \text{as }x\to \infty\text{ with }\theta\geq\theta_{0}>0\\ &F(x,\theta)= 0\quad \text{on }\theta=0. \end{align*} The above system can be solve using the $\texttt{NDSolve}$ package in Mathematica. The details are in Chapter 6 of this paper. Now, I want to solve it using finite difference method instead. However, there are some issues I encountered:

$1$. The first issue is on forming the iteration matrix. I can replace $\partial^{2}F/\partial \theta^{2}$ with the centered difference \begin{align*} \frac {F(x,\theta+h)-2F(x,\theta)+F(x,\theta-h)}{h^2} \end{align*} and $\partial F/\partial \theta$ with either forward difference \begin{align*} \frac {F(x,\theta+h)-F(x,\theta)}{h} \end{align*} or backward difference \begin{align*} \frac {F(x,\theta)-F(x,\theta-h)}{h}. \end{align*} However, the resulting iteration matrix will involve both $\theta$ and $x$. Say if I set the initial condition to be: \begin{align*} F(x_{0},\theta)=\begin{cases} \Phi\left(\frac {x_{0}-\cot^{2}\theta}{\sqrt{\left(4/\beta\right)\cot \theta}}\right)\quad &0\leq \theta\leq \pi/2\\ 1\quad &\theta\geq \pi/2 \end{cases} \end{align*} at the initial time $x_{0}=10$. I will perform the iteration until the final time $x_{1}=-10$, and I will focus on the interval between $\theta=0$ and $\theta_{1}=2\pi$. I also set $\Delta x=0.01$ and $\Delta \theta=0.1$. If I use forward Euler to step forward, since the iteration matrix will involve both $\theta$ and $x$, I don't know what values $\theta$ and $x$ should take in each iteration.

$2$. My second concern is on the boundary conditions. The first one is trivial, which is just $F(x,\omega)=0$ at $\theta=0$. The second one is a little bit tricky since according to what Alexander said in his paper, if we are interested in values $\theta\leq k\pi$, it seems wise to use $\theta_{1}=(k+1)\pi$ at least. Therefore, if I want to focus on the interval $[0,2\pi]$, then I need to set the other boundary condition at $\theta_{1}=3\pi$, which means that I'm actually focusing on the interval $[0,3\pi]$ instead. But I guess this should be fine?

I will keep working on this problem throughout this summer and any advises will be appreciated!

Answer 1

I am not sure if their paper is right, after I do the change of variable, I think the equation should be \begin{align*} \frac {\partial F}{\partial x}+\left(\frac {2}{\beta}\sin^{4} \theta\right)\frac {\partial^{2} F}{\partial \theta^{2}}+\left(\left(x+\frac {2}{\beta}\sin 2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right)\frac {\partial F}{\partial \theta}=0,\quad 0<\theta<\pi,\quad x<10 \end{align*} \begin{align*} F(10,\theta)=\begin{cases} \Phi\left(\frac {10-\cot^{2}\theta}{\sqrt{\left(4/\beta\right)\cot \theta}}\right)\quad &0\leq \theta\leq \pi/2\\ 1\quad &\pi/2\leq \theta\leq \pi \end{cases} \end{align*} \begin{equation*} F(x,0)=0,\quad x<10 \end{equation*} where we use $-\cot(\pi)=\infty$.

However, they might do the extension to $F(x,\theta)$ and force the $F(x,\theta)$ is a constant in $\theta$ when $\theta>\pi$. Anyway, if you want to numerically solve the following equation: \begin{align*} \frac {\partial F}{\partial x}+\left(\frac {2}{\beta}\sin^{4} \theta\right)\frac {\partial^{2} F}{\partial \theta^{2}}+\left(\left(x+\frac {2}{\beta}\sin 2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right)\frac {\partial F}{\partial \theta}=0,\quad 0<\theta<\infty,\quad x<10 \end{align*} \begin{align*} F(10,\theta)=\begin{cases} \Phi\left(\frac {10-\cot^{2}\theta}{\sqrt{\left(4/\beta\right)\cot \theta}}\right)\quad &0\leq \theta\leq \pi/2\\ 1\quad &\pi/2\leq \theta \end{cases} \end{align*} \begin{equation*} F(x,0)=0,\quad x<10 \end{equation*}

You can use the finite-difference method.

To use the finite-difference method, we first need to set up grid points: For example, you can define $$ x^0=10,\quad x^n=x^0-n\Delta x=10-n\Delta x $$ and $$ \theta_0=0,\quad \theta_m=\theta_m+m\Delta\theta\,. $$ where $m,n\in\mathbb{N}$. Then we denote the approximation of $F(x^n,\theta_m)$ by $F^n_m$.

Since you want solution at $x=-10$ with $\theta\in[0,2\pi]$, we denote denote $M=\frac{2\pi}{\Delta \theta}$ the number of grid points and $N=\frac{20}{\Delta x}$ the steps you need to run. Because it only has one side boundary condition $\theta=0$, we should use forward difference. More specifically, from step $n$ to $n+1$, according to the equation, we can get $F^{n+1}_m$ by solving \begin{equation} \begin{aligned} &\frac{F^{n+1}_m-F^{n}_m}{\Delta x}=\left(\frac {2}{\beta}\sin^{4} \theta_m\right)\frac {F^n_{m+1}-2F^n_m+F^n_{m-1}}{\Delta \theta}\\ &+\left(\left(x^n+\frac {2}{\beta}\sin 2\theta_m\right)\sin^{2}\theta_m-\cos^{2}\theta_m\right)\frac {F^n_{m+1}-F^n_m}{\Delta \theta}=0 \end{aligned}\,,\quad 0<m<M,\ 0\leq n\leq N-1\tag{1} \end{equation} And we always let $F^{n+1}_0=0$. However, there is a problem. From (1), you can notice that we haven't updated $F^{n+1}_M$. For example, you don't know how to calculate $F^1_{M}$ since you need to know $F^0_{M+1}$ if you want to use (1).

To overcome this problem, we need to use the initial condition is defined for all $\theta\geq 0$. According to the initial condition, we know the value of $$ F^0_m,\quad \forall 0\leq m\leq N+M\,. $$ Then, in each step, instead of (1), we can use \begin{equation} \begin{aligned} &\frac{F^{n+1}_m-F^{n}_m}{\Delta x}=\left(\frac {2}{\beta}\sin^{4} \theta_m\right)\frac {F^n_{m+1}-2F^n_m+F^n_{m-1}}{\Delta \theta}\\ &+\left(\left(x^n+\frac {2}{\beta}\sin 2\theta_m\right)\sin^{2}\theta_m-\cos^{2}\theta_m\right)\frac {F^n_{m+1}-F^n_m}{\Delta \theta}=0 \end{aligned}\,,\quad 0<m<M+N-n,0\leq n\leq N-1\tag{2} \end{equation} where we notice that $m$'s upper bound is changed to $M+N-n$. The reason is similar to before. Since we know $F^0_m$ for $0\leq m\leq M+N$. Using (1), we can obtain $F^1_m$ for $0\leq m\leq M+N-1$. Iteratively, we can obtain $F^n_m$ for $0\leq m\leq M+N-n$. Using (2), we can always update $F^{n}_M$ when $n\leq N$. Thus, we can obtain an approximation to $F(-10,\theta)$ for $\theta\in[0,2\pi]$. I think their paper also means similar idea by focusing on a larger interval.