Pular para conteúdo

Lista 05 · Newton e secante

Acompanha o capítulo 5. Métodos rápidos, que partem de um chute: a raiz cúbica da calculadora, um diodo num circuito, a órbita de um satélite. Confira sempre se \(f(x)\) ficou perto de zero.

Abrir no Colab

Baixar o caderno

Como esta lista funciona

Abra o caderno no Colab, escreva a sua solução no lugar do # TODO e rode a célula de teste logo abaixo. Ela mesma diz o que passou e o que falhou — você não precisa esperar correção de ninguém.

O código de cada exercício também está aqui embaixo, caso você prefira colar num editor local. Os testes comparam números com tolerância.

Exercícios marcados com ✏️ são à mão, no papel, com 4 casas decimais — como na parte em papel das avaliações.

Exercício 01

✏️ À mão. Queremos \(\sqrt[3]{10} = 2{,}154435\ldots\), a raiz de \(f(x) = x^3 - 10\).

a) Com \(f'(x) = 3x^2\) e o chute \(x_0 = 2\), faça 3 iterações de Newton (6 casas decimais), preenchendo \(x_i\), \(f(x_i)\), \(f'(x_i)\) e \(x_{i+1}\).

b) Faça 2 iterações da secante com \(x_{-1} = 2{,}0\) e \(x_0 = 2{,}2\).

c) Em cada método, quantos algarismos certos há na última iteração? Qual converge mais depressa?

Exercício 02

Escreva newton(f, df, x0, tol, max_iter), que aplica Newton-Raphson a partir de x0, com a derivada df dada, e para quando \(\varepsilon_a\) (%) ficar abaixo de tol ou depois de max_iter iterações.

def g(x):
    return x**3 - 10

def dg(x):
    return 3 * x**2

newton(g, dg, 2, 1e-8, 50)    # -> 2.154434690...

def newton(f, df, x0, tol, max_iter):
    """Raiz de f por Newton-Raphson a partir de x0; para quando eps_a (%) < tol."""
    x = x0
    # TODO: laço: guarda o anterior, x = x - f(x) / df(x), break quando eps_a < tol
    return x


# --- testes (não altere) ---
def g_ex02(x):
    return x**3 - 10


def dg_ex02(x):
    return 3 * x**2


r = newton(g_ex02, dg_ex02, 2, 1e-8, 50)
assert abs(r - 2.154434690031884) < 1e-9, f"veio {r}"
r = newton(g_ex02, dg_ex02, 2, 1e-8, 1)
assert abs(r - 13 / 6) < 1e-12, "com uma iteração só: 2 - (8 - 10) / 12"
print("Exercício 02: todos os testes passaram!")

Exercício 03

Escreva secante(f, x0, x1, tol, max_iter), que aplica o método da secante a partir dos chutes x0 e x1:

\[ x_{i+1} = x_i - \frac{f(x_i)\,(x_{i-1} - x_i)}{f(x_{i-1}) - f(x_i)}. \]

Cuidado com a atualização: depois de calcular o novo ponto, o "anterior" passa a ser o atual, e o atual passa a ser o novo — nessa ordem.

def secante(f, x0, x1, tol, max_iter):
    """Raiz de f pela secante a partir de x0 e x1; para quando eps_a (%) < tol."""
    x_ant = x0
    x = x1
    # TODO: laço: calcula o novo ponto, x_ant = x, x = novo, break quando eps_a < tol
    return x


# --- testes (não altere) ---
def g_ex03(x):
    return x**3 - 10


r = secante(g_ex03, 2.0, 2.2, 1e-8, 50)
assert abs(r - 2.154434690031884) < 1e-9, f"veio {r}"
r = secante(g_ex03, 2.0, 2.2, 1e-8, 1)
assert abs(r - 2.151057401812689) < 1e-9, "com uma iteração só"
print("Exercício 03: todos os testes passaram!")

Exercício 04

Como a calculadora faz. Escreva raiz_cubica(a), que calcula \(\sqrt[3]{a}\) (para \(a > 0\)) por Newton em \(f(x) = x^3 - a\), sem usar **(1/3). Comece com o chute \(x_0 = a\) e pare quando a mudança entre duas iterações for menor que \(10^{-12}\) vezes o valor atual (no máximo 100 iterações).

raiz_cubica(27)      # -> 3.0
raiz_cubica(1000)    # -> 10.0

Quantas iterações foram necessárias para \(a = 1000\)? E para \(a = 10^{9}\)? Por que o chute \(x_0 = a\) é ruim para números grandes?

def raiz_cubica(a):
    """Raiz cúbica de a > 0 por Newton, como uma calculadora."""
    x = a
    # TODO: laço (até 100): novo = x - (x**3 - a) / (3 * x**2); se a mudança
    # for menor que 1e-12 * novo, devolva novo; senão x = novo
    return x


# --- testes (não altere) ---
assert abs(raiz_cubica(27) - 3) < 1e-10
assert abs(raiz_cubica(1000) - 10) < 1e-10
assert abs(raiz_cubica(2) - 1.2599210498948732) < 1e-12
assert abs(raiz_cubica(1e9) - 1000) < 1e-8, "funciona também com números grandes"
print("Exercício 04: todos os testes passaram!")

Exercício 05

Eletrônica — o diodo. Uma fonte de \(V_s\) volts alimenta um resistor de \(R\) ohms em série com um diodo. A corrente no resistor, \((V_s - V)/R\), tem de ser igual à corrente no diodo, dada pela equação de Shockley, \(I_s\,(e^{V/V_T} - 1)\). A função corrente_residual(V, Vs, R) (já escrita) calcula a diferença entre as duas: a tensão \(V\) no diodo é a raiz dela.

Escreva tensao_no_diodo(Vs, R), que acha essa raiz por Newton com a derivada calculada pela diferença central (com h = 1e-6), a partir do chute \(V = 0{,}7\) V. Pare quando a mudança for menor que \(10^{-12}\) (no máximo 100 iterações).

tensao_no_diodo(5, 1000)    # -> ≈ 0.5742 V

Por que o chute 0,7 V é bom, e um chute de 5 V seria péssimo? (Calcule \(e^{5/0{,}0259}\).)

import numpy as np


def corrente_residual(V, Vs, R):
    """Corrente no resistor menos corrente no diodo, com tensão V no diodo. JÁ ESCRITA."""
    Is = 1e-12       # corrente de saturação do diodo (A)
    VT = 0.025852    # tensão térmica a 27 °C (V)
    return (Vs - V) / R - Is * (np.exp(V / VT) - 1)


def tensao_no_diodo(Vs, R):
    """Tensão (V) no diodo de um circuito fonte + resistor + diodo, por Newton
    com derivada numérica."""
    V = 0.7
    h = 1e-6
    # TODO: laço: f, df pela central, novo = V - f / df, pare se a mudança < 1e-12
    return V


# --- testes (não altere) ---
V = tensao_no_diodo(5, 1000)
assert abs(V - 0.5741915026813174) < 1e-9, f"veio {V}"
assert abs(tensao_no_diodo(12, 470) - 0.6181293397517915) < 1e-9
print("Exercício 05: todos os testes passaram!")

Exercício 06

Astronomia — a equação de Kepler. Escreva anomalia_excentrica(M, e), que resolve \(M = E - e\,\sin E\) para \(E\) por Newton, com \(f(E) = E - e\sin E - M\), \(f'(E) = 1 - e\cos E\) e o chute \(E_0 = M\). Pare quando a mudança for menor que \(10^{-12}\) (no máximo 50 iterações).

anomalia_excentrica(1.0, 0.0167)    # a órbita da Terra -> ≈ 1.01418
anomalia_excentrica(0.5, 0.74)      # um satélite Molniya -> ≈ 1.1858

import numpy as np


def anomalia_excentrica(M, e):
    """Solução E da equação de Kepler, M = E - e*sen(E), por Newton."""
    E = M
    # TODO: laço de Newton; pare quando a mudança for menor que 1e-12
    return E


# --- testes (não altere) ---
import numpy as np


E = anomalia_excentrica(1.0, 0.0167)
assert abs(E - 0.0167 * np.sin(E) - 1.0) < 1e-12, "E tem de satisfazer a equação de Kepler"
assert abs(anomalia_excentrica(0.5, 0.74) - 1.1858441224708631) < 1e-10
assert anomalia_excentrica(2.0, 0.0) == 2.0, "órbita circular: E = M"
print("Exercício 06: todos os testes passaram!")

Exercício 07

Newton sem derivar à mão. Escreva newton_numerico(f, x0, tol, max_iter), igual ao newton do ex02, mas sem receber a derivada: calcule-a em cada iteração pela diferença central com \(h = 10^{-6}\).

Teste na função do paraquedista do capítulo 4 (massa que leva a 36 m/s em 4 s), partindo de 100 kg. A resposta é 142,7376 kg.

def newton_numerico(f, x0, tol, max_iter):
    """Newton com a derivada pela diferença central (h = 1e-6)."""
    x = x0
    h = 1e-6
    # TODO: em cada iteração, a derivada pela central e o passo de Newton
    return x


# --- testes (não altere) ---
import numpy as np



def paraquedista_ex07(m):
    return np.sqrt(9.81 * m / 0.25) * np.tanh(np.sqrt(9.81 * 0.25 / m) * 4) - 36


def cubo_ex07(x):
    return x**3 - 10


assert abs(newton_numerico(paraquedista_ex07, 100, 1e-8, 50) - 142.7376331) < 1e-5
assert abs(newton_numerico(cubo_ex07, 2, 1e-8, 50) - 2.154434690031884) < 1e-8
print("Exercício 07: todos os testes passaram!")

Exercício 08

O melhor dos dois mundos. A bisseção é segura; Newton é rápido. Escreva newton_protegido(f, df, a, b, tol), que:

  1. começa no ponto médio de [a, b] (onde f troca de sinal);
  2. a cada iteração, calcula o passo de Newton; se ele cair fora de (a, b), usa o ponto médio no lugar;
  3. atualiza o intervalo com o teste de sinal, como na bisseção (a raiz continua presa);
  4. para quando \(|f(x)|\) ficar menor que tol (no máximo 200 iterações).

Teste em \(x^{10} - 1\) no intervalo \([0; 1{,}3]\) — onde Newton puro, partindo de \(0{,}5\), foi atirado para 51,65.

def newton_protegido(f, df, a, b, tol):
    """Newton que nunca sai de [a, b]: se o passo de Newton cair fora do
    intervalo, usa o ponto médio (um passo de bisseção)."""
    x = (a + b) / 2
    # TODO: laço (até 200): passo de Newton; se cair fora de (a, b), ponto médio;
    # pare se abs(f(novo)) < tol; senão, teste de sinal para atualizar a ou b
    return x


# --- testes (não altere) ---
def g_ex08(x):
    return x**10 - 1


def dg_ex08(x):
    return 10 * x**9


def h_ex08(x):
    return x**3 - 2 * x + 2


def dh_ex08(x):
    return 3 * x**2 - 2


r = newton_protegido(g_ex08, dg_ex08, 0, 1.3, 1e-12)
assert abs(r - 1) < 1e-10, f"a raiz de x**10 - 1 é 1 (veio {r})"
r = newton_protegido(h_ex08, dh_ex08, -3, 0, 1e-12)
assert abs(h_ex08(r)) < 1e-10, "x**3 - 2x + 2 tem raiz entre -3 e 0 (perto de -1.77)"
print("Exercício 08: todos os testes passaram!")