Departamento de Física - Faculdade de Ciências e Tecnologia da Universidade de Coimbra

Física Computacional - Ficha 2 - Zeros de Funções

Rafael Isaque Santos - 2012144694 - Licenciatura em Física


In [1]:
import numpy as np
%matplotlib inline

Define-se a nossa função:

$f(x) = \sin (x)$

e a sua derivada:

$f'(x) = \cos (x)$


In [2]:
f = lambda x: np.sin(x)
df = lambda x: np.cos(x)
my_stop = 1.e-4
my_nitmax = 100000
my_cdif = 1.e-6

Método da Bisecção


In [3]:
def bi(a, b, fun, eps, nitmax):
    c = (a + b) / 2
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        if fun(a)*fun(c) < 0: b = c
        else: a = c
        c = (a + b) / 2
        it += 1
    return it, c, fun(c)

In [4]:
bi(2, 4, f, my_stop, my_nitmax)


Out[4]:
(11, 3.1416015625, -8.9089102066436886e-06)

Método da Falsa Posição:


In [5]:
def regfalsi(a, b, fun, eps, nitmax):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        if fun(a) * fun(c) < 0: b = c
        else: a = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c)

In [6]:
regfalsi(2, 4, f, my_stop, my_nitmax)


Out[6]:
(3, 3.1415903579556947, 2.2956340984908623e-06)

Método de Newton-Raphson:


In [7]:
def newtraph(c0, fun, dfun, eps, nitmax):
    c = c0
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        c = c - fun(c) / dfun(c)
        it += 1
    return it, c, fun(c)

In [8]:
newtraph(2, f, df, my_stop, my_nitmax)


Out[8]:
(6, 3.1415926536808043, -9.1011076112782062e-11)

Método da Secante


In [9]:
def secant(a, b, fun, eps, nitmax):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        a = b
        b = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c)

In [10]:
secant(2, 4, f, my_stop, my_nitmax)


Out[10]:
(3, 3.1415903579556947, 2.2956340984908623e-06)

In [11]:
def bi2(a, b, fun, eps, nitmax, cdif):
    c = (a + b) / 2
    c_prev = a
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        if fun(a)*fun(c) < 0: b = c
        else: a = c
        c_prev = c
        c = (a + b) / 2
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)

In [12]:
bi2(2, 4, f, my_stop, my_nitmax, my_cdif)


Out[12]:
(21,
 3.1415929794311523,
 -3.2584135910528161e-07,
 3.141592025756836,
 6.2783295730092144e-07,
 9.5367431640625e-07)

In [13]:
def regfalsi2(a, b, fun, eps, nitmax, cdif):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    c_prev = c + cdif/2
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        if fun(a) * fun(c) < 0: b = c
        else: a = c
        c_prev = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)

In [14]:
regfalsi2(2, 4, f, my_stop, my_nitmax, my_cdif)


Out[14]:
(5,
 3.1415926535897931,
 1.2246467991473532e-16,
 3.1415926536048882,
 -1.5094913867333564e-11,
 -1.5095036332013478e-11)

In [15]:
def newtraph2(c0, fun, dfun, eps, nitmax, cdif):
    c = c0
    c_prev = c + cdif/2
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        c_prev = c
        c = c - fun(c) / dfun(c)
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)

In [16]:
newtraph2(2, f, df, my_stop, my_nitmax, my_cdif)


Out[16]:
(7,
 3.1415926535897931,
 1.2246467991473532e-16,
 3.1415926536808043,
 -9.1011076112782062e-11,
 -9.1011198577461982e-11)

In [17]:
def secant2(a, b, fun, eps, nitmax, cdif):
    c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
    c_prev = c + cdif/2
    it = 1
    while not((np.abs(fun(c)) < eps and np.abs(c - c_prev) < cdif) or it > nitmax):
        a = b
        b = c
        c_prev = c
        c = (a * fun(b) - b * fun(a)) / (fun(b) - fun(a))
        it += 1
    return it, c, fun(c), c_prev, fun(c_prev), (c - c_prev)

In [18]:
secant2(2, 4, f, my_stop, my_nitmax, my_cdif)


Out[18]:
(5,
 3.1415926535897931,
 1.2246467991473532e-16,
 3.1415926536048882,
 -1.5094913867333564e-11,
 -1.5095036332013478e-11)

Exercício 2

$f(x) = x + e^{-k x^{2}}$


In [107]:
from scipy.misc import derivative


def newtraphd(c0, fun, eps, nitmax):
    c = c0
    dfun = lambda x: derivative(fun, x, 0.0001)
    it = 1
    while np.abs(fun(c)) > eps and it < nitmax:
        c = c - (fun(c) / dfun(c))
        it += 1
    return it, c, fun(c), dfun(c)

In [97]:
f2 = lambda x, k: x + np.e ** (-k * x**2) * np.cos(x)

In [105]:
f2_k1 = lambda x: f2(x, 1)
newtraph(0, f2_k1, df2_k1, 1e-4, my_nitmax)


Out[105]:
(15, -0.5883681212998364, 7.0172442686589065e-05)

In [122]:
f2_k50 = lambda x: f2(x, 50)
for i in range(1, 10+1): print(newtraphd(0, f2_k50, 1e-4, i))


(1, 0, 1.0, 0.99999999999988987)
(2, -1.0000000000001101, -1.0000000000001101, 0.99999999999988987)
(3, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(4, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)
(5, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(6, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)
(7, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(8, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)
(9, 1.1013412404281553e-13, 1.0000000000001101, 0.99999999998934275)
(10, -1.0000000000106573, -1.0000000000106573, 0.99999999999988987)

In [123]:
for i in range(1, 10+1): print(newtraphd(-0.1, f2_k50, 1e-4, i))


(1, -0.1, 0.50350053278289919, 7.0955553139445682)
(2, -0.17095998981128269, 0.057578346024847515, 4.9465471855514309)
(3, -0.18260009832839749, 0.003047567616852459, 4.424209558854808)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)
(4, -0.18328893722183393, 1.0527388237396851e-05, 4.3936547460581012)