# How to solve PDEs for function $\Psi(z,\bar{z})$ dependent on independent variables $z,\bar{z}$? (Wirtinger derivatives)

For the past few days, I have been struggling to convey to mathematica to solve a PDE that is in terms of the independent variables $$(z,\bar{z})$$. I know mathematica supports solving PDEs with respect to $$(x,y)$$, but that conversion would be very very tedious.

The PDE I am looking at is an eigenvalue equation. I'm looking for the plane wave solutions of $$(-\Delta^{(H)}+\epsilon \hat{Q})\psi(z,\bar{z})= E \psi(z,\bar{z})$$ where $$\Delta^{(H)}= -(1-|z|^2)^2 \cdot \underbrace{4\frac{\partial}{\partial z}\frac{\partial}{\partial \bar{z}}}_{=\Delta}$$ is called the hyperbolic Laplacian, $$\Delta = \nabla^2$$ refers to the standard Laplacian in $$\mathbb{R}^2$$ and $$\epsilon>0$$ is a small parameter.

Denote by $$\frac{\partial}{\partial z}:= \partial_{z}\ ;\ \frac{\partial}{\partial \bar{z}}:= \bar{\partial_{z}}$$

The Operator $$\hat{Q}$$ is as follows $$\hat{Q}= \partial_{z}^2 \left((1-|z|^2)^3\partial_{z}\right)+\bar{\partial_{z}}^2 \left((1-|z|^2)^3\bar{\partial_{z}}\right)$$

That is the full problem I'm trying to solve but I wanted to start with a test case first. I am currently struggling to have mathematica consider $$z,\bar{z}$$ independently yet retain the fact that $$z,\bar{z}\in\mathbb{C}$$ and then solve the following eigenvalue problem.

Consider the simpler eigenvalue equation for which I know the solutions

$$-\Delta^{(H)} \psi=\varepsilon_{K}\psi$$ With $$K=ke^{i\beta}\in\mathbb{C}$$, The eigenfunctions have been analytically found to be $$\psi(z)=\left(\frac{1-|z|^2}{|1-ze^{-i\beta}|^2}\right)^{\frac{1+ik}{2}}$$ that satisfy the eigenvalue equation $$-\Delta^{(H)} \psi=\varepsilon_{K}\psi$$ with eigenvalue $$\varepsilon_{K}=-(1+k^2)$$.

I have verified that this is indeed the eigenfunctions however, It seems I am having trouble telling mathematica to solve this equation as a test case.

I have tried

eqn = -(1 - Abs[z]^2)^2*4*
D[D[u[z, Conjugate[z]], Conjugate[z]], z];
{vals, funs} =
NDEigensystem[{eqn},
u[z, Conjugate[z]], z \[Element] Disk[], 10];


but it is of no avail. I wanted to ask if I could get any help on this matter. Any help would be greatly appreciated!

I will only discuss the test case. An external reference is this PDF by Michael Stone. Some formulas for the Poincare disk model are on this Wikipedia page.

Translation to $$x,y$$ coordinates. Here $$z = x+iy$$. Then

DeltaH[f_] := -(1-x^2-y^2)^2*Laplacian[f,{x,y}];


One can see that it is invariant under rotation about the origin, but there is actually a bigger symmetry group, coming from the symmetries of the Poincare disk model, see the PDF above.

Eigenfunctions that OP has given. In $$x,y$$ coordinates these are

psi[k_,beta_] := ((1-x^2-y^2)/(1+x^2+y^2-2*x*Cos[beta]-2*y*Sin[beta]))^((1+I*k)/2);
lambda[k_,beta_] := -(1+k^2);


Let us check that these are eigenfunctions

- DeltaH[psi[k,beta]] - lambda[k,beta]*psi[k,beta] // Simplify
(* gives 0 *)


The operator is formally selfadjoint on some Hilbert space using an $$L^2$$-type inner product with measure $$dx\,dy / (1-x^2-y^2)^2$$. The functions psi[k,beta] are not in that Hilbert space, they have infinite norm. Therefore "eigenfunctions" is to be understood in a loose sense here.

Numerical computation of the eigenvalues. I will use this code:

lambdaN[a_,h_]:=First[NDEigensystem[{
-(1-x^2-y^2)^a*Laplacian[u[x,y],{x,y}],
DirichletCondition[u[x,y]==0,True]},
u[x,y],{x,y}\[Element]Disk[],6,
Method->{"PDEDiscretization"->{"FiniteElement",
"MeshOptions"->{"MaxCellMeasure"->h}}}]]//Chop[#,10^(-2)]&;


Here h determines a mesh size, the smaller the better. For DeltaH we must set a=2, but we will also consider the standard Laplacian a=0 and the intermediate value a=1 for comparison.

Standard Laplacian a=0:

lambdaN[0,0.01]
lambdaN[0,0.001]
lambdaN[0,0.0001]
(*
{5.78323,14.6827,14.6827,26.3784,26.3791,30.4787}
{5.78319,14.682,14.682,26.3747,26.3747,30.4713}
{5.78319,14.682,14.682,26.3746,26.3746,30.4713}
*)


The result is quite stable under reduction of mesh size. Compare with known exact values:

{BesselJZero[0,1],BesselJZero[1,1],BesselJZero[2,1],BesselJZero[0,2]}^2//N
(* {5.78319,14.682,26.3746,30.4713} *)


Intermediate case a=1:

lambdaN[1,0.01]
lambdaN[1,0.001]
lambdaN[1,0.0001]
(*
{4.00024,8.00574,8.00574,12.0312,12.0397,16.0613}
{4.00001,8.00017,8.00018,12.001,12.0012,16.002}
{4.,8.,8.,12.,12.,16.0001}
*)


We can again compare with exact results: The eigenfunctions include $$(1-r^2) r^{|k|} e^{ik \varphi}$$ with integer $$k$$, whose eigenvalue is $$4(1+|k|)$$. We have used polar coordinates $$x+iy = re^{i\varphi}$$.

The case of interest a=2:

lambdaN[2,0.01]
lambdaN[2,0.001]
lambdaN[2,0.0001]
(*
{2.50649,2.50649,2.97186,3.11274,3.53068,3.53068}
{1.5831,1.95856,1.96664,2.19987,2.22956,2.42164}
{1.66045,1.66154,1.79644,1.79762,1.90218,1.90218}
*)


We get some numbers, but they are not stable under changing the mesh size. This is expected: This operator has continuous spectrum. One would not expect NDEigensystem to give a meaningful result here. Perhaps someone familiar with NDEigensystem can explain how these numbers are computed.

Formulas. Since OP has mentioned this, here are formulas for the derivatives with respect to $$z$$ in terms of the partial derivatives with respect to $$x$$ and $$y$$:

Dz[f_]:=1/2*(D[f,x]-I*D[f,y]);
DzConjugate[f_]:=1/2*(D[f,x]+I*D[f,y]);

• Thank so much for your input! This is insightful, however does not answer the main question of solving the PDE with $(z,\bar{z})$ as variables since this is just switching into $(x,y)$ and solving it. Commented Aug 5, 2022 at 18:32