分析解决方案与解决方案之间的差异

Question

因此，我正在尝试求解微分方程$ \ frac {d^2y} {dx^2} = -y（x）$受边界条件y（0）= 0和y（1）= 1的约束，分析解决方案是y（x）= sin（x）/sin（1）。

我正在使用三点模板来近似双衍生物。 通过这些方式获得的曲线至少应在边界上匹配，但是即使在边界处，我的解决方案也有很小的差异。 我正在附上代码，请告诉我有什么问题。

import numpy as np
import scipy.linalg as lg
from scipy.sparse.linalg import eigs
from scipy.sparse.linalg import inv
from scipy import sparse
import matplotlib.pyplot as plt

a = 0
b = 1
N = 1000
h = (b-a)/N
r = np.arange(a,b+h,h)

y_a = 0
y_b = 1

def lap_three(r):
    h = r[1]-r[0]
    n = len(r)
   
    
    M_d = -2*np.ones(n)
   #M_d = M_d + B_d
    O_d = np.ones(n-1)
    mat = sparse.diags([M_d,O_d,O_d],offsets=(0,+1,-1))
   #print(mat)
    return mat


def f(r):
    h = r[1]-r[0]
    n = len(r)
    return -1*np.ones(len(r))*(h**2)

def R_mat(f,r):
    r_d = f(r)
    R_mat = sparse.diags([r_d],offsets=[0])
   #print(R_mat)
    return R_mat

#def R_mat(r):
 #   M_d = -1*np.ones(len(r))
    
    
def make_mat(r):
    main = lap_three(r) - R_mat(f,r)
    return main

main = make_mat(r)
main_mat = main.toarray()
print(main_mat)
'''
eig_val , eig_vec = eigs(main, k = 20,which = 'SM')

#print(eig_val)
Val = eig_vec.T
plt.plot(r,Val[0])
'''
main_inv = inv(main)
inv_mat = main_inv.toarray()
#print(inv_mat)
#print(np.dot(main_mat,inv_mat))
n = len(r)
B_d = np.zeros(n)
B_d[0] = 0
B_d[-1] = 1
#print(B_d)
#from scipy.sparse.linalg import spsolve
A = np.abs(np.dot(inv_mat,B_d))
plt.plot(r[0:10],A[0:10],label='calculated solution')

real = np.sin(r)/np.sin(1)

plt.plot(r[0:10],real[0:10],label='analytic solution')
plt.legend()
#plt.plot(r,real)
#plt.plot(r,A)

'''diff = A-real
plt.plot(r,diff)'''

Answer 1

无法保证arange（a，b+h，h）的最后一点，它主要是b，但在某些情况下也可以是B+H。更好地使用

r,h = np.linspace(a,b,N+1,retstep=True)

线性系统由中间位置的方程组成r [1]，...，r [n-1]。这些是n-1方程，因此您的矩阵大小太大了。

您可以通过在m_d中包含H^2术语来保持矩阵构造的短路。

如果使用稀疏矩阵，也可以使用稀疏求解器a = spsolve（main，b_d）。

因此，组成系统的方程式是

A[k-1] + (-2+h^2)*A[k] + A[k+1] = 0

右侧的向量，因此需要包含值-a [0]和-a [n]。这应该清除标志问题，无需用绝对值作弊。

解决方案向量a从开始构建的，r [1：-1]对应。由于邮政没有值0和n内部的值，因此也没有区别。

---

PS：这里不涉及放松，因为这不是迭代方法。也许您的意思是有限的差异方法。