0

我的代码遇到了二阶耦合 PDE 问题的问题。为了澄清，我正在用有限差分解决这个问题，即使它不需要它，因为我有一个与这个类似的更普遍的问题，我想检查代码是否合理。

\begin{equation}
    -\partial^2_r u_{\uparrow}(r) -\frac{1}{r} \partial_r u_{\uparrow}(r) + u_{\uparrow}(r)
    + \partial_r u_{\downarrow}(r)
    +\frac{u_{\downarrow}(r)}{r} 
    - \Delta (r)  u_{\downarrow}(r) = 0
    \end{equation}

\begin{equation}
    -\partial_r u_{\uparrow}(r)
    + \Delta (r) u_{\uparrow}(r)
    -\partial^2_r u_{\downarrow}(r) -\frac{1}{r} \partial_r u_{\downarrow}(r) + \frac{u_{\downarrow}(r)}{r^2} - u_{\downarrow}(r) = 0
    \end{equation}

因为 up(r) 和 down(r) 是耦合的，它最初没有给我一个三对角矩阵，这使得很难施加边界条件 (up(0) =1, down(0) = up(40) = down (40) = 0) 在 RHS 中。我知道答案将是振幅减小的振荡，但我的解决方案会增加振幅。我试图重新排列行和列以使其有点对角线。它仍然给我同样的错误答案。

def sol_fin(ran):
            
    r = np.linspace(0, 40, ran)
    Δ = 0.1*np.heaviside(r-1,1)
    Δr = np.abs(r[2]-r[1])
    
    M = np.zeros([2*ran,2*ran])
    
    for i in range(1, ran-1):

        M[i-1,i] = -1/(Δr**2) + 1/(2*r[i]*Δr) #up_i_eq1
        M[i,i] = 2/(Δr**2)   + 1  #up_i+1_eq1
        M[i+1,i] = -1/(Δr**2) - 1/(2*r[i]*Δr) #up_i+1_eq1

        M[i-1,i+ran] =  0 - 1/(2*Δr)#down_i_eq1
        M[i,i+ran] = -Δ[i]  + 1/r[i] #down_i+1_eq1
        M[i+1,i+ran] = 1/(2*Δr) #10**-15 #down_i+2_eq1

        M[i+ran-1,i] =  0 + 1/(2*Δr)#up_i_eq2
        M[i+ran,i] =  Δ[i]  #up_i+1_eq2
        M[i+1+ran,i] = -1/(2*Δr) #up_i+2_eq2

        M[ran+i-1,i+ran] = -1/(Δr**2) + 1/(2*r[i]*Δr)#down_i_eq2
        M[ran+i,i+ran] = 2/(Δr**2)  +1/(r[i]**2) - 1 #down_i+1_eq2
        M[ran+i+1,i+ran] = -1/(Δr**2) - 1/(2*r[i]*Δr)#down_i+2_eq2

#16 ghost points
    #up1
    ##initials
    M[0,0] = M[1,1] #exact
    M[1,0] = -1/(Δr**2) - 1/(2*r[1]*Δr)  #approx
    
    ##finals
    M[ran-1,ran-1] =  2/(Δr**2)   + 1  #approx
    M[ran-2,ran-1] = -1/(Δr**2) + 1/(2*r[i-1]*Δr) #approx 
    
    #down1
    ##initials
    M[0,ran] = -Δ[1]  + 1/r[1] #approx
    M[1,ran] = M[2,ran+1] #exact 
    
    
    ##finals
    M[ran-2,2*ran-1] = M[ran-3,2*ran-2] #exact
    M[ran-1,2*ran-1] = -Δ[ran-1]  + 1/r[ran-1] #approx
    
    #up2
    ##initials
    M[ran,0] = 0.1 # exact
    M[ran+1,0] = M[ran+2,1] # exact
    
    ##finals
    M[2*ran-1,ran-1] = 0.1 # exact
    M[2*ran-2,ran-1] = M[2*ran-3,ran-2] # exact
   
    #down2
    ##initials
    M[ran,ran] = 2/(Δr**2)  +1/(r[1]**2) - 1 #approx
    M[ran+1,ran] =  -1/(Δr**2) - 1/(2*r[1]*Δr)#approx
    
    ##finals
    M[2*ran-1,2*ran-1] =  2/(Δr**2)  +1/(r[ran-1]**2) - 1 #approx
    M[2*ran-2,2*ran-1] = -1/(Δr**2) + 1/(2*r[ran-1]*Δr) #approx
    
#reorganizing lines and columns to make it somewhat tridiagonal    
    
    N = np.zeros([2*ran,2*ran])
    for i in range(0, ran): #change lines
        N[2*i,:] = M[i,:]
        N[2*i+1,:] = M[ran+i,:]
        N[2*ran-1,:] = M[2*ran-1,:]

    N1 = np.zeros([2*ran,2*ran])
    for i in range(0, ran): #change columns - adding regular order in pairs and the rest in the other
        N1[:,2*i] = N[:,i]
        N1[:,2*i+1] = N[:,ran+i]
        N1[:,2*ran-1] = N[:,2*ran-1]
    # b    
        b = np.zeros([2*ran])
        b[0] = 1 #I guess??
    
    # solution:
    solA = np.linalg.solve(N1, b)
    print(solA)
    #solA = solA / (np.linalg.norm(solA) + 10**-16) #normalization
    
    sol_ups = np.linspace(10**-6, 40, ran-1)
    sol_downs = np.linspace(10**-6, 40, ran-1)
        
    for i in range(0, ran-1):
            sol_ups[i] = solA[2*i]
        
    for i in range(0, ran-1):
            sol_downs[i] = solA[2*i+1]
        
    r_plot = np.linspace(10**-6, 40, ran-1)
    
    plt.plot(Jayx_plot, Jayy_plot_a,'--', label='up-j', color='r') 
    plt.plot(Jayx_plot, Jayy_plot_b,'--', label='down-j', color='b')
    plt.plot(r_plot, sol_ups, label='up', color='r'), plt.xlabel("r"),
    plt.plot(r_plot, sol_downs, label='down', color='b'), plt.legend(),
    
    plt.ylim(-0.7,0.7)
    plt.xlim(0, 42)
    plt.ylabel("Psi(r)"), plt.axhline(y=0,color='black')
    plt.title("Spinor components") #for r = {}".format(r_in)) 
    plt.show()
    
sol_fin(1000)

这是我的代码：

4

0 回答 0