Discordance entre solution analytique et solution par méthode de relaxation

J'essaie donc de résoudre l'équation différentielle $\frac{d^2y}{dx^2} = -y(x)$ sous réserve des conditions aux limites y(0) = 0 et y(1) = 1, la solution analytique est y(x) = sin(x)/sin(1).

J'utilise un gabarit à trois points pour approximer la dérivée double. Les courbes obtenues de ces manières doivent correspondre au moins aux limites, mais mes solutions présentent de petites différences même aux limites. Je joins le code, s'il vous plaît dites-moi ce qui ne va pas.

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)'''

Solution du problème

Il n'y a aucune garantie de ce que sera le dernier point arange(a,b+h,h), ce sera principalement b, mais pourrait dans certains cas aussi être b+h. Meilleure utilisation

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

Le système linéaire se compose des équations pour les positions médianes r[1],...,r[N-1]. Ce sont des N-1équations, donc la taille de votre matrice est trop grande.

You could keep the matrix construction shorter by including the h^2 term already in M_d.

Si vous utilisez des matrices creuses, vous pouvez également utiliser le solveur creux A = spsolve(main, B_d).

Les équations qui composent le système sont

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

Le vecteur de droite doit donc contenir les valeurs -A[0]et -A[N]. Cela devrait résoudre le problème de signe, pas besoin de tricher avec la valeur absolue.

Le vecteur solution Acorrespond, tel que construit dès le départ, à r[1:-1]. Comme il n'y a pas de valeurs pour les positions 0et Nà l'intérieur, il ne peut pas non plus y avoir de différence.

PS : Il n'y a pas de relaxation impliquée ici, avant tout parce qu'il ne s'agit pas d'une méthode itérative. Peut-être que vous vouliez dire une méthode de différence finie.