Programação numérica

Cálculo de uma série

Cálculo de uma série

Podemos exprimir a função exponencial no ponto \(x\) como uma soma infinita: \[ e^x = 1 + x + \frac{x^2}{2!} + \cdots + \frac{x^n}{n!} + \cdots\,, \] onde \(n!\) é o fatorial de \(n\).

Vamos calcular \(e^x\) aproximamente somando um número finito de termos da série acima.

Quanto mais termos somarmos, melhor será a aproximação.

Soma truncada

Seja \(S_n\) a soma dos termos até ordem \(n\) da série: \[ S_n = 1 + x + \frac{x^2}{2!} + \cdots + \frac{x^n}{n!} \]

Vamos definir uma função

double soma_termos(double x, int n);

que calcula e retorna a soma \(S_n\) para \(x\) e \(n\) dados.

Calcular os termos da série

\[\begin{aligned} S_n &= 1 + x + \frac{x^2}{2!} + \cdots + \frac{x^n}{n!} \\[2ex] &= \underbrace{\frac{x^0}{0!}}_{\text{termo 0}} + \underbrace{\frac{x^1}{1!}}_{\text{termo 1}} + \underbrace{\frac{x^2}{2!}}_{\text{termo 2}} + \cdots + \underbrace{\frac{x^n}{n!}}_{\text{termo}~n} \\[2ex] &= T_0 + T_1 + T_2 + \cdots + T_n \\[2ex] T_k &= \frac{x^k}{k!} \end{aligned} \]

Calcular os termos da série

O primeiro termo é: \[ T_0 = \frac{x^0}{0!} = 1 \]

Podemos obter o termo seguinte \(T_{k+1}\) a partir de \(T_k\) sem calcular potências ou fatoriais: \[ T_{k+1} = T_k \times \frac{x}{k+1}\]

Implementação

double soma_termos(double x, int n) {
   double soma = 0.0;   /* acumulador */
   double termo = 1.0;  /* termo inicial */

   for (int i = 0; i <= n; i++) {
       soma = soma + termo;
       termo = termo * x/(i+1);
       /* atualizar termo i -> termo i+1 */
       /* NB: double/int -> double */
   }
   return soma;
}

Testar a função

#include <stdio.h>
#include <math.h>  // biblioteca matemática

int main(void) {
  printf("S5 = %.15f\n", soma_termos(1.0,5));
  printf("S10= %.15f\n", soma_termos(1.0,10));
  printf("S20= %.15f\n", soma_termos(1.0,20));
  printf("e = %.15f\n", exp(1.0));
                         // função biblioteca
  return 0;
}

Execução

$ gcc -Wall -o serie serie.c -lm
$ ./serie 
S5   = 2.716666666666666
S10  = 2.718281801146385
S20  = 2.718281828459046
e    = 2.718281828459045


A opção -lm indica ao compilador para ligar com a biblioteca de funções matemáticas (exp, sin, cos, fabs, etc.).

NB: no GCC em Intel X86 podemos omitir esta opção.

Outros critérios de paragem

Segunda versão

double soma_eps(double x, double eps) {
   double soma = 0.0;   // soma acumulada
   double termo = 1.0;  // termo inicial

   for(int i = 0; fabs(termo) > eps; i++) {
       // fabs(...): o valor absoluto
       soma = soma + termo;
       termo =  termo * x/(i+1);
   }
   return soma;
}

Exemplo de uso

#include <stdio.h>
#include <math.h>

int main(void) {
  printf("S1=%.10f\n", soma_eps(1.0, 1E-3));
  printf("S2=%.10f\n", soma_eps(1.0, 1E-6));
  printf("e=%.10f\n", exp(1.0)); 
        /* função da biblioteca math */
  return 0;
}

Execução

$ gcc -Wall -o serie serie.c -lm
$ ./serie 
S1=2.7180555556
S2=2.7182815256
e =2.7182818285

Cálculo de potências

Cálculo de potências

Solução simples

\[ x^n = \underbrace{x \times x \times \cdots \times x}_{n~\text{fatores}} \]

double potencia(double x, int n) {
   double r = 1.0;  /* resultado */
   /* repetir n vezes;
      i: contador de iterações */
   for (int i=0; i<n; i++) {
     r = r*x;
   }
   return r;
}

Exemplo

Quantas multiplicações para calcular potencia(3.0, 4)?

  x = 3.0;
  r = 1.0;
  r = r*x; //  3.0
  r = r*x; //  9.0
  r = r*x; //  27.0
  r = r*x; //  81.0

Exemplo (2)

Podemos calcular \(3^4\) usando apenas duas multiplicações:

   x = 3.0;    // x é 3.0
   x = x*x;    // x é 9.0
   x = x*x;    // x é 81.0

Vamos generalizar esta solução para qualquer \(n\geq 0\).

Propriedades de potências

\[ \begin{aligned} 3^4 = 3^{2\cdot 2} &= \left(3^2\right)^2 \\ & = \left(3\cdot 3\right)^2 \end{aligned} \]

Mais geralmente, para expoentes pares:

\[ \begin{aligned} x^{2n} &= \left(x^2\right)^n \\ & = (x\cdot x)^n \end{aligned} \]

Propriedades de potências (2)

\[ \begin{aligned} 3^5 = 3^{1+4} &= 3 \cdot 3^4 \\ &= 3 \cdot 3^{2\cdot 2} \\ &= 3 \cdot \left(3^2\right)^2 \\ & = 3 \cdot (3\cdot 3)^2 \end{aligned} \]

Mais geralmente, para expoentes ímpares:

\[\begin{aligned} x^{2n+1} &= x\cdot x^{2n} \\ &= x\cdot (x \cdot x)^n \end{aligned} \]

Propriedades de potências (3)

Combinando os dois casos:

\[ \begin{aligned} x^{2n} &= (x\cdot x)^n \\ x^{2n+1} &= x\cdot (x \cdot x)^n \end{aligned} \]

Temos ainda:

\[ x^0 = 1 \]

Estas propriedades dão-nos um algoritmo recursivo para calcular \(x^n\).

Algoritmo recursivo

\[ \begin{aligned} P(x,\, 0) &= 1 \\ P(x,\, n) &= P(x\cdot x,\, \frac{n}{2}), & (\text{se}~ n>0~\text{e par} ) \\ P(x,\, n) &= x \cdot P(x\cdot x,\, \frac{n-1}{2}), & (\text{se}~ n>0~\text{e impar}) \end{aligned} \]

Implementação

double potencia(double x, int n) {
   if(n == 0)  
       return 1.0;   /* caso base */
   else if(n%2 == 0) 
       /* o expoente é par */
       return potencia(x*x, n/2);
   else 
        /* o expoente é ímpar */
       return x * potencia(x*x, (n-1)/2);
}

Análise

Quantas multiplicações para calcular \(x^n\)?

Análise

Quando \(n\) aumenta, \(\log_2 n\) é muito menor do que \(n\):

\(n\) \(\log_2 n\)
2 1
8 3
32 5
512 9
1024 10

Por exemplo: para \(n=1000\) efetuamos:

Otimização

Implementação otimizada

double potencia(double x, int n) {
  double r = 1.0;
  while(n > 0) {
     if (n%2 != 0)
        r = r*x;
     x = x*x;
     n = n/2;
  }
  return r;
}