Computação Numérica em C

Transforme matemática em código funcional. Implemente os algoritmos que você vai encontrar com mais frequência em computação científica e de engenharia.

Estatística (Média, Variância, Desvio Padrão)

Calcular estatísticas básicas de um conjunto de dados — de longe a operação mais comum em relatórios científicos e de engenharia.
Média: x̄ = (1/n) Σxi Variância: s² = (1/n) Σxi² - x̄² Desvio padrão: s = √s²
#include <stdio.h>
#include <math.h>

int main(void) {
    double data[] = {72.5, 85.3, 91.0, 68.2, 77.8, 83.1, 90.5, 65.9};
    int n = 8;
    double sum = 0, sum2 = 0;

    for (int i = 0; i < n; i++) {
        sum  += data[i];       // soma
        sum2 += data[i] * data[i]; // soma dos quadrados
    }

    double avg = sum / n;
    double var = sum2 / n - avg * avg; // variância = E[X^2] - (E[X])^2
    double sd  = sqrt(var);

    printf("N:       %d\n", n);
    printf("Média:   %.2f\n", avg);
    printf("Var:     %.2f\n", var);
    printf("Std:     %.2f\n", sd);
    return 0;
}
Ponto-chave: a variância pode ser calculada em uma única passagem como "média dos quadrados menos o quadrado da média". Não precisa de dois laços.

Método de Newton (Busca Numérica de Raízes)

Um método numérico para achar uma solução de f(x) = 0. Ele usa repetidamente o ponto em que a reta tangente cruza o eixo x como próxima aproximação.
xn+1 = xn - f(xn) / f'(xn)
Exemplo: ache a solução positiva de x² - 2 = 0 (ou seja, √2 ≈ 1,41421356...).
#include <stdio.h>
#include <math.h>

double f(double x)  { return x*x - 2; }       // f(x) = x^2 - 2
double df(double x) { return 2*x; }            // f'(x) = 2x

int main(void) {
    double x = 2.0;  // chute inicial
    double EPS = 1e-10;

    for (int i = 0; i < 100; i++) {
        double x_new = x - f(x) / df(x);
        printf("iter %2d: x = %.15f  err = %e\n", i+1, x_new, fabs(x_new - x));
        if (fabs(x_new - x) < EPS) break;
        x = x_new;
    }
    printf("sqrt(2) = %.15f\n", x);
    return 0;
}
iter 1: x = 1.500000000000000 err = 5.00e-01 iter 2: x = 1.416666666666667 err = 8.33e-02 iter 3: x = 1.414215686274510 err = 2.45e-03 iter 4: x = 1.414213562374690 err = 2.12e-06 iter 5: x = 1.414213562373095 err = 1.59e-12 sqrt(2) = 1.414213562373095
15 dígitos de precisão em apenas 5 iterações. O método de Newton converge quadraticamente, dobrando aproximadamente o número de dígitos corretos a cada passo.
Ressalva: um chute inicial ruim pode não convergir. Cuidado com pontos onde f'(x) = 0.

Método de Monte Carlo (Estimando Pi)

Uma técnica de simulação que usa amostragem aleatória para responder uma pergunta de forma probabilística. Jogue pontos aleatórios num quadrado 1×1 e estime pi pela fração que cai dentro do círculo unitário.
pi ≈ 4 × (pontos dentro do círculo / total de pontos)
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(void) {
    int N;
    printf("Tentativas: "); scanf("%d", &N);
    srand((unsigned)time(NULL));

    int inside = 0;
    for (int i = 0; i < N; i++) {
        double x = (double)rand() / RAND_MAX; // aleatório em [0,1]
        double y = (double)rand() / RAND_MAX;
        if (x*x + y*y <= 1.0) inside++;     // dentro do círculo?
    }

    double pi = 4.0 * inside / N;
    printf("pi ≈ %.6f  (N=%d)\n", pi, N);
    printf("real: 3.141593\n");
    return 0;
}
Precisão vs. número de tentativas:
Tentativas NEstimativa (ex.)Erro
1.0003,156~0,01
10.0003,1432~0,001
100.0003,14052~0,001
1.000.0003,141284~0,0003
10.000.0003,141552~0,00004
A convergência é lenta (o erro encolhe como 1/√N), mas o método se estende facilmente a problemas em alta dimensão e é simples de programar. É muito usado em simulações de física e em engenharia financeira.

Integração Numérica (Regra do Trapézio & de Simpson)

Quando uma função não pode ser integrada analiticamente, você ainda consegue aproximar a integral numericamente.

Regra do Trapézio

∫ab f(x)dx ≈ h × [ f(a)/2 + f(x₁) + f(x₂) + ... + f(b)/2 ] (h = (b-a)/n)
#include <stdio.h>
#include <math.h>

double f(double x) { return sin(x); } // integrando

// Regra do trapézio
double trapezoid(double a, double b, int n) {
    double h = (b - a) / n;
    double sum = (f(a) + f(b)) / 2.0;
    for (int i = 1; i < n; i++)
        sum += f(a + i * h);
    return h * sum;
}

// Regra de Simpson (maior precisão)
double simpson(double a, double b, int n) {
    double h = (b - a) / n;
    double sum = f(a) + f(b);
    for (int i = 1; i < n; i++)
        sum += (i % 2 == 0) ? 2*f(a+i*h) : 4*f(a+i*h);
    return h / 3.0 * sum;
}

int main(void) {
    printf("Valor real da integral de sin(x) de 0 a pi = 2.0\n\n");
    int ns[] = {10, 100, 1000};
    for (int i = 0; i < 3; i++) {
        printf("n=%4d  trap=%.10f  Simpson=%.10f\n",
            ns[i], trapezoid(0, M_PI, ns[i]), simpson(0, M_PI, ns[i]));
    }
    return 0;
}
Comparação de precisão: a regra de Simpson (precisão de quarta ordem) é ordens de magnitude mais precisa que a do trapézio.

Sistemas Lineares (Eliminação Gaussiana)

Resolva um sistema Ax = b usando eliminação direta seguida de substituição reversa.
2x + y - z = 8 -3x - y + 2z = -11 -2x + y + 2z = -3
#include <stdio.h>
#include <math.h>
#define N 3

void gauss(double a[N][N+1], double x[N]) {
    // eliminação direta
    for (int i = 0; i < N; i++) {
        for (int j = i+1; j < N; j++) {
            double ratio = a[j][i] / a[i][i];
            for (int k = i; k <= N; k++)
                a[j][k] -= ratio * a[i][k];
        }
    }
    // substituição reversa
    for (int i = N-1; i >= 0; i--) {
        x[i] = a[i][N];
        for (int j = i+1; j < N; j++)
            x[i] -= a[i][j] * x[j];
        x[i] /= a[i][i];
    }
}

int main(void) {
    double a[N][N+1] = {
        { 2,  1, -1,   8},
        {-3, -1,  2, -11},
        {-2,  1,  2,  -3}
    };
    double x[N];
    gauss(a, x);
    for (int i = 0; i < N; i++)
        printf("x[%d] = %.4f\n", i, x[i]);
    return 0;  // resposta: x=2, y=3, z=-1
}

Método da Bisseção (Delimitando Raízes)

Se f(a) e f(b) têm sinais opostos, existe uma raiz entre eles. Divida o intervalo pela metade repetidamente. É mais lento que o método de Newton, mas com convergência garantida.
#include <stdio.h>
#include <math.h>

double f(double x) { return x*x*x - x - 1; } // x^3 - x - 1 = 0

double bisection(double a, double b, double eps) {
    while (b - a > eps) {
        double mid = (a + b) / 2.0;
        if (f(a) * f(mid) <= 0)
            b = mid;  // raiz na metade esquerda
        else
            a = mid;  // raiz na metade direita
    }
    return (a + b) / 2.0;
}

int main(void) {
    double root = bisection(1.0, 2.0, 1e-10);
    printf("raiz de x^3 - x - 1 = 0: %.10f\n", root);
    printf("checagem: f(raiz) = %e\n", f(root));
    return 0;
}
Bisseção: lenta (convergência linear), mas segura — sem precisar se preocupar com chute inicial.
Método de Newton: rápido (convergência quadrática), mas exige a derivada e um bom chute inicial.

Equações Diferenciais Ordinárias (Euler method & Runge-Kutta)

Equações da forma dy/dx = f(x, y) podem ser resolvidas passo a passo a partir de um valor inicial. Movimento de projétil, modelos epidêmicos, circuitos elétricos, cinética química — todos têm esse formato.

Problema teste

dy/dx = -y
y(0) = 1
Exato: y(x) = e-x
Integramos até x=5 com três métodos e comparamos. O valor exato é e-5 ≈ 0,00673795.

① Euler method (1ª ordem)

O esquema mais simples — avança usando a inclinação no ponto atual.
yn+1 = yn + h · f(xn, yn)

② Euler melhorado / Heun (2ª ordem)

Primeiro prediz, depois corrige usando a inclinação no ponto predito.
k1 = f(xn, yn)
k2 = f(xn+h, yn+h·k1)
yn+1 = yn + h · (k1 + k2) / 2

③ Runge-Kutta 4 (RK4, 4ª ordem)

Avalia a inclinação quatro vezes (início, dois pontos médios e fim) e tira uma média ponderada. O solver de EDO de uso geral mais utilizado.
k1 = f(xn, yn)
k2 = f(xn + h/2, yn + h/2 · k1)
k3 = f(xn + h/2, yn + h/2 · k2)
k4 = f(xn + h, yn + h · k3)
yn+1 = yn + h/6 · (k1 + 2·k2 + 2·k3 + k4)

Programa de comparação

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

// dy/dx = f(x, y)
double f(double x, double y) {
    (void)x;
    return -y;
}

double euler(double x0, double y0, double x_end, double h) {
    double x = x0, y = y0;
    while (x < x_end) {
        y = y + h * f(x, y);
        x = x + h;
    }
    return y;
}

double heun(double x0, double y0, double x_end, double h) {
    double x = x0, y = y0;
    while (x < x_end) {
        double k1 = f(x,     y);
        double k2 = f(x + h, y + h * k1);
        y = y + h * (k1 + k2) / 2;
        x = x + h;
    }
    return y;
}

double rk4(double x0, double y0, double x_end, double h) {
    double x = x0, y = y0;
    while (x < x_end) {
        double k1 = f(x,         y);
        double k2 = f(x + h/2,   y + h/2 * k1);
        double k3 = f(x + h/2,   y + h/2 * k2);
        double k4 = f(x + h,     y + h * k3);
        y = y + h / 6 * (k1 + 2*k2 + 2*k3 + k4);
        x = x + h;
    }
    return y;
}

int main(void) {
    double exact = exp(-5.0);
    printf("exato y(5) = %.8f\n\n", exact);
    printf("   h    | Euler       | Heun        | RK4\n");
    printf("--------+-------------+-------------+-------------\n");
    double hs[] = {0.5, 0.1, 0.01};
    for (int i = 0; i < 3; i++) {
        double h = hs[i];
        printf(" %5.3f  | %.8f  | %.8f  | %.8f\n",
               h, euler(0, 1, 5, h),
                  heun(0, 1, 5, h),
                  rk4(0, 1, 5, h));
    }
    return 0;
}
$ gcc ode_compare.c -o ode_compare -lm $ ./ode_compare exato y(5) = 0.00673795 h | Euler | Heun | RK4 --------+-------------+-------------+------------- 0.500 | 0.00097656 | 0.01000214 | 0.00679947 0.100 | 0.00515378 | 0.00681497 | 0.00673797 0.010 | 0.00657202 | 0.00673874 | 0.00673795
Ordem de precisão: dividir h por 10 reduz o erro em cerca de 10× (Euler, 1ª ordem), 100× (Heun, 2ª) e 10000× (RK4, 4ª). RK4 já bate 6 dígitos com h=0,1, enquanto Euler mal chega a 2 dígitos mesmo com h=0,01.

Extensão para sistemas (oscilador de mola)

Movimento de projétil, Lotka-Volterra e o modelo SIR envolvem várias variáveis acopladas, e o RK4 se estende naturalmente. Aqui resolvemos um oscilador de mola m·x'' = -k·x reescrevendo como duas equações de primeira ordem:
dx/dt = v, dv/dt = -(k/m)·x
#include <stdio.h>

#define K 1.0
#define M 1.0

// lado direito: retorna dx/dt e dv/dt
void deriv(double x, double v, double *dx, double *dv) {
    *dx = v;
    *dv = -(K / M) * x;
}

void rk4_step(double *x, double *v, double h) {
    double k1x, k1v, k2x, k2v, k3x, k3v, k4x, k4v;
    deriv(*x,              *v,              &k1x, &k1v);
    deriv(*x + h/2*k1x, *v + h/2*k1v, &k2x, &k2v);
    deriv(*x + h/2*k2x, *v + h/2*k2v, &k3x, &k3v);
    deriv(*x + h*k3x,    *v + h*k3v,    &k4x, &k4v);
    *x += h/6 * (k1x + 2*k2x + 2*k3x + k4x);
    *v += h/6 * (k1v + 2*k2v + 2*k3v + k4v);
}

int main(void) {
    double x = 1.0, v = 0.0;   // amplitude 1, solto do repouso
    double h = 0.05;
    printf("t,x,v\n");
    for (double t = 0; t <= 20; t += h) {
        printf("%.2f,%.4f,%.4f\n", t, x, v);
        rk4_step(&x, &v, h);
    }
    return 0;
}
Verificação de energia: em teoria E = ½·k·x² + ½·m·v² deve ser constante. O Euler simples desvia e a energia cresce; o RK4 a mantém essencialmente constante mesmo em intervalos longos de tempo. É uma das razões pelas quais o RK4 é tão usado.
Veja também
Os exemplos de projétil, SIR e Lotka-Volterra na página Construa Simulações ficam muito mais precisos quando você troca a atualização de Euler pelo RK4 — geralmente com menos passos.

Problemas de Desafio

Coloque o que aprendeu para funcionar nos problemas a seguir.
C1. Raiz cúbica pelo método de Newton
Ache x tal que x³ = a. Aplique o método de Newton com f(x) = x³ - a e f'(x) = 3x². Verifique que a=27 dá x=3.
C2. Integral definida por Monte Carlo
Calcule ∫₀¹ e-x² dx com Monte Carlo. (Faz parte da integral gaussiana — não tem forma fechada.)
Dica: para x aleatório em [0,1], calcule f(x) = e-x². A média desses valores aproxima a integral.
C3. Distribuição de frequência & histograma
Gere 100 números aleatórios em [0,1], construa uma distribuição de frequência com bins de largura 0,1 e imprima um histograma em texto.
C4. Intervalo de confiança
Estime pi via Monte Carlo N vezes, depois calcule um intervalo de confiança de 95% (x̄ ± 2σ/√N).
C5. Crivo de primos (Crivo de Eratóstenes)
Enumere todos os primos ≤ N com o Crivo de Eratóstenes. Use um array e confirme a complexidade O(N log log N).
C6. Multiplicação de matrizes como função
Escreva uma função que calcula C = AB para matrizes A e B de n×n. Verifique o resultado em matrizes 3×3.
Compartilhe este artigo
Compartilhar no X

Quiz de Revisão

Teste seu entendimento desta aula!

Q1. Qual a principal razão de erros de ponto flutuante ocorrerem em cálculos com double?

O binário só consegue representar valores com um número finito de dígitos
A CPU não é rápida o suficiente
O compilador tem bugs

Valores como 0,1 são dízimas em binário, então são armazenados como aproximações e introduzem erro de arredondamento.

Q2. Qual é a ideia central do método de Newton?

Aproximar a solução iterativamente usando retas tangentes da função
Buscar o valor com uma árvore binária
Estimar a solução por tentativas aleatórias

O ponto em que a reta tangente à estimativa atual cruza o eixo x vira a próxima estimativa, convergindo rapidamente para a solução.

Q3. Qual algoritmo é comumente usado para somar elementos de array com alta precisão?

Soma de Kahan
Quicksort
Algoritmo de Dijkstra

A soma de Kahan é uma técnica clássica que rastreia os bits de baixa ordem perdidos na adição de ponto flutuante e os reincorpora como termo de compensação.