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!
Best Answer
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.