# Solving a PDE with tricky boundary conditions

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!

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{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 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{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 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.

• Hi, thank you for your answer. I have two questions. The first one is I wonder how did you get $M+N$ to be the upper bound of $m$ and subtract $n$ from it in each iteration? The second one is I wonder if $x$ should be $x^n$ in each iteration. Commented Jul 25, 2021 at 2:40
• @Jiexiong687691 For the first question, the upper bound for $m$ in $n$-th step is $(M+N-n)$, the reason is stated in (1). For the second question, yes, I fixed the typo. Commented Jul 25, 2021 at 14:34
• @Jiexiong687691 I guess I didn't answer your first question clear. At the initial condition, we can choose any upper bound for $m$ at $n=0$ since we know F(10,θ) for every θ. Then in (1), we can see that if you obtain $F^{n+1}_m$, we need to know $F^n_{m+1}$. That's why we need to subtract $n$ from the upper bound in $n$-th step. Hope this can help you. Commented Jul 25, 2021 at 19:35