Hello,

I have just started learning Julia, and I am attempting to use it to solve a one-dimensional eigenvalue problem in my field, plasma physics (specifically, the tearing mode, if anyone is interested). The problem involves a very thin boundary layer around the center of the domain, which might be the source of the difficulties I’m having.

Since I’m only interested in the growing mode, i.e. the eigenvalue of the problem with the largest real part, after looking around I thought to use KrylovKit to do this. However, it seems to have quite a bit of trouble converging, even when my resolution is sufficient to get a good solution using brute force (just using Julia’s eigen()). Here is my code:

```
f(x)=tanh(x)
ddf(x)=-2*tanh(x)*(sech(x))^2
function vargrid(dxmin,xlim,eps)
if eps > 0.0
n = convert(Int64,floor(log(xlim*eps/dxmin+1)/log((1+eps))))+1
else
n = convert(Int64,xlim/dxmin)+1
end
dxmax = dxmin*(1+eps)^n
x=zeros(2*n-1)
dx=zeros(2*n)
dx[n+1]=dxmin
dx[n]=dxmin
for i in (n+1):2*n-1
dx[i+1]=dx[i]*(1+eps)
x[i]=x[i-1]+dx[i]
end
for i in n-1:-1:1
dx[i]=dx[2*n+1-i]
x[i]=-x[2*n-i]
end
return x,dx
end
function tearingVG(S,K,ar,dxmin,xlim,eps)
x,dx=vargrid(dxmin,xlim,eps)
N=length(x)
DD=zeros(N,N)
for i in 1:N, j in 1:N
i-j==0 && ( DD[i,j] = -2.0*(1.0+dx[i+1]/dx[i])/(dx[i+1]^2+dx[i]*dx[i+1]) )
(i-j)==-1 && ( DD[i,j] = 2.0/(dx[i+1]^2+dx[i]*dx[i+1]) )
(i-j)==1 && ( DD[i,j] = 2.0*(dx[i+1]/dx[i])/(dx[i+1]^2+dx[i]*dx[i+1]) )
end
DD[1,1]=DD[1,1]+2.0*((dx[2]/dx[1])/(1+K*dx[1]))/(dx[2]^2+dx[1]*dx[2])
DD[N,N]=DD[N,N]+2.0/(1+K*dx[N+1])/(dx[N+1]^2+dx[N]*dx[N+1])
L=DD-I*K^2
A=zeros(2*N,2*N)
B=zeros(2*N,2*N)
A=complex(A)
B=complex(B)
B[1:N,1:N]=1.0*I(N)
B[N+1:end,N+1:end]=L
A[1:N,1:N]=L/S - im*ar*diagm(f.(x))*K
A[1:N,N+1:end]=im*diagm(f.(x))*K
A[N+1:end,1:N]=im*diagm(f.(x))*K*L-im*diagm(ddf.(x))*K
A[N+1:end,N+1:end]=-im*ar*diagm(f.(x))*K*L+im*ar*diagm(ddf.(x))*K
return eigen(A,B)
end
function tearing_KrylovKit(S,K,ar,dxmin,xlim,eps)
x,dx=vargrid(dxmin,xlim,eps)
N=length(x)
DD=zeros(N,N)
for i in 1:N, j in 1:N
i-j==0 && ( DD[i,j] = -2.0*(1.0+dx[i+1]/dx[i])/(dx[i+1]^2+dx[i]*dx[i+1]) )
(i-j)==-1 && ( DD[i,j] = 2.0/(dx[i+1]^2+dx[i]*dx[i+1]) )
(i-j)==1 && ( DD[i,j] = 2.0*(dx[i+1]/dx[i])/(dx[i+1]^2+dx[i]*dx[i+1]) )
end
DD[1,1]=DD[1,1]+2.0*((dx[2]/dx[1])/(1+K*dx[1]))/(dx[2]^2+dx[1]*dx[2])
DD[N,N]=DD[N,N]+2.0/(1+K*dx[N+1])/(dx[N+1]^2+dx[N]*dx[N+1])
L=DD-I*K^2
A=zeros(2*N,2*N)
B=zeros(2*N,2*N)
A=complex(A)
B=complex(B)
B[1:N,1:N]=1.0*I(N)
B[N+1:end,N+1:end]=L
A[1:N,1:N]=L/S - im*ar*diagm(f.(x))*K
A[1:N,N+1:end]=im*diagm(f.(x))*K
A[N+1:end,1:N]=im*diagm(f.(x))*K*L-im*diagm(ddf.(x))*K
A[N+1:end,N+1:end]=-im*ar*diagm(f.(x))*K*L+im*ar*diagm(ddf.(x))*K
C=inv(B)*A
return eigsolve(C,1,:LR,maxiter=1000)
end
```

`vargrid`

sets up the grid, while the other two functions use either `eigen`

or KrylovKit’s `eigsolve`

to solve the problem. Typical values are

S ~ 10^4-10^8

K ~ 1-10^-4

ar=0 (for now)

dxmin=10^-3 - 10^-5

xlim=10

eps=0.02

The typical values of the eigenvalue I’m looking for are small but positive and real, while typically all other eigenvalues have negative real parts. I am not experienced with Krylov methods, so any advice would be greatly appreciated; I thought that this would work quite well since the required eigenvalue is on the edge of the spectrum, but I’m sure there is some subtlety I’m overlooking. Thanks in advance!