선생님의 글
@dc1394
라고 안내했고, 실제 공간차분법으로 줄리아로 풀어봤다.

Juria 버전

Julia 1.5.3

미분 방정식

구해야 할 미분 방정식은 수소 원자에 대한 Schr≠dinger 방정식이다

\left(-\frac{1}{2} \frac{d^2}{dr^2}-\frac{1}{r}\frac{d}{dr} -\frac{1}{r} \right) \psi(r) = E \psi(r)

의 지름 좌표의 미분 방정식입니다.경계 조건은

\psi(r_c) = 0

.

미분산자의 차분화

그러면 미분산자를 차분한다.
$f(r i+a) 및 $f(r i-a)$r주변의 테일러 전개:

f(r_i+a) = f(r_i) + \frac{df}{dr}|_{r = r_i}  a + \frac{1}{2} \frac{d^2f}{dr^2}|_{r = r_i}  a^2 +\cdots

f(r_i-a) = f(r_i) -\frac{df}{dr}|_{r=r_i}  a + \frac{1}{2} \frac{d^2f}{dr^2}|_{r = r_i}   a^2 + \cdots

1층은 미분이다

\frac{df}{dr}|_{r = r_i} = \frac{f(r_i+a) -f(r_i-a) }{2a}

2 단계 미분

\frac{d^2f}{dr^2}|_{r = r_i} = \frac{f(r_i+a) -2f(r_o)+ f(r_i-a)}{a^2}

.
공간 차이는 $i=1,2,\cdots,N달러

r_i = a i - \frac{a}{2}

.시스템 반지름은 $R=aN달러입니다.이 때 첫 번째 점의 좌표는 $r입니다.1=a/2$, 마지막 점의 좌표는 $rN=a-a/2달러입니다.그럼 $i=1달러일 때 1층 미분은

\frac{df}{dr}|_{r = a/2} = \frac{f(3a/2) -f(-a/2) }{2a}

단, $f(-a/2)달러는 정의되지 않았습니다.단, 달러 r<0$의 구역은 회전 대칭성을 가지기 때문에 $r$r와 같은 값이어야 합니다.즉, $f(-a/2)=f(a/2)달러.
$i=N달러일 경우 $f(aN+a/2)달러의 값이 필요합니다.현재 디스크 영역을 고려하고 있기 때문에 $r>R달러가 0으로 풀리면 $f(aN+a/2)=0달러가 훨씬 넘을 수 있습니다.
1 단계 미분과 2 단계 미분

\frac{df}{dr}|_{r = r_i} = \sum_j c_{ij} f(r_j)

\frac{d^2f}{dr^2}|_{r = r_i} = \sum_j d_{ij} f(r_j)

... 면
계수 $c{ij} 달러 및 $d{ij}달러가 $i=1달러일 때
$$
c_{1j} =\frac{1}{2a}\left[\delta_{2j} -\delta_{1j}\right]
$$
$$
d_{1j} =\frac{1}{a^2}\left[\delta_{2j} - 2\delta_{1j}+\delta_{1j}\right]
$$
$i=N$
$$
c_{Nj} =\frac{1}{2a}\left[ -\delta_{N-1 j}\right]
$$
$$
d_{Nj} =\frac{1}{a^2}\left[ - 2\delta_{N j}+\delta_{N-1 j}\right]
$$
그 외에
$$
c_{ij} =\frac{1}{2a}\left[\delta_{i+1 j} -\delta_{i-1 j}\right]
$$
$$
d_{ij} =\frac{1}{a^2}\left[\delta_{i+1 j} - 2\delta_{i j}+\delta_{i-1 j}\right]
$$
.
하면, 만약, 만약...

function calc_cij(i,j,a,N)
    cij = 0.0
    if i==1
        cij = ifelse(j==2,1.0,0.0)        
        cij += ifelse(j==1,-1.0,0.0)
    elseif i==N
        cij = ifelse(j==N-1,-1.0,0.0)
    else
        cij = ifelse(j==i+1,1.0,0.0)
        cij += ifelse(j==i-1,-1.0,0.0)
    end
    cij = cij/(2a)
    return cij
end

function calc_dij(i,j,a,N)
    dij = 0.0
    if i==1
        dij += ifelse(j==2,1.0,0.0)
        dij += ifelse(j==1,-1.0,0.0) #-2+1
    elseif i==N
        dij += ifelse(j==N,-2.0,0.0)        
        dij += ifelse(j==N-1,1.0,0.0)
    else
        dij += ifelse(j==i+1,1.0,0.0)
        dij += ifelse(j==i,-2.0,0.0)
        dij += ifelse(j==i-1,1.0,0.0)
    end
    dij = dij/(a^2)
    return dij
end

.

행렬 생성

이렇게 되면 미분산자가 분화될 수 있다, 한미어턴 행렬

function make_Hr(a,N,n)
    mat_Hr = zeros(Float64,N,N)
    for i in 1:N
        r = a*i -a/2
        mat_Hr[i,i] = -1/r
        for dr in -1:1
            j = i + dr
            if 1 <= j <= N
                cij = calc_cij(i,j,a,N)
                dij = calc_dij(i,j,a,N)
                mat_Hr[i,j] += -(dij/2 + cij/r)
            end
        end
    end
    return mat_Hr
end

.
결과 보기 코드는 다음과 같습니다.

using LinearAlgebra 
using Plots
N = 1000
rc = 30
rcs = [30,60,120]
nmax = 10
for i = 1:length(rcs)
    rc = rcs[i]
    a = rc/N
    n = 0
    mat_Hr = make_Hr(a,N,n)

    ε,ψ = eigen(mat_Hr)
    ε = sort(ε)
    if i==1
        plot(ε[1:nmax],label="rc = $rc")
    else
        plot!(ε[1:nmax],label="rc = $rc")
    end

end
savefig("eigen.png")

네.

결실

얻은 결과는...
유한원법으로 수소 원자의 Schr≠dinger 방정식을 수치 구해(C++와 Juria의 소스 코드)
네.
선생님의 글

의 엄밀한 비교를 rc=120이면 좋겠다.