Numerical Study of the Josephson Effect Using the Usadel Equation

Superconductors can carry electric current without resistance due to the formation of Cooper pairs. When two superconductors are separated by a thin normal metal, quantum interference effects cause a supercurrent to flow — even without an applied voltage. This phenomenon is known as the Josephson effect.

In this notebook, we simulate the Josephson effect by numerically solving the Usadel equation, which describes superconducting correlations in diffusive systems.

Goal

To understand how phase differences between superconductors give rise to supercurrents in a normal metal through the proximity effect.

Approach

• Reformulate the Usadel equation as a boundary value problem
• Implement a custom Runge-Kutta method for initial value problems
• Use scipy.solve_bvp to solve the full system numerically
• Compute the local density of states and the supercurrent
• Analyze how physical parameters affect the current-phase relationship

This group project was completed as part of the course TMA4320 – Scientific Computing at NTNU, Spring 2025.

Oppgave 1a)

vi har til å begynne med $\frac{d^2 y}{dx^2} = -4 \sin(2x), \quad y(0)=0, \quad y'(0)=2.$
Vi integrerer så med hensyn på y og får:

$\frac{d y}{dx} = -2 \cos(2x) + C_0$

etter å ha integrert enda en gang med hensyn på y får vi til slutt:

$y = \sin(2x) + C_0 x + C_1$

Vi setter deretter inn for initialbetingelser:

$\quad y(0)=0$ og $\quad y'(0)=2$ gir:

$\sin(0) + C_0 \cdot 0 + C_1 = 0$ og $2 \cos(0) + C_0 = 2$

dette gir $C_0 = 0$ og $C_1 = 0$

dermed blir løsningen $y = \sin(2x)$

Oppgave 1b)

Systemet vi starter med er: $\frac{d^2 y}{dx^2} = -4 \sin(2x), \quad y(0)=0, \quad y'(0)=2.$

vi starter med å definere $y_1 = y(x)$ og $y_2 = \frac{d}{dx} y$ fpr å gjøre dette om til et system av førsteordens differensiallikninger. Vi definerer videre vektoren y.
$\vec{y(x)} = \begin{pmatrix} y_1(x) \\ y_2(x) \end{pmatrix}$ Som vi kan derivere og gir $\frac{d}{dx} \vec{y} = \begin{pmatrix} y_1'(x) \\ y_2'(x) \end{pmatrix}$ da vi har $\frac{d^2 y}{dx^2} = -4 \sin(2x)$ kan vi skrive det hele som :

$\frac{d}{dx} \vec{y} = \begin{pmatrix} y(x) \\ -4\sin(2x) \end{pmatrix}$

vi uttrykker dette videre som:

$\frac{d}{dx} \vec{y} = A \cdot \vec{y(x)} + \vec{f(x)}$, der A er en matrise

Systemet vårt blir derfor som følger:

$\vec{y(x)} = \begin{pmatrix} y_1(x) \\ y_2(x) \end{pmatrix}$, $A = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$ og $\vec{f(x)} = \begin{pmatrix} 0 \\ -4\sin(2x) \end{pmatrix}$

Og med initialbetingelser: $\quad y(0)=0, \quad y'(0)=2.$ får vi:

$\vec{y(0)} = \begin{pmatrix} 0 \\ 2 \end{pmatrix}$

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp
from scipy.integrate import simpson
import time
import warnings
from IPython.display import display, Latex
warnings.filterwarnings("ignore")
import pandas as pd

plt.rcParams.update({
     "font.family": "serif",             
    "mathtext.fontset": "cm",
    "font.size": 12,                
    "axes.titlesize": 16,           
    "axes.labelsize": 16,           
    "xtick.labelsize": 10,          
    "ytick.labelsize": 10,          
    "legend.fontsize": 10,          
    "lines.linewidth": 2,           
    "axes.grid": True,              
    "grid.linestyle": "--",         
    "grid.alpha": 0.6, 
})

Oppgave 1c)

'''
Runge kutta funksjonen Funksjonen starter med initialbetingelsene, beregner de nødvendige k-verdiene for hvert steg, 
oppdaterer løsningen og justerer steglengden dynamisk for å holde den lokale feilen under den gitte toleransen. 
Resultatet er tilnæringer av y(x) og y'(x), og en oversikt over steglengden som funksjon av x. 

'''

tol = 10**(-7)
alpha = 0.8
h0 = 0.01
def y_anal(x):
    return np.sin(2 * x)

def y_doubleDerivative(x):
    return -4*np.sin(2*x)

def f(x,y):
    f1 = y[-1]
    f2 = y_doubleDerivative(x)
    return np.array([f1,f2])


def rungeKuttaMethod(x_init, x_end, y_init,h0, f,tol,alpha):
    h = h0
    n = 0
    n_tot = 0
    x_new = 0

    y_vals = [y_init]
    x_vals = [x_init]
    h_list = []
    k1 = np.array(f(x_vals[n],y_vals[n]))

    while (x_vals[n] < x_end):
        n_tot += 1
        h = min(h, x_end-x_vals[n])
        k2 = np.array(f((x_vals[n]+(1/2)*h),(y_vals[n]+((1/2)*h*k1))))
        k3 = np.array(f((x_vals[n]+(3/4)*h),(y_vals[n]+((3/4)*h*k2))))
        y_new = y_vals[n] +(1/9)*h*(2*k1+3*k2+4*k3)

        k4 = np.array(f((x_vals[n]+h),y_new))
        z_new = y_vals[n] + (1/24)*h*(7*k1+6*k2+8*k3+3*k4)

        est = np.linalg.norm(y_new-z_new)

        if est < tol:
            x_new = x_vals[n] + h
            n += 1
            y_vals.append(y_new)
            x_vals.append(x_new)
            k1 = k4
        
        h = alpha * h * (tol/est)**(1/3)
        h_list.append(h)

    return np.array(x_vals),np.array(y_vals),n, np.array(h_list), n_tot


x,y,n,h_list, n_tot = rungeKuttaMethod(0,2*np.pi,[0,2],h0,f,tol,alpha)

plt.plot(x,y[:,0], label = r"$y(x)$")
plt.plot(x,y[:,1], label = r"$y'(x)$")
plt.title(r"Løsning for $\frac{d^2 y}{dx^2} = -4 \sin(2x), \quad y(0)=0, \quad y'(0)=2.$ med Rk3(2)")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.show()

plt.plot(x,h_list)
plt.xlabel("X")
plt.ylabel("Steglengde")
plt.title("Steglengde som funksjon av x")
plt.show()

Steglengden tilpasser seg i forhold til hvordan hvor raskt løsningen endrer seg. Vi observerer at der løsningen endrer seg raskt vil feilen ofte bli større og som et resultat vil steglegnden bli mindre. Mens når løsningen endrer seg sakte blir lokal feil mindre og steglengden blir større.

Når funksjonen har et maksimum eller minimum er den førstederiverte nær null, men den andrederiverte er stor (funksjonen krummer mest her). Store verdier av den andrederiverte betyr at funksjonen endrer retning raskt, noe som krever en kortere steglengde. Steglengden er derfor minst ved ekstremalpunkter og størst ved nullpunkter til funksjonen.

oppgave 1d)

def analytisk(x):
    return np.sin(2*x)

tol_list = np.linspace(1e-1,1e-13,100)
avg_error_list = np.zeros(len(tol_list))

for i in range(len(tol_list)):
    x_vals,y_vals,y_deriv_vals,h_vals,n = rungeKuttaMethod(0,2*np.pi,[0,2],h0,f,tol_list[i],alpha)
    avg_error_list[i] = np.abs(y_vals[-1,0] - analytisk(x_vals[-1]))

plt.loglog(tol_list, avg_error_list)
plt.xlabel("Toleranse")
plt.ylabel("Feil")
plt.title("Feil som funksjon av toleranse")
plt.show()

Feilen øker som en funksjon av økt toleranse. Dette gir i stor grad mening da vi ved en høyere toleranse lar funksjonen vår ta "dårlige steg". slik vil man da ofte ta for "lange" steg og akkumulere seg en større feil. ved en lavere toleranse evaluerer vi enda flere verdier og tar mye flere steg, men da med en mye mindre feil som man kan se. Modellen tar da korte og "gode" steg.

alpha_list = np.linspace(0.5,0.95,100)
n_accept_low = []
n_denial_low = []
n_total_low = []
n_accept_high = []
n_denial_high = []
n_total_high = []

for i in alpha_list:
    x, y, n, h_list, n_tot = rungeKuttaMethod(0,2*np.pi,[0,2],h0,f,1e-7,i)
    n_accept_low.append(n)
    n_denial_low.append(n_tot - n)
    n_total_low.append(n_tot)
    x_high, y_high, n_high, h_list_high, n_tot_high = rungeKuttaMethod(0,2*np.pi,[0,2],h0,f,1e-1,i)
    n_accept_high.append(n_high)
    n_denial_high.append(n_tot_high - n_high)
    n_total_high.append(n_tot_high)

plt.plot(alpha_list, n_accept_low, label = "aksepterte steg")
plt.plot(alpha_list, n_denial_low, label = "forkastet steg")
plt.plot(alpha_list, n_total_low, label = "totale steg")
plt.legend()
plt.title("Antall steg som en funksjon av alpha med tol = 1e-7")
plt.ylabel("Antall steg")
plt.xlabel("Alpha")
plt.grid()
plt.show()

plt.plot(alpha_list, n_accept_high, label = "aksepterte steg")
plt.plot(alpha_list, n_denial_high, label = "forkastet steg")
plt.plot(alpha_list, n_total_high, label = "totale steg")
plt.legend()
plt.title("Antall steg som en funksjon av alpha med tol = 1e-1")
plt.ylabel("Antall steg")
plt.xlabel("Alpha")
plt.grid()
plt.show()

Her har vi plottet antall steg som en funksjon av alfa, men med forskjellige toleranser. Den øverste med toleranse på $e-7$ og den andre på $e-1$. Som vi kan se på den første ser vi at en lav toleranse tvinger modellen til å ta gode steg og at alfa derfor ikke påvirker om noen steg blir forkastet. Men når vi øker toleransen ser vi til slutt at høyere alfa verdier bidrar til at steg blir forkastet da en høy toleranse godtar at modellen vår tar "dårlige" steg.

Vi ser jo også selvførgelig at antall steg totalt går ned med en høyere alpha verdi. Dette gir veldig mening da steglengden er direkte proporsjonal med alpha

Oppgave 1e)

def secant_method(g,z0,z1, tol):
    z =[z0,z1]
    z.append((z[-2]*g(z[-1])-z[-1]*g(z[-2]))/(g(z[-1])-g(z[-2])))

    while (abs(z[-2] - z[-1]) > tol):
        z.append((z[-2]*g(z[-1])-z[-1]*g(z[-2]))/(g(z[-1])-g(z[-2])))
    return z


def funksjon(x):
    return  x + np.sin(x)  + np.cos(x)

rot = secant_method(funksjon, 2,3, 10**-7)

#tar ikke med de initielle gjettene, plotter fra verdi 2 og utover, som er de første beregnede gjettene.

plt.plot(np.arange(len(rot)-2), rot[2:], label= f'Root: {np.round(rot[-1],5)}')
plt.title("Sekant Metoden")
plt.ylabel("Beregnet rot")
plt.xlabel("Antall Iterasjoner")
plt.legend()
plt.show()

Plottet viser hvordan vi starter med to initiale gjetninger, og hvordan hver iterasjon gir en ny tilnærming til roten. Metoden fortsetter inntil differansen mellom to påfølgende gjetninger er mindre enn en gitt toleranse. Vi kan observere at etter åtte iterasjoner konvergerer metoden og stabiliseres.

beregnet_rot = rot[-1]
beregnet_verdi = funksjon(beregnet_rot)
print(f'Verdi for funksjon ved å sette inn beregnet rot')
print(f'y = {beregnet_verdi}')
print("------------------------------------")

x_verd = np.linspace(-5,10,100)
y_verd = []
for i in x_verd:
    y_verd.append(funksjon(i))
plt.plot(x_verd,y_verd, label = "Y = x + sin(x)  + cos(x)")
plt.axhline(y = 0, color = "green", linestyle = "--", label = " Y = 0")
plt.axvline(beregnet_rot, color = "r", linestyle = "--", label = f"X = {np.round(beregnet_rot,3)}")
plt.title("sekantmetoden for Y = x + sin(x)  + cos(x)")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.show()

Verdi for funksjon ved å sette inn beregnet rot
y = -1.1102230246251565e-16
------------------------------------

Vi setter inn den beregnede roten fra sekantmetoden, og ser at resultatet er svært nært null. Siden vi ikke kan finne det eksakte nullpunktet algebraisk, må vi nøye oss med en approksimasjon. Det viser seg imidlertid at denne tilnærmingen er svært god, og sekantmetoden ser ut til å fungere bra. I tillegg kan vi også se beregningen grafisk.

oppgave 1f)

'''
Koden finner den riktige initialverdien for et randverdi-problem. 
Dette ved hjelp av en kombinasjon av Runge-Kutta og sekantmetoden, 
og visualiserer deretter løsningen for hver iterasjon.
'''

tol = 1e-7

def g(b):
    rk = rungeKuttaMethod(0,2*np.pi,[0,b],h0,f,tol,alpha)
    return rk[1][-1][0]

b_vals = secant_method(g,1,4,tol)

for i in range(len(b_vals)):
    x_vals,y_vals,y_deriv_vals,h_vals,n = rungeKuttaMethod(0,2*np.pi,[0,b_vals[i]],h0,f,tol,alpha)
    plt.plot(x_vals,y_vals[:,0], label = f'b nr. {i+1}: {np.round(b_vals[i],3)}')

plt.title("Løsning etter hver iterasjon av rotfunksjon")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.show()

fra plottet ser vi at hvis vi velger feil b-verdi så vil ikke løsningen stemme, men etterhvert som sekantmetoden finner rett b-verdi tilsvarer dette den samme løsningen vi fant både analytisk og ved rungekutta metoden.

oppgave 1g)

'''
Koden finner den riktige initialverdien for et randverdi-problem. 
Dette ved hjelp av en kombinasjon av Runge-Kutta og sekantmetoden.
Plotter løsning for riktig b verdi.
'''
tol = 1e-7

def ydd(x):
    return np.sin(x)

'her skriver vi om 2. ordenslikninga om til 1. orden, som i oppgave 1b). metoden er den samme'
def f(x,y):
    f1 = y[1]
    f2 = y[0] + ydd(x)
    return np.array([f1,f2])

def g(b):
    rk = rungeKuttaMethod(0,12,[0,b],h0,f,tol,alpha)
    return rk[1][-1][0]

start_time_rk32 = time.time()
b = secant_method(g,-1,1,tol)
b_best = b[-1]

rk1 = rungeKuttaMethod(0,12,[0,b_best], h0, f, tol, alpha)
end_time_rk32 = time.time()
time_rk32 = end_time_rk32 - start_time_rk32

plt.plot(rk1[0], rk1[1][:,0], color = "darkred")
plt.title(r"løsning for  $ \frac{d^2 y}{dx^2} = y(x) + sin(x)$")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

Ved hjelp av å løse runge kutta for ulike b verdier og deretter bruke sekantmetoden på randbetingelsene kommer vi fram til den riktige løsningen.

Oppgave 1h)

def boundary_conditions(y_l, y_r):
    return np.array((y_l[0],y_r[0]))

x_vals = rk1[0]
y_init = np.zeros((2, x_vals.size))
start_time_scipy = time.time()
sol = solve_bvp(f,boundary_conditions,x_vals, y_init)
end_time_scipy = time.time()
time_scipy = end_time_scipy - start_time_scipy

plt.plot(sol.x, sol.y[0], label = r"$scipy$")
plt.plot(x_vals, rk1[1][:,0],linestyle = "--", label = r"$RK3(2)$")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Løsning med RK3(2) og Scipy")
plt.legend(loc = "upper left")
plt.show()
print("---------------------------")
print(f"Tid RK32: {np.round(time_rk32,2)} s")
print(f"Tid Scipy: {np.round(time_scipy,2)} s")
print("---------------------------")

---------------------------
Tid RK32: 2.93 s
Tid Scipy: 0.68 s
---------------------------

Begge metodene løser ligningen tilfredsstillende innenfor den gitte toleransen. Samtidig ser vi at SciPy-metoden er betydelig raskere – nesten fem ganger raskere i dette tilfellet. Derfor vil løsningen med SciPy være å foretrekke, og ved mer komplekse problemer kan denne tidsforskjellen bli enda mer markant.

error = abs(rk1[1][:,0]-sol.y[0])
plt.plot(x_vals,error)
plt.xlabel("x")
plt.ylabel("Differanse")
plt.title("Differanse mellom RK3(2) og SciPy som funksjon av x")
plt.show()

Ut fra plottet ser vi at de to numeriske metodene gir nesten identiske verdier ved start- og sluttpunktene, noe som er forventet siden begge benytter de samme betingelsene. Videre er løsningene svært like i de fleste ekstremalpunkter, med unntak av det siste, der differansen er størst. For øvrig er metodene mest forskjellige mellom ekstremalpunktene.

Oppgave 2a)

'''
Tar inn en kompleks matrise og returnerer denne som en reel vektor
Og motsatte funksjon som tar inn en reel vektor og returerer en kompleks matrise
'''
def flatten_complex_matrix(matrix):
    matrix = matrix.flatten()
    matrix_real = np.real(matrix)
    matrix_imag = np.imag(matrix)
    matrix = np.concatenate((matrix_real,matrix_imag))
    return matrix
def unflatten_complex_matrix(matrix):
    N = int(len(matrix)//2)
    matrix_real = matrix[:N]
    matrix_imag = matrix[N:]

    complex_matrix = matrix_real + 1j * matrix_imag
    return complex_matrix.reshape((2,2))

Validering av funksjoner:

matrise = np.array([[1 + 1j, 2 + 2j],
              [3 + 3j, 4 + 4j]])
print(matrise)
print("--------------------------")
matrise_flat = flatten_complex_matrix(matrise)

print(matrise_flat)
print("--------------------------")
matrise_original = unflatten_complex_matrix(matrise_flat)

print(matrise_original)

[[1.+1.j 2.+2.j]
 [3.+3.j 4.+4.j]]
--------------------------
[1. 2. 3. 4. 1. 2. 3. 4.]
--------------------------
[[1.+1.j 2.+2.j]
 [3.+3.j 4.+4.j]]

Oppgave 2b)

'''
Funksjon som tar inn 4 vektorer og deretter setter disse sammen til en lang vektor
Andre funksjon gjør motsatte prosessen ved slicing
'''
def make_long_vec(m1,m2,m3,m4):
    vector = np.concatenate((m1,m2,m3,m4))
    return vector

def unmake_long_vec(vector):
    m1 = vector[0:8]
    m2 = vector[8:16]
    m3 = vector[16:24]
    m4 = vector[24:32]
    return m1,m2,m3,m4

Validering av funksjoner:

m1 = np.array([1,2,3,4,5,6,7,8])
m2 = np.array([9,10,11,12,13,14,15,16])
m3 = np.array([17,18,19,20,21,22,23,24])
m4 = np.array([25,26,27,28,29,30,31,32])
print(m1)
print(m2)
print(m3)
print(m4)
print("----------------------------------------------------------------")
long_vec = make_long_vec(m1,m2,m3,m4)
print(long_vec)
print("----------------------------------------------------------------")
m1,m2,m3,m4 = unmake_long_vec(long_vec)
print(m1)
print(m2)
print(m3)
print(m4)

[1 2 3 4 5 6 7 8]
[ 9 10 11 12 13 14 15 16]
[17 18 19 20 21 22 23 24]
[25 26 27 28 29 30 31 32]
----------------------------------------------------------------
[ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
 25 26 27 28 29 30 31 32]
----------------------------------------------------------------
[1 2 3 4 5 6 7 8]
[ 9 10 11 12 13 14 15 16]
[17 18 19 20 21 22 23 24]
[25 26 27 28 29 30 31 32]

Oppgave 2c)

'''
Vi tar inn 4  2x2 matriser. Vi gjør om hver av matrisene til en vektor med lengde 8. 
Til slutt setter vi sammen de 4 vektorene til en lang vektor med lengde 32.

På den  mosatte funksjonen tar vi inn en lang vektor med lengde 32
Denne splitter vi til 4 mindre vektorer med lengde 8 som settes sammen til hver sin 2x2 matrise. 
'''

def make_unknown_vector(matrix1,matrix2,matrix3,matrix4):
    vec1 = flatten_complex_matrix(matrix1)
    vec2 = flatten_complex_matrix(matrix2)
    vec3 = flatten_complex_matrix(matrix3)
    vec4 = flatten_complex_matrix(matrix4)

    unknown_vector = make_long_vec(vec1,vec2,vec3,vec4)
    return unknown_vector

def unmake_unknown_vector(unknown_vector):
    vec1,vec2,vec3,vec4 = unmake_long_vec(unknown_vector)

    matrix1 = unflatten_complex_matrix(vec1)
    matrix2 = unflatten_complex_matrix(vec2)
    matrix3 = unflatten_complex_matrix(vec3)
    matrix4 = unflatten_complex_matrix(vec4)

    return matrix1,matrix2,matrix3,matrix4

Validering:

m1 = np.array([[1+1j,2+2j],[3+3j,4+4j]])
m2 = np.array([[5+5j,6+6j],[7+7j,8+8j]])
m3 = np.array([[9+9j,10+10j],[11+11j,12+12j]])
m4 = np.array([[13+13j,14+14j],[15+15j,16+16j]])
print(m1)
print("---")
print(m2)
print("---")
print(m3)
print("---")
print(m4)
print("----------------------------")
unknown_vector = make_unknown_vector(m1,m2,m3,m4)
print(unknown_vector)
print("----------------------------")
m5,m6,m7,m8 = unmake_unknown_vector(unknown_vector)
print(m5)
print("---")
print(m6)
print("---")
print(m7)
print("---")
print(m8)

[[1.+1.j 2.+2.j]
 [3.+3.j 4.+4.j]]
---
[[5.+5.j 6.+6.j]
 [7.+7.j 8.+8.j]]
---
[[ 9. +9.j 10.+10.j]
 [11.+11.j 12.+12.j]]
---
[[13.+13.j 14.+14.j]
 [15.+15.j 16.+16.j]]
----------------------------
[ 1.  2.  3.  4.  1.  2.  3.  4.  5.  6.  7.  8.  5.  6.  7.  8.  9. 10.
 11. 12.  9. 10. 11. 12. 13. 14. 15. 16. 13. 14. 15. 16.]
----------------------------
[[1.+1.j 2.+2.j]
 [3.+3.j 4.+4.j]]
---
[[5.+5.j 6.+6.j]
 [7.+7.j 8.+8.j]]
---
[[ 9. +9.j 10.+10.j]
 [11.+11.j 12.+12.j]]
---
[[13.+13.j 14.+14.j]
 [15.+15.j 16.+16.j]]

oppgave 2d)

'''
funksjon som tar inn 32 element lang vektor, gjør det om til matriser og regner ut deriverte. 
Før den returneres til vektor form igjen.
'''
def derivate_vec(vec,eps):
    gamma,gamma_tilde,omega,omega_tilde = unmake_unknown_vector(vec)

    I = np.eye(2)
    N = np.linalg.inv((I-(gamma@gamma_tilde)))
    N_tilde = np.linalg.inv((I-(gamma_tilde@gamma)))

    deriv_gamma = omega
    deriv_gamma_tilde = omega_tilde

    deriv_omega = -2j*complex(eps,delta)*gamma - (2*omega@N_tilde@gamma_tilde@omega)
    deriv_omega_tilde = -2j*complex(eps,delta)*gamma_tilde - (2*omega_tilde@N@gamma@omega_tilde)

    deriv_vec = make_unknown_vector(deriv_gamma,deriv_gamma_tilde,deriv_omega,deriv_omega_tilde)
    return deriv_vec

oppgave 2e)

'''
Funksjon som tar inn en vektor x  med m elementer og en  32xm matrise
Den returnerer den deriverte av vektoren i hver posisjon x. 
Den returerer derfor en 32 x m matrise
'''
def make_deriv_matrix(x_pos,matrix):
    new_matrix = np.zeros_like(matrix)
    for i in range(len(x_pos)):
        new_matrix[:,i] = derivate_vec(matrix[:,i],eps)
    return new_matrix

oppgave 2f)

'''
Funksjonen tar to 32-komponent vektorer (v_left og v_right). De 
inneholder Riccati-matrisene og deres deriverte ved venstre og høyre kant av den normale metallen.
Den beregner venstresidene og returnerer disse som en enkelt reell 32-komponent vektor, slik at den er kompatibel med solve_bvp.
'''
def boundry_condition(v_left,v_right):
    gamma_left,gamma_tilde_left,omega_left,omega_tilde_left = unmake_unknown_vector(v_left)
    gamma_right,gamma_tilde_right,omega_right,omega_tilde_right = unmake_unknown_vector(v_right)

    N_left = np.linalg.inv((np.eye(2)-(gamma_l @ gamma_t_l)))
    N_tilde_left = np.linalg.inv((np.eye(2)-(gamma_t_l @ gamma_l)))
    N_right = np.linalg.inv((np.eye(2)-(gamma_r @ gamma_t_r)))
    N_tilde_right = np.linalg.inv((np.eye(2)-(gamma_t_r @ gamma_r)))

    left_riccati = omega_left + (1/(zeta*l)*(np.eye(2)-(gamma_left@gamma_t_l))) @ N_left @ (gamma_l-gamma_left)
    left_riccati_tilde = omega_tilde_left + (1/(zeta*l)*(np.eye(2)-(gamma_tilde_left@gamma_l))) @ N_tilde_left @ (gamma_t_l-gamma_tilde_left)
    right_riccati = omega_right - (1/(zeta*l)*(np.eye(2)-(gamma_right@gamma_t_r))) @ N_right @ (gamma_r-gamma_right)
    right_riccati_tilde = omega_tilde_right - (1/(zeta*l)*(np.eye(2)-gamma_tilde_right@gamma_r))@N_tilde_right @ (gamma_t_r-gamma_tilde_right)

    vec = make_unknown_vector(left_riccati,left_riccati_tilde,right_riccati,right_riccati_tilde)
    return vec

Funksjonen som beregner randbetingelsene returnerer en vektor med 32 elementer. Denne er kompatibel med matrisen som returneres av funksjonen som beregner de deriverte for alle x-verdiene, ettersom hver kolonne i denne matrisen har samme lengde som randbetingelsesvektoren.

Oppgave 2g)

'''
sjekk_lik_null er en funksjon som sjekker om det er numerisk lik 0
Det kjøres en simulering for et normalt metall imellom to normale metaller
Det forventes at det ikke skal skje noe interessant
Simuleringen kjøres med solve_bvp
'''
def sjekk_lik_null(sol):
    if np.allclose(sol, 0, atol=0):
        print("Alle komponentene er numerisk null.")
    else:
        print("Det finnes komponenter som ikke er null.")

zeta = 3
delta = 0.01
m = 101
l = 1

gamma_l = np.zeros((2,2))
gamma_r = np.zeros((2,2))
gamma_t_l = np.zeros((2,2))
gamma_t_r = np.zeros((2,2))

x = np.linspace(0,1,m)
y = np.zeros((32,m))

eps = 0
sol_1 = solve_bvp(make_deriv_matrix,boundry_condition,x,y)

eps = 1
sol_2 = solve_bvp(make_deriv_matrix,boundry_condition,x,y)

eps = 2
sol_3 = solve_bvp(make_deriv_matrix,boundry_condition,x,y)

print(r"Sjekk for $\epsilon = 0$")
print(sol_1.y)
sjekk_lik_null(sol_1.y)
print("--------------------------------------")
print(r"Sjekk for $\epsilon = 1$")
print(sol_2.y)
sjekk_lik_null(sol_2.y)
print("--------------------------------------")
print(r"Sjekk for $\epsilon = 2$")
print(sol_3.y)
sjekk_lik_null(sol_3.y)
print("--------------------------------------")

Sjekk for $\epsilon = 0$
[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]]
Alle komponentene er numerisk null.
--------------------------------------
Sjekk for $\epsilon = 1$
[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]]
Alle komponentene er numerisk null.
--------------------------------------
Sjekk for $\epsilon = 2$
[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]]
Alle komponentene er numerisk null.
--------------------------------------

Oppgave 2h)

'''
Tettheten av tilstander er beregnet ved hjelp av greens function og density function.
Der definisjoner er gitt i oppgavetekst
 '''
def green_function(x,vec):
    gamma,gamma_tilde,omega,omega_tilde = unmake_unknown_vector(vec)
    N = np.linalg.inv((np.eye(2)-gamma@gamma_tilde))
    N_tilde = np.linalg.inv((np.eye(2)-gamma_tilde@gamma))

    green = np.block([[2*N-np.eye(2),2*N@gamma],[-2*N_tilde@gamma_tilde,-2*N_tilde+np.eye(2)]])

    return green

def density(x,solution):
    density_list = np.zeros(len(x))
    p_hat_three = np.block([[np.eye(2),np.zeros((2,2))],[np.zeros((2,2)),-np.eye(2)]])

    for i in range(len(x)):
        density_list[i] = np.real(np.trace(p_hat_three@green_function(x[i],solution.y[:,i])))/4

    return density_list

density_1 = density(x,sol_1)
density_2 = density(x,sol_2)
density_3 = density(x,sol_3)

plt.plot(x,density_1,label="Eps = 0")
plt.title("Tetthet av tilstander for $eps = 0$")
plt.xlabel("X")
plt.ylabel("Tetthet av tilstander")
plt.legend()
plt.show()
plt.plot(x,density_2,label="Eps = 1")
plt.title("Tetthet av tilstander for $eps = 1$")
plt.xlabel("X")
plt.ylabel("Tetthet av tilstander")
plt.legend()
plt.show()
plt.plot(x,density_3,label="Eps = 2")
plt.title("Tetthet av tilstander for $eps = 2$")
plt.xlabel("X")
plt.ylabel("Tetthet av tilstander")
plt.legend()
plt.show()

Vi har et system med et vanlig metall mellom to andre vanlige metaller. Da forventes det ingen variasjon i den elektroniske tilstandstettheten. Den vil være konstant og når vi normaliserer med $D_0$ vil den være lik 1.

Oppgave 2i)

'''
Gjentagelse av forrige oppgave men nå med en superconductor på hver side
Grensebetingelene er som i oppgaveteksten og er definert under.
'''
def bold_v(sigma,eps):
    return np.arctanh(sigma/complex(eps,delta))

def s(sigma,eps):
    return np.sinh(bold_v(sigma,eps))

def c(sigma,eps):
    return np.cosh(bold_v(sigma,eps))

def make_bc(sigma,eps,phi_L,phi_R):
    gamma_l = np.block([[0,s(sigma,eps)/(1+c(sigma,eps))],[s(-sigma,eps)/(1+c(-sigma,eps)),0]]) * np.e**(1j*phi_L)
    gamma_r = np.block([[0,s(sigma,eps)/(1+c(sigma,eps))],[s(-sigma,eps)/(1+c(-sigma,eps)),0]]) * np.e**(1j*phi_R)
    gamma_t_l = np.block([[0,s(-sigma,eps)/(1+c(-sigma,eps))],[s(sigma,eps)/(1+c(sigma,eps)),0]]) * np.e**(-1j*phi_L)
    gamma_t_r = np.block([[0,s(-sigma,eps)/(1+c(-sigma,eps))],[s(sigma,eps)/(1+c(sigma,eps)),0]]) * np.e**(-1j*phi_R)

    return gamma_l,gamma_r,gamma_t_l,gamma_t_r

delta=0.01
l = 1
m = 101
eps = 2
x = np.linspace(0,l,m)
y = np.zeros((32,m))
gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,0,0) #her phi_r og phi_l = 0 og sigma = +-1
sol_4 = solve_bvp(make_deriv_matrix,boundry_condition,x,y)

Oppgave 2j)

dens = density(x,sol_4)
plt.plot(x,dens)
plt.title("Tetthet av tilstander")
plt.xlabel("X")
plt.ylabel("Tetthet av tilstander")
plt.axvline(x = 0, color = "green", linestyle = "--")
plt.axvline(x = 1, color = "green", linestyle = "--")
plt.show()

Vi observerte at når randbetingelsene ble satt til null, var tettheten av tilstander (density of states) jevnt fordelt. Når vi deretter introduserte ikke-trivielle randbetingelser for superlederen, ser vi at superlederens egenskaper påvirker metallet i midten, noe som fører til en økt tetthet av tilstander nær grenseflatene mot superlederne. Dette er i samsvar med forventningene fra figur 1, som viser at antallet tilgjengelige tilstander er høyest nær grenseflatene og lavest i midten av metallet.

Oppgave 2k)

'''
Her tetthet av tilstander for forskjellige verdier av epsilon, evaluert i miden av det vanlige metallet
'''
m = 101
l = 1
eps_list = np.linspace(2,0,m)
x = np.linspace(0,l,m)
y = np.zeros((32,m))
density_list1 = np.zeros(m)
density_list2 = np.zeros(m)
density_list3 = np.zeros(m)

for i in range(m):
    eps = eps_list[i]
    gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,0,0)
    sol_1 = solve_bvp(make_deriv_matrix,boundry_condition,x,y, max_nodes = m)
    density_list1[i] = density(x,sol_1)[int(m/2)]
    y = sol_1.y

l = 0.5
y = np.zeros((32,m))
x = np.linspace(0,l,m)
for i in range(m):
    eps = eps_list[i]
    gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,0,0)
    sol_2 = solve_bvp(make_deriv_matrix,boundry_condition,x,y, max_nodes = m)
    density_list2[i] = density(x,sol_2)[int(m/2)]
    y = sol_2.y

l = 2
y = np.zeros((32,m))
x = np.linspace(0,l,m)
for i in range(m):
    eps = eps_list[i]
    gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,0,0)
    sol_3 = solve_bvp(make_deriv_matrix,boundry_condition,x,y, max_nodes = m)
    density_list3[i] = density(x,sol_3)[int(m/2)]
    y = sol_3.y

plt.plot(eps_list,density_list1,label="l = 1")
plt.plot(eps_list,density_list2,label="l = 0.5")
plt.plot(eps_list,density_list3,label="l = 2")
plt.title("Tetthet av tilstander i midten av metallet")
plt.xlabel(r"$\epsilon$")
plt.ylabel("Tetthet av tilstander")
plt.legend()
plt.show()

For kortere lengder så ser vi et betydelig minigap ved mindre energier. Dette viser til at superlederegenskapene "lekker" inn i det vanlige metallet. Det kan gi mening at dette gapet er størst for kortere lengder da større del av egenskapene "lekker" dypere inn i metallet. Vi ser at gapet minker i tråd med økning i l. I tillegg blir også gapkantene mindre markante. Dette understreker videre at superlederegenskapene i mindre grad når midten av metallet ved økende lengde. vi ser også at når epsilon blir større at tilstandstettheten går tilbake til normalen da den går mot 1 uavhengig av lengde.

oppgave 2l)

def integrand(vec):
    gamma,gamma_tilde,omega,omega_tilde = unmake_unknown_vector(vec)
    N = np.linalg.inv((np.eye(2)-(gamma@gamma_tilde)))
    N_tilde = np.linalg.inv((np.eye(2)-(gamma_tilde@gamma)))

    deriv_N = N @ (omega @ gamma_tilde + gamma @ omega_tilde) @ N
    deriv_N_tilde = N_tilde @ (omega_tilde @ gamma + gamma_tilde @ omega) @ N_tilde

    deriv_green = 2*np.block([[deriv_N,(N@omega)+(deriv_N@gamma)],[(-1*N_tilde@omega_tilde)-(deriv_N_tilde@gamma_tilde),-deriv_N_tilde]])
    green = np.block([[2*N-np.eye(2),2*N@gamma],[-2*N_tilde@gamma_tilde,-2*N_tilde+np.eye(2)]])

    p_hat_three = np.block([[np.eye(2),np.zeros((2,2))],[np.zeros((2,2)),-np.eye(2)]]) 
    integrand = np.real(np.trace(p_hat_three @ (green @ deriv_green - deriv_green @ green)))
    return integrand

l = 1
m = 101
eps_list = np.array([2,1.5,1,0.5,0])
integrand_list = np.zeros(m)

for j in range(len(eps_list)):
    y = np.zeros((32,m))
    x = np.linspace(0,l,m)
    eps = eps_list[j]
    gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,0,0)
    sol = solve_bvp(make_deriv_matrix,boundry_condition,x,y)
    for i in range(m):
        integrand_list[i] = integrand(sol.y[:,i])
    
    plt.plot(x,integrand_list,label=f"Eps = {eps}")
    plt.title("Strøm integrand for forskjellige eps verdier")
    plt.ylabel("Strøm")
    plt.xlabel("Tversnitt")
    plt.legend()
    integrand_list = np.zeros(m)
plt.show()

Strøm integranden for forskjellige epsilon verdier er tilnærmet lik null da det er i størrelsesordener på $e-14$ og mindre. Ifølge teoridelen er strøm integranden relatert til faseforskjell. Da denne er null vil naturligvis strøm integranden også være tilnærmet lik eller lik 0.

oppgave 2m)

l = 1
m = 101

eps_list = np.linspace(2,0,m)
integrand_list = np.zeros(m)
for i in range(m):
    eps = eps_list[i]
    gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,1,0)
    sol = solve_bvp(make_deriv_matrix,boundry_condition,x,y, max_nodes = m)
    integrand_list[i] = integrand(sol.y[:,int(m/2)])
    y = sol.y

plt.plot(eps_list,integrand_list)
plt.title(r"Strøm integrand som funksjon av energi med $\Delta \phi = 1$")
plt.xlabel(r"$\epsilon$")
plt.ylabel("Strøm integrand")
plt.xlabel(r"$\epsilon$")
plt.show()

Figuren viser hvordan bidraget til strømintegranden varierer med energien $\epsilon$ når faseforskjellen $\Delta \phi = 1$. Vi ser at integranden skifter fortegn, noe som indikerer at ulike energier gir forskjellig bidrag til strømmen.

l = 1
m = 101
eps_list = np.array([2,1.5,1,0.5,0])
integrand_list = np.zeros(m)

for j in range(len(eps_list)):
    y = np.zeros((32,m))
    x = np.linspace(0,l,m)
    eps = eps_list[j]
    gamma_l,gamma_r,gamma_t_l,gamma_t_r = make_bc(1,eps,1,0)
    sol = solve_bvp(make_deriv_matrix,boundry_condition,x,y)
    for i in range(m):
        integrand_list[i] = integrand(sol.y[:,i])
    diff = abs(abs(max(integrand_list)) - abs(min(integrand_list)))
    print(rf"Absolutt differanse med eps = {eps}")
    print(diff)
    print("---------------------------------------")
    plt.plot(x,integrand_list,label=f"Eps = {eps}")
    plt.title("Strøm integrand for forskjellige eps verdier")
    plt.ylabel("Strøm")
    plt.xlabel("Tversnitt")
    plt.legend()
    integrand_list = np.zeros(m)
plt.show()

Absolutt differanse med eps = 2.0
1.02292954831662e-06
---------------------------------------
Absolutt differanse med eps = 1.5
4.11123583926587e-07
---------------------------------------
Absolutt differanse med eps = 1.0
0.00019115650854573119
---------------------------------------
Absolutt differanse med eps = 0.5
8.974326382027442e-06
---------------------------------------
Absolutt differanse med eps = 0.0
3.626307587545341e-12
---------------------------------------

I oppgave 2l hadde vi ingen faseforskjell mellom superlederne. Nå har vi en faseforskjell på 1, derfor induseres det en strøm imellom superlederene. Vi kan også konkludere med at strøm integraden er tilnærmet konservativ. Dette er understreket ved at maksimal differanse mellom strøm integranden er svært liten over tversnittet av metallet. Dette er printet ut og er visulisert for 5 forskjellige epsilon verdier ovenfor.

oppgave 2n)

m = 101
phase_diff = np.linspace(0, 2*np.pi, 10)
eps_list = np.linspace(2, 0, m)
result_list = np.zeros(len(phase_diff))
l = 1

for j, phi in enumerate(phase_diff):
    integrand_vals = np.zeros(m)
    
    for i, eps in enumerate(eps_list):
        gamma_l, gamma_r, gamma_t_l, gamma_t_r = make_bc(1, eps, phi, 0)
        

        x_vals = np.linspace(0, l, m)
        y_init = np.zeros((32, m))
        sol = solve_bvp(make_deriv_matrix, boundry_condition, x_vals, y_init)
        
        integrand_values = np.array([integrand(sol.y[:, k]) for k in range(m)])
        
        integrand_vals[i] = - simpson(integrand_values, eps_list)
    result_list[j] = simpson(integrand_vals, eps_list)

plt.plot(phase_diff, result_list, 'o-')
plt.xlabel(r"Faseforskjell $\Delta \phi$")
plt.ylabel(r"$I/I_0$")
plt.title("Strøm som funksjon av faseforskjell")
plt.tight_layout()
plt.show()

data = []
for i, strom in enumerate(result_list):
    fasevinkel = np.round((phase_diff[i]), 2)
    strom_rundet = np.round(strom, 2)
    data.append([fasevinkel, strom_rundet])

data = pd.DataFrame(data, columns = ["Fasevinkel", "Strøm"])

data

	Fasevinkel	Strøm
0	0.00	0.00
1	0.70	0.76
2	1.40	1.27
3	2.09	1.32
4	2.79	0.72
5	3.49	-0.72
6	4.19	-1.32
7	4.89	-1.27
8	5.59	-0.76
9	6.28	0.00

Vi observerer at når vi plotter $\frac{I}{I_0}$ som funksjon av faseforskjellen $\Delta \phi$ følger den en sinuskurve.

Vi observerer at strømmen er null ved $\Delta \phi = 0, \pi, 2\pi$ og maksimal ved $\Delta\phi = \frac{\pi}{2}, \frac{3\pi}{2}$

Dette kan tyde på at strømmen er proporsjonal med $\sin(\Delta \phi)$

Fra oppgaven er strømmen $\frac{I}{I_0}$ definert som det negative integralet av strømintegranden over epsilon [0,2], fra plottet ser vi at for $\sin(\Delta \phi = 1)$ er strømmen ca lik 1. Dette gir også mening fra det vi plottet i oppgave 2m), Der vi kan se at arealet under kurven er negativ og ser ut til å være lik ca samme verdi. Med fortegnet fra definisjonen stemmer dette godt.