4
$\begingroup$

I am trying to solve numerically the GL equations in cylindrical coordinates. I already know what the shape of the solutions to these equations should look like, but I am not able to get the correct numerical shape. I am trying to solve the problem using the Shooting method. The code is the following:

R := 10
k := 0.1
h0 := 0.6
nlde1 = {f''[r] + 1/r f'[r] - k^2 f[r] (f[r]^2 - 1 + (a[r] - 1/(r*k))^2) == 0, f[R] == 1, f'[R] == 0};
nlde2 = {a''[r] + 1/r a'[r] - 1/r^2 a[r] - f[r]^2 (a[r] - 1/(r*k)) == 0, a[R] == 0, a'[R] == 0};
sol = NDSolve[{nlde1, nlde2}, {f, a}, {r, 0, R},Method -> "Shooting", "StartingInitialConditions" -> {f[0] == 0, f'[0] == 0.4, a[0] == 0, a[0] == h0}}]

From the physics of the problem I know that:

  • the order parameter $f$ must be $0$ at $r=0$ and $1$ at $r\rightarrow\infty$,
  • $f'> 0$ at $r=0$ and $f'= 0$ at $r\rightarrow\infty$,
  • the vector potential $a$ must be $0$ at $r=0$ and $0$ at $r\rightarrow\infty$,
  • $a' > 0$ at $r=0$ and $a'=0$ at $r\rightarrow\infty$.

The main problem I think comes from the fact that for $r=0$ both the equations have a singularity and NDSolve does not like that. Do you have some ideas on how to solve this problem?

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ What is GL? Is it Ginzburg-Landau? Or something else? $\endgroup$
    – Alexei Boulbitch
    Feb 17, 2016 at 9:22
  • $\begingroup$ Yes Ginzburg-Landau $\endgroup$
    – Mattia
    Feb 17, 2016 at 21:35
  • $\begingroup$ What does your boundary condition at r=0 describe then? What physics stays behind? $\endgroup$
    – Alexei Boulbitch
    Feb 18, 2016 at 8:36
  • $\begingroup$ The GL equations written in cylindrical coordinates are used to describe a single non interacting vortex. Both the order parameter and the vector potential go to zero at the center of a vortex. In my case such calculation is useful since I can calculate the lower critical field as a function of the GL parameter k. Analytical methods allows you to calculate Hc1(k) only in the limit k>>1 $\endgroup$
    – Mattia
    Feb 18, 2016 at 18:15
  • $\begingroup$ It sounds doubtful, since as soon as you put f(r=0)=0 as a boundary condition, you impose a strong limitation onto the system forcing such a structure. In contrast the Abrikosov vortexes should emerge spontaneously. Have a look into Landau&Lifshitz Statistical Physics pt. 2, Chapter 47. To me it sounds as a vortex that exists there inspite of the system due to some reason that is not described by the equations. $\endgroup$
    – Alexei Boulbitch
    Feb 19, 2016 at 8:29

1 Answer 1

Reset to default
5
$\begingroup$

The goal of this question, as I understand it, is the find a solution that vanishes at r = 0 and has a bounded, monotonic asymptotic solution at infinity. The first step is to identify the asymptotic solution. The asymptotic solution is chosen to eliminate nonlinearities:

a[r] = 1/(r*k)
f[r] = 1

With these as boundary conditions at large r and f[r0] = a]r0] = 0 at the origin, the ODE system becomes

R = 12; k = 1/10; r0 = 10^-5;
nlde1 = {f''[r] + 1/r f'[r] - k^2 f[r] (f[r]^2 - 1 + (a[r] - 1/(r*k))^2) == 0, 
    f[r0] == 0, f[R] == 1};
nlde2 = {a''[r] + 1/r a'[r] - 1/r^2 a[r] - f[r]^2 (a[r] - 1/(r*k)) == 0, 
    a[r0] == 0, a[R] == 1/(R k)};
sol = NDSolve[{nlde1, nlde2}, {f, a}, {r, r0, R}, Method -> { "Shooting", 
    "StartingInitialConditions" -> {f[r0] == 0, f'[r0] == 0.208, 
    a[r0] == 0, a'[r0] == 0.938}}]

The integration is started at r0 = 10^-5 to avoid the singularity in the equations at the origin. That this is a valid approximation can be seen by varying r0.

The result is,

Plot[{Evaluate[{f[r], a[r]} /. sol], 1/(k r)}, {r, r0, R}, AxesLabel -> {r, "f, a"}, 
    PlotRange -> {0, 2}, PlotStyle -> {Blue, Orange, Directive[Black, Dashed]}]

enter image description here

(The dashed curve, 1/(k r), is superimposed to show that a assumes its asymptotic form for r > 8.) Some experimentation with "StartingInitialConditions" is necessary to obtain a solution for R as large as 12.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you very much bbgodfrey! Exactly what I was looking for!! $\endgroup$
    – Mattia
    Feb 17, 2016 at 21:36
  • $\begingroup$ Do not want to bother you, but I am now using this code to find solutions for different values of k. My approach is to change k and acting on the initials conditions of the shooting method to find the analytical solution. The process is though very slow since the correct initial conditions have to be found point by point. Do you know a way to make the resolution more stable? I was trying with "StiffnessSwitching" as a sub method of "Shooting" buth seems nothing change. $\endgroup$
    – Mattia
    Feb 18, 2016 at 18:24
  • $\begingroup$ @Mattia Solving such problems as this usually is impossible. Fortunately, this one is merely very difficult. I obtained the solution in my answer by solving the equations for R = 2, determining a'[r0] and f[r0` from the answer and using it as a first guess for initial conditions for, say, R = 5 and proceeding this way until I got to R = 12.2 from which I could proceed no further. Not very satisfying but the best I could do. Good luck. $\endgroup$
    – bbgodfrey
    Feb 19, 2016 at 1:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.