Solving the Bloemendal-Virag equation with ultraspherical spectral method

I’m trying to solve the following PDE over x\in \left[-10,\frac {13}{\sqrt{\beta}}\right] and \theta\in \left[0,\pi\right] using the ultraspherical spectral method:

\begin{equation} \frac {\partial H}{\partial x}+\left(\frac {2}{\beta}\sin^{4} \theta\right)\frac {\partial^{2} H}{\partial \theta^{2}}+\left[\left(x+\frac {2}{\beta}\sin 2\theta\right)\sin^{2}\theta-\cos^{2}\theta\right]\frac {\partial H}{\partial \theta}=0, \end{equation}

with boundary condition

\begin{align*} H(x,0)=0, \end{align*}

and the approximate asymptotic initial condition

\begin{align*} H\left(\frac {13}{\sqrt{\beta}},\theta\right)=\begin{cases} \Phi\left(\frac {\frac {13}{\sqrt{\beta}}-\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*}

Here \Phi denotes the standard normal distribution function. Following the sample codes in this Notebook, my code for this problem is the following:

β=2;
dθ = 0 ..π; 
dx = -10 ..13/sqrt(β);
d = dθ×dx;
θ,x = Fun(d);
Dθ = Derivative(d,[1,0]);
Dx = Derivative(d,[0,1]);
Φ = xx -> erf.(xx/sqrt(2))/2 + 0.5;
g = (t) -> Φ( (13/sqrt(β) - cot(t).^2)./sqrt.(4/β * cot.(t)) );
h0 = (t) -> t < pi/2 ? g(t) : 1.0;
u0 = Fun(θ->h0(θ),dθ);
u =\([ldirichlet(dx)⊗I; ldirichlet(dθ)⊗I;
Dx+((2/β)*(sin(θ))^4)*(Dθ^2)+((x+(2/β)*sin(2θ))*sin(θ)^2-cos(θ)^2)*Dθ],
[u0; 0; 0];tolerance=1E-4);

When I run the code, I got the following error:

AssertionError: length(order) == 2

I wonder how to fix it. Thanks.

1 Like

Looks like there are some bugs here. Meanwhile, inserting the tensor products manually gets you some distance, but there’s still some error in the solve that needs to be looked into.

β=2;
dθ = 0 ..π;
dx = -10 ..13/sqrt(β);
d = dθ×dx;
θ = Fun(dθ);
x = Fun(dx);
Dθ = Derivative(space(θ)) ⊗ I;
Dx = I ⊗ Derivative(space(x));
Φ = xx -> erf.(xx/sqrt(2))/2 + 0.5;
g = (t) -> Φ( (13/sqrt(β) - cot(t).^2)./sqrt.(4/β * cot.(t)) );
h0 = (t) -> t < pi/2 ? g(t) : 1.0;
u0 = Fun(θ->h0(θ),dθ);
u =\([ldirichlet(dx)⊗I; ldirichlet(dθ)⊗I;
Dx + (Multiplication((2/β)*(sin(θ))^4) ⊗ I) * (Dθ^2) +
	(((I ⊗ Multiplication(x)) + (Multiplication((2/β)*sin(2θ)) ⊗ I)) * (Multiplication(sin(θ)^2) ⊗ I)
		- (Multiplication(cos(θ)^2) ⊗ I)) * Dθ],
			[u0; 0; 0];tolerance=1E-4);

Could you file an issue with ApproxFun.jl?

1 Like

Thanks! I’ll do that.