Sistemas Lineares

Resolução de Sistema Triangular Superior

Algorítmo

Dados:

Matriz \(A_{ij}\) triangular superior, \(B_i\) vetor.

Saída:

Vetor de com as soluções

Obtem a dimensionalidade \(n\) de \(A\) e \(B\)
\(x_n=\frac{b_n}{a_{nn}},\qquad soma=0\),
Para \(k=n-1\) até 1 (passo -1):
1. \(soma=b_k\)
2. Para \(j=k+1\) até \(n\):
  1. \(soma=soma-a_{kj}x_j\)
3. \(x_k=\frac{soma}{a_{kk}}\)

Implementação Python

Salve no arquivo sislinlib.py


# Rotina para resolução de sistema Triangular Superior.
def ts(A,B):
    """
       Resolve o sistema linear AX=B, onde A é uma matriz triangular superior. 
       Entrada:
       A - matrix triangular superior
       B - vetor coluna

       Saida: 
       X - vetor coluna com a solução
    """
    

    # Copia os argumentos para poder manipula-los
    # Isso é necessário em Python quando os argumentos são vetores e matrizes.
    a=A.copy()
    b=B.copy()

    # Obtem dimensionalidade de B
    n=B.size

    # x será do mesmo tamanho que B mas inicializado com zeros
    x=B.copy()
    x[:]=0

    # Em python, vetores e matrizes tem índices começando a partir do 0.
    # ou seja, x(i) será x(0), x(1), ..., x(n-1)
    
    # A função "range(a,b)" gera uma sequência de números de a até b-1 !
    # Para se criar uma sequência 1,2,3,4 temos range(1,5)
    # A função "range(b,a,-1)" gera uma sequência de b descendo até a+1 !
    # Para se criar uma sequência 4,3,2,1 temos range(4,0,-1)

    x[n-1]=b[n-1]/a[n-1,n-1];
    for k in range(n-1,0,-1):
        soma=b[k-1]
        for j in range(k+1,n+1):
            soma=soma-a[k-1,j-1]*x[j-1]
        x[k-1]=soma/a[k-1,k-1]

    return x

Rotina Principal


import numpy as np

import sislinlib as sl

# Teste 
A=np.array([[2,-1,1],[0,1,2],[0,0,1]])
B=np.array([2,3,1])

x=sl.ts(A,B)

Eliminação Gaussiana

Este método transforma um sistema linear geral em um sistema triangular superior.

Algorítmo

Dados:

Matriz \(A_{ij}\).
Vetor \(B_i\).

Saída:

Matriz \(AT_{ij}\) triangular superior.
Vetor \(BT_i\) vetor de constantes do novo sistema.

Obtem a dimensionalidade \(n\) de \(A\) e \(B\)
Para \(k=1\) até \(n-1\):
1. \(w=|a_{kk}|\) (Aplica o pivoteamento - usa a linha com maior pivot)
2. Para \(j=k\) até \(n\):
  1. Se \(|a_{jk}|>w\):
  2. \(w=|a_{jk}|\) e \(r=j\)
3. Para \(j=1\) até \(n\): (Troca linhas \(k\) e \(r\))
  1. \(v_j=a_{kj}\)
  2. \(a_{kj}=a_{rj}\)
  3. \(a_{rj}=v_j\)
4. Para \(i=k+1\) até \(n\): (Calcula novos coeficientes A e B)
  1. \(m=\frac{a_{ik}}{a_{kk}}\)
  2. \(b_i=b_i-m_{ik}b_k\)
  3. Para \(j=k+1\) até \(n\):
  4. \(a_{ij}=a_{ij}-m a_{kj}\)

Nota: Porque o pivoteamento funciona. A equação matricial é da forma: \[ AX=B,\] ou em forma de somatório fica \[ \Sigma_j A_{ij}X_j = B_i,\] ou seja, para cada linha \(i\) da matrix \(A\) e do vetor \(B\), temos a somatória acima. Se trocarmos quem uma linha por outra, por exemplo, trocarmos os números da linha 1 com a linha 3, de ambos \(A\) e \(B\), a soma ainda vai ser a mesma para os mesmos valores de \(X\).

Implementação Python

Salve no arquivo sislinlib.py


# Eliminação Gaussiana.
def eg(A,B):
    """Usa o método da Eliminação Gaussiana para 
    transformar o sistema linear AX=B em um sistema Triangular superior,
    retornando a nova matriz At, triangular superior, e o vetor Bt, 
    tal que AtX=Bt tem a mesma solução."""

    # Necessário em Python quando trabalhando com matrizes e vetores
    a=A.copy()
    b=B.copy()

    # Obtem número de equacões a partir do tamanho do vetor b
    n=B.size
   
    # Loop da primeira até a penúltima linha (de 1 até n-1)
    for klin in range(1,n-1+1):
        #######################################
        # Aplica Pivoteamento
        w=abs(a[klin-1,klin-1])
        for jlin in range(klin,n+1):
           if abs(a[jlin-1,klin-1])>w:
              w=abs(a[jlin-1,klin-1])
              rlin=jlin
        # Terminado, troca as linhas "k" e "r"
        for jlin in range(1,n+1):
           vj=a[klin-1,jlin-1]
           a[klin-1,jlin-1]=a_[rlin-1,jlin-1]
           a[rlin-1,jlin-1]=vj
        #########################################
        # Calcula novos valores de A e B
        for ilin in range(klin+1,n+1):
            m=a[ilin-1,klin-1]/a[klin-1,klin-1]
            b[ilin-1]=b[ilin-1]-m*b[klin-1]
            for jlin in range(k+1,n+1):
               a[ilin-1,jlin-1]=a[ilin-1,jlin-1]-m*a[klin-1,jlin-1]

    return a,b

Rotina Principal


import numpy as np

import sislinlib as sl

a=np.array([[1,0,1],[1,1,0],[2,3,1]])
b=np.array([0,1,1]);

print('A=',a)
print('B=',b)

# Executa Eliminação Gaussiana:
a1,b1=sl.eg(a,b)

print('A1=',a)
print('B1=',b)

# Resolve o sistema triangular
x=sl.ts(a1,b1)

print('X=',x);

Zeros de fuções não-lineares

Gráficos


import matplotlib.pyplot at plt
import numpy as np

###############
# Define a funcao a ser estudada
def f(x):
    return x**3-2*x-5

################
# Faz uma figura da funcao
x = np.arange(-5,5,.1);
plt.plot(x,f(x));
plt.grid()
plt.ion()
plt.show()

Método da Bisseção

Algorítmo

Dados

A função \(f(x)\),
Os pontos \(a\) e \(b\) tais que \(f(a)f(b)<0\),
Uma tolerância máxima \(tol\)

Retorno:

Estimativa para a raiz: \(x^*\).

Para \(k=0\) até 50: (itera 50 vezes no maximo)
1. \(x=(a+b)/2\) (calcula ponto medio)
2. Se \(f(a)f(b)<0\), então: (escolhe lado onde o zero se encontra)
  1. \(b=x\)
3. Caso contrário:
  1. \(a=x\)
4. Se \(|b-a|<tol\), então: (verifica se chegou no limite)
  1. \(x^*=x\).
  2. Retorne \(x^*\).
Emite aviso: numero maximo de iteracoes atingidas
Retorna \(x^*=x\).

Implementação Python

Salvar em um arquivo chamado zeroslib.py


#####################
# Método da Bissetriz
#####################

def bissetriz(f, a, b, tol):
    """Aplica o método da Bissetriz para encontra o zero da função f(x).

       f  - função de uma variável real: y=f(x).
       a  - ponto 'x' mais a esquerda.
       b  - ponto 'x' mais a direita.
       tol- valor de tolerância.
    """

    for j in range(0,50):
        x=(a+b)/2
        if f(a)*f(x)<0:
            b=x
        else:
            a=x

        #print([a,b])

        if abs(b-a)<tol:
            return (b+a)/2

    print('Cheguei no final de 50 iterações...');
    return (b+a)/2

Rotina Principal

import matplotlib.pyplot as plt
import numpy as np
import zeroslib as zz


###############
# Define a funcao a ser estudada

def f(x):
    return x**3-2*x-5


################
# Faz uma figura da funcao
x=np.arange(-5,5,.1);
plt.plot(x,f(x))
plt.grid()
plt.ion()
plt.show()

################
print('Teste da Bissetriz')
zero = zz.bissetriz(f, 2, 3,1e-5)
print('x=',zero,'f(x)=',f(zero))

Método de Newton

Algorítmo

Dados:

Ponto inicial \(x_0\).
A função \(f(x)\) que se quer encontrar uma raiz.
A derivada da função \(f'(x)\).
Tolerância \(e\).
Variável de seleção entre dois (ou mais) métodos de parada \(imet\)

Retorno:

Extimativa \(x^*\) do zero da função \(f(x)\).

Para \(k=1\) até 50: (itera 50 vezes no maximo)
1. \(x_1=x_0-\frac{f(x_0)}{f'(x_0)}\) (calcula o novo x )
2. Se \(imet=1\) então:
  1. Se \(|f(x_1)|<e\), então \(x^*=x_1\). RETORNE.
3. Se \(imet=2\):
  1. Se \(|x_1-x_0|<e\), então \(x^*=x_1\). RETORNE.
4. \(x_0=x_1\). (novo x torna-se o inicial)
Emite aviso: numero maximo de iteracoes atingidas
Retorna \(x^*=x_1\).

Implementação Python

Salvar em no mesmo arquivo chamado zeroslib.py

####################################
# Método da Newton
def newton(f,df,x0,tol, imet):
    """Aplica o método de Newton para encontra o zero da função f(x).

       f  - função de uma variável real: y=f(x).
       df - derivada da função f(x): df/dx
       x0 - ponto 'x' inicial.
       tol- valor de tolerância.
       imet- Seleciona qual método de parada usar:
             imet=1 - para quando |f(x)|<e
             imet=2 - para quando |x(n)-x(n-1)|<e
    """
    for k in range(0,50):

        x1 = x0-f(x0)/df(x0)

        #print('x=',x0,'dx=',x1-x0)

        if imet==1:
            if abs(f(x1))<tol:
                return x1
        if imet==2:
            if abs(x1-x0)<tol:
                return x1
        x0 = x1
    print('Cheguei no final das 50 iterações...');
    return x0

Rotina Principal


import matplotlib.pyplot as plt
import numpy as np
import zeroslib as zz


###############
# Define a funcao a ser estudada

def f(x):
    return x**3-2*x-5

################
# Define a derivada de f para o metodo de newton
def df(x):
    return 3*x**2-2

################
# Metodo de Newton
################
# Teste no valor absoluto da função
print('Newton - teste |f|<e')
zero = zz.newton(f, df, 3, 1e-5, 1)
print('x=',zero,'f(x)=',f(zero))

################
# Test na tamanho do passo
print('Newton - teste |dx|<e')
zero = zz.newton(f, df, 3, 1e-5, 2)
print('x=',zero,'f(x)=',f(zero))

Método da Secante

Este método é equivalente ao método de Newton, quando voce substitui a derivada analítica por uma derivada numérica, sem tomar o limite:

\(f'(x_1)=\lim_{x_1->x_0}\frac{f(x_1)-f(x_0)}{x_1-x_0} \rightarrow f'(x_1)\approx \frac{f(x_1)-f(x_0)}{x_1-x_0}\).

Algorítmo

Dados:

Valor inicial \(x_0\) e um ponto próximo \(x_1\).
A função \(f(x)\).
Variável de seleção entre dois (ou mais) métodos de parada \(imet\)

Retorno:

Extimativa \(x^*\) do zero da função \(f(x)\).

Para \(k=1\) até 50: (itera 50 vezes no maximo)
1. \(f_0=f(x_0)\), \(f_1=f(x_1).\)
2. \(df = (f_1-f_0)/(x_1-x_0).\)
3. \(x_2=x_1-\frac{f(x_1)}{df}\) (calcula o novo x )
4. Se \(imet=1\) então:
  1. Se \(|f(x_1)|<e\), então \(x^*=x_1\). RETORNE.
5. Se \(imet=2\):
  1. Se \(|x_1-x_0|<e\), então \(x^*=x_1\). RETORNE.
6. \(x_0=x_1\). (novo x torna-se o inicial)
Emite aviso: numero maximo de iteracoes atingidas
Retorna \(x^*=x_1\).

Implementação Python

Salvar em no mesmo arquivo chamado zeroslib.py



####################################
# Método da Newton-Secante
def newtonsec(f,df,x0,x1,tol, imet):
    """Aplica o método de Newton para encontra o zero da função f(x).

       f  - função de uma variável real: y=f(x).
       df - derivada da função f(x): df/dx
       x0 - ponto 'x' inicial.
       x1 - ponto 'x' próximo a x0
       tol- valor de tolerância.
       imet- Seleciona qual método de parada usar:
             imet=1 - para quando |f(x)|<e
             imet=2 - para quando |x(n)-x(n-1)|<e
    """
    for k in range(0,50):

        f0 = f(x0)
        f1 = f(x1)
        df = (f1-f0)/(x1-x0)
        x2 = x1-f(x1)/df

        #print('x=',x0,'dx=',x1-x0)

        if imet==1:
            if abs(f(x2))<tol:
                return x2
        if imet==2:
            if abs(x2-x1)<tol:
                return x2
        x0 = x1
        x1 = x2
    print('Cheguei no final das 50 iterações...');
    return x0

Rotina Principal

import matplotlib.pyplot as plt
import numpy as np
import zeroslib as zz


###############
# Define a funcao a ser estudada

def f(x):
    return x**3-2*x-5 

################
# Define a derivada de f para o metodo de newton
def df(x):
    return 3*x**2-2




################
# Metodo Newton-Secante
################
# Test na tamanho do passo
print('Newtonsec - teste |dx|<e')
zero = zz.newtonsec(f, df, 3, 3.1, 1e-5, 2)
print('x=',zero,'f(x)=',f(zero))

Método de Newton para sistema de equações não-lineares

Implementação Python

Salve no arquivo zeroslib.py


##############
# Newton para sistemas não-lineares
##############


def newtonsys(F,J,x0,n,tol):

    # Cria vetor e matriz com zeros.
    f = np.zeros(n)
    jac = np.zeros((n,n))

    # Loop de iterações (max 50)
    for k in range(0,50):

        # Calcula valores de f e j
        for i in range(0,n):
            f[i]=F[i](x0)
            for j in range(0,n):
                jac[i,j]=J[i][j](x0)

        # Resolve Jv=-F
        # Ax=B <-> A=J, B=-F
        # Transforma B em triangular (Eliminação Gaussiana)
        A,B=eg(jac,-f)
        # Resolve o sistema trangular.
        v=ts(A,B)

        x1 = x0+v

        print(x1)
        if abs(v).max()<tol:
            return x1

        x0 = x1

    print('Cheguei no final das 50 iterações...');
    return x1

Rotina Principal


import zeroslib as zz


# Define a função f1
def f1(x):
    return -x[0]*(x[0]+1)+2*x[1]-18

# Define a função f2
def f2(x):
    return (x[0]-1)**2+(x[1]-6)**2 - 25

# Define a derivada  df1/dx1
def j11(x):
    return -2*x[0]-1

# Define a derivada  df1/dx2
def j12(x):
    return 2

# Define a derivada  df2/dx1
def j21(x):
    return 2*x[0]-2

# Define a derivada  df2/dx2
def j22(x):
    return 2*x[1]-12


# Junta as funções em um "vetor de funções"
F = [f1, f2]

# Junta as derivadas em uma "matriz de funções"
J = [ [j11, j12], [j21, j22]]


# Aplica o método de newton no sistema acima
zero = zz.newtonsys(F,J, (3,12), 2, 1e-6)