/ramiro/demec/ufpe

Método direto sequencial

Introdução

Vamos mostrar aqui a implementação mais convencional possível, que terá um desempenho ruim, com todos os laços em Python e o acesso individual a cada elemento dos vetores. Não há realmente muito o que fazer em relação a isto com esta implementação.

Vamos também comparar os tempos de execução desta implementação com uma implementação de um solver para matrizes em banda, disponível na bibliteca SciPy.

Implementação Clássica

A implementação é praticamente uma tradução direta do algoritmo para a linguagem Python.

Em primeiro lugar precisamos reorganizar o arquivo de suporte, eliminando as coisas relacionadas à iteração de Gauss-Seidel. Nada de muito excitante aqui.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
"""
tridisup01.py

This just organizes a few support functions for the Direct Poisson Solver.
"""

import logging
import numpy       as     np
from   math        import pi

__format = "%(asctime)s: %(message)s"
logging.basicConfig(format=__format, level=logging.INFO, datefmt="%H:%M:%S")

def info(*args, **kwargs):
    """Directly logs info calls, 99% of what we're doing."""
    logging.info(*args, **kwargs)

def zero_rhs(x, L):
    """Sets up a null vector for the RHS of the equation"""
    return np.zeros_like(x)

def sin_rhs(x, L, A):
    """Sets up a rhs that produces a single period of a senoidal function with
       amplitude A as the solution."""
    f = 2*pi/L
    return -f*f*A*np.sin(f*x)

def parab_rhs(x, L, A):
    """
    Untested rhs for parabolic solution.
    """
    return np.full_like(x, 2*A/(L*L))

def cub_rhs(x, L, A):
    """Sets up rhs to have a cubic solution with zeros on 0, L/2, L/3
       and value A on L.
       Don't forget to set up the correct boundary conditions."""
    return -9*(A*L - 3*A*x)/(L*L*L)

A classe que define o problema agora também inclui a montagem da matriz de coeficientes e sua triangularização com o algoritmo de Thomas.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
"""

"problem01.py"

Computational domain definitions.
"""
import numpy as np
from tridisup01 import zero_rhs

#%%
class TriDiag:
    """
    A tridiagonal matrix with diagonals stored as vectors o size N.
    Columns will be named 'a', 'b', 'c', as is traditional.
    """
    def __init__(self, size):
        self._a = np.full(size, 1.0)
        self._b = np.full(size, -2.0)
        self._c = np.full(size, 1.0)
        self._a[0] = 0.0
        self._c[-1] = 0.0

    def solve(self, u, f):
        """
        Solves the system of equations with the given rhs.
        Returns the solution.
        If "keep" is true, do not change the coeficient matrix and rhs.
        """
        a = self._a
        b = self._b
        c = self._c
        for i in np.arange(1, b.size):
            fac = a[i]/b[i-1]
            b[i] -= c[i-1]*fac 
            f[i] -= f[i-1]*fac
        u[-1] = f[-1]/b[-1]
        for i in np.arange(b.size-2, -1, -1):
            u[i] = (f[i] - c[i]*u[i+1])/b[i]

#%%
class Problem:
    def __init__(self, L, nnode, U0=None, UL=None, calc_rhs=zero_rhs, **kwargs):
        """
        L: float      -- Length of the domain
        nnode: int    -- Number of nodes
        U0, UL: float -- Dirichlet boundary conditions
        calc_rhs function(x, L, *args) -- computes the rhs of the equation
        """
        self.x = np.linspace(0, L, nnode)
        self.dx = self.x[1] - self.x[0]
        self.dx2 = self.dx*self.dx
        self.u = np.zeros_like(self.x)
        self.L = L
        self.nn = nnode
        if U0 is not None: self.u[0] = U0 
        if UL is not None: self.u[-1] = UL 
        self.f = calc_rhs(self.x, self.L, **kwargs)

    def solve(self):
        trid = TriDiag(self.x.size-2)
        self.f *= self.dx2
        self.f[1] -= self.u[0]
        self.f[-2] -= self.u[-1]
        trid.solve(self.u[1:-1], self.f[1:-1])

Como a solução do problema foi basicamente trazida para dentro da classe que define o problema, o programa de simulação fica muito simples.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
"""
Direct solve with tridiagonal matrix algorithm
"""
import matplotlib.pyplot as plt
from problem02 import Problem
import tridisup01 as td

td.info("Setting up problems.")
A = 4.0
L = 3.0

nn = 31

probs_rbv = { 'linear': Problem(L, nn, UL=A, calc_rhs=td.zero_rhs),
          'sinusoidal': Problem(L, nn, calc_rhs=td.sin_rhs, A=A),
           'parabolic': Problem(L, nn, calc_rhs=td.parab_rhs, A=A),
               'cubic': Problem(L, nn, UL=A, calc_rhs=td.cub_rhs, A=A)}

for name, p in probs_rbv.items():
    p.solve()
    td.info("Solving {} problem.".format(name))

fig = 0
for name, p in probs_rbv.items():
    plt.figure(fig)
    plt.title("Solution for problem {}".format(name))
    plt.plot(p.x, p.u, 'b')
    fig += 1

Método direto do SciPy

Não há