hripz

n=5
A=np.tril(np.random.rand(n,n))
x = np.ones(n)
b=np.dot(A,x)
descente(A,b)

#**Q2.** Écrire la solution $x$ du système $Ax=b$ lorsque $A$ est cette fois-ci triangulaire supérieure, en fonction des coefficients de $A$ et de $b$, en résolvant successivement les équations depuis la dernière jusqu'à la première (on dit qu'on résout le système $Ax=b$ par *remontée*).
def remontée(A,b):
    ''' 
    Solve a linear lower triangular system A x = b 

    PARAMETERS 
    ----------
    A : numpy array of size N x N
    b : numpy array of size N

    RETURNS
    -------

    x numpy array of size N, the exact solution of A x = b.

    '''
   
    n = np.size(b)
    x = np.zeros(np.shape(b))
    
    if np.max(abs(np.tril(A,-1))) < 10**(-14): # Si sa partie triang. sup stricte = 0, i.e. si A tr. inférieure
        x[n-1] = b[n-1]/A[n-1,n-1]                    # Dans ce cas on résout Ax=b en descendant le système 
        for i in range(n-2, -1, -1):
            x[i] = (b[i]-np.vdot(A[i,i+1:],x[i+1:]))/A[i,i]
        return x
    else:
        print("A n'est pas triangulaire sup")
C=np.triu(np.random.rand(n,n))
x = np.ones(n)
d=np.dot(C,x)
remontée(C,d)

#**Q3.** Modifier votre fonction ``descente`` en une fonction que vous appelerez ``remonte_descente`` qui permet la résolution du système $Ax=b$ lorsque $A$ est triangulaire inférieure ou triangulaire supérieure. Votre fonction devra tester si la matrice $A$ est triangulaire supérieure ou inférieure.
def remontée_descente(A,b):
    n = np.size(b)
    x = np.zeros(np.shape(b))
    
    if np.max(abs(np.triu(A,1))) < 10**(-14): # Si sa partie triang. sup stricte = 0, i.e. si A tr. inférieure
        x[0] = b[0]/A[0,0]                    # Dans ce cas on résout Ax=b en descendant le système 
        for i in range(n-1):
            x[i+1] = (b[i+1]-np.vdot(A[i+1,0:i+1],x[0:i+1]))/A[i+1,i+1]
        return x
    elif np.max(abs(np.tril(A,-1))) < 10**(-14): # Si sa partie triang. sup stricte = 0, i.e. si A tr. inférieure
        x[n-1] = b[n-1]/A[n-1,n-1]                    # Dans ce cas on résout Ax=b en descendant le système 
        for i in range(n-2, -1, -1):
            x[i] = (b[i]-np.vdot(A[i,i+1:],x[i+1:]))/A[i,i]
        return x

print(remontée_descente(A,b), remontée_descente(C,d))