Praktikum Metode Numerik - Pertemuan 2

Kembali ke Metode Numerik

OUTLINE

SymPy

Dalam pembelajaran metode numerik, seringkali kita perlu membandingkan hasil aproksimasi kita dengan nilai yang sesungguhnya. Seringkali pula, sebenarnya nilai yang sesungguhnya itu dapat kita peroleh (karena kita masih dalam tahap belajar; penerapan metode numerik di dunia nyata adalah pada kasus di mana nilai eksak tidak dapat diperoleh).

Hasil perhitungan eksak (seperti perhitungan menggunakan aljabar biasa atau ilmu kalkulus) juga disebut hasil perhitungan analitik atau simbolik. Istilah “analitik” bisa dianggap antonim dari istilah “numerik”.

Di Python, ada module/package bernama SymPy (symbolic Python) yang dapat melakukan perhitungan simbolik, seperti menghitung turunan, yang misalnya digunakan di metode Newton.

(Fun fact: aplikasi/package di komputer yang dapat melakukan perhitungan simbolik disebut Computer Algebra System (CAS). Beberapa contoh CAS adalah SymPy, Wolfram Mathematica, dan Maple.)

Mari kita import sympy:

import sympy

Seperti untuk NumPy dan tabulate, apabila terjadi error karena sympy tidak ditemukan, artinya package sympy belum terinstall, dan bisa di-install menggunakan pip install sympy (atau dengan tanda seru: !pip install sympy)

pip install sympy
Requirement already satisfied: sympy in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (1.11.1)
Requirement already satisfied: mpmath>=0.19 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from sympy) (1.2.1)
Note: you may need to restart the kernel to use updated packages.

Tentunya, penggunaan SymPy melibatkan variabel. Misalnya, kita ingin melakukan perhitungan simbolik dengan variabel \(x\). Kita perlu memberitahu SymPy, dengan syntax seperti berikut:

x = sympy.symbols("x")

Artinya, kita baru saja memberitahu SymPy bahwa, pada string apapun yang dijumpai oleh SymPy, huruf “x” perlu dianggap sebagai simbol, atau lebih tepatnya sebagai variabel.

Perhatikan pula bahwa kode di atas adalah assignment ke variabel pemrograman yang juga bernama x. Dengan demikian, untuk ke depannya, variabel x yang kita ketik di mana saja pada program kita akan dianggap sebagai variabel “x” oleh SymPy.

Dengan variabel x tersebut, kita dapat mendefinisikan suatu expression (ekspresi atau kalimat matematika), misal \(5x^4\), seperti berikut:

polinom = 5 * (x ** 4) / 2
print(polinom)
5*x**4/2

SymPy memiliki fitur pprint (pretty print), yaitu menampilkan suatu ekspresi secara cantik atau indah, layaknya seperti kita tulis di kertas:

sympy.pprint(polinom)
   4
5⋅x 
────
 2

Untuk melakukan diferensiasi atau menghitung turunan (dalam hal ini secara simbolik/analitik), gunakan sympy.diff:

turunan = sympy.diff(polinom, x)
sympy.pprint(turunan)
    3
10⋅x

dengan begitu, SymPy menghitung turunan dari ekspresi polinom yang kita berikan itu, terhadap variabel x. Sebenarnya, mengetik sympy.diff(polinom) saja sudah cukup, tapi lebih lengkap lebih baik.

Sejauh ini, semua ekspresi yang kita jumpai masih berbentuk simbol/tulisan, sehingga kita belum bisa men-substitusi variabel x dengan sembarang nilai. Misalnya kita ingin menjadikan ekspresi di atas sebagai suatu fungsi func(x), di mana kita bisa memasukkan nilai x apapun dan mendapatkan hasil. Caranya adalah menggunakan sympy.lambdify:

func = sympy.lambdify(x, turunan)
print(func(5))
1250

Pada syntax lambdify di atas, kita perlu memberitahu SymPy terlebih dahulu, variabel apa yang digunakan pada ekspresi tersebut; barulah kita tuliskan ekspresinya. Dalam hal ini, kita mengetik sympy.lambdify(x, turunan) karena sedang menggunakan variabel x untuk ekspresi turunan yang ingin kita ubah menjadi fungsi yang bisa di-substitusi nilai x nya.

Fungsi hasil lambdify sudah bisa digunakan seperti fungsi lainnya pada Python. Bahkan, kita bisa mencampur penggunaan SymPy dengan NumPy (maupun package lainnya). Contohnya, setelah tadi memperoleh func(x) dari SymPy:

import numpy as np
arr = np.array([2, 3, 5, 10])
print(func(arr))
[   80   270  1250 10000]

Seperti NumPy, SymPy juga memiliki fungsi sin, cos, log, exp dll, sehingga kita bisa melakukan perhitungan analitik yang melibatkan fungsi-fungsi tersebut.

g = x**2 * sympy.cos(x) + sympy.exp(-5*x)
print("Fungsinya:")
sympy.pprint(g)

gp = sympy.diff(g, x)
print("Turunannya:")
sympy.pprint(gp)
Fungsinya:
 2           -5⋅x
x ⋅cos(x) + ℯ    
Turunannya:
   2                          -5⋅x
- x ⋅sin(x) + 2⋅x⋅cos(x) - 5⋅ℯ

Meskipun kita bisa saja melakukan, misalnya, from sympy import cos, hal tersebut tidak disarankan, apalagi ketika program kita juga menggunkaan NumPy dengan from numpy import cos atau bahkan from numpy import *. Alasannya, dengan begitu, program bisa menjadi membingungkan, karena tidak ada pembeda antara cos dari NumPy (numerik) dengan cos dari SymPy (analitik/simbolik).

Namun, kalau Anda berhati-hati dan hanya melakukan hal tersebut untuk salah satu package saja, silakan.

Menariknya, SymPy bisa jadi lebih unggul daripada NumPy untuk beberapa perhitungan yang melibatkan akurasi tinggi, terutama untuk perhitungan yang sebenarnya bersifat analitik. Misalnya, kita tahu bahwa \(\sin(\pi) = 0\). Menurut SymPy,

print("Menurut SymPy, sin(pi) = " + str(sympy.sin(sympy.pi)))
Menurut SymPy, sin(pi) = 0

karena SumPy menghitung nilai sin dari \(\pi\) secara analitik, yaitu tanpa perlu menghitung nilai \(\pi\) (karena nilainya sudah jelas nol berdasarkan sifat fungsi sin). Sedangkan, NumPy mengaproksimasi nilai \(\pi\) terlebih dahulu, barulah hasil aproksimasi tersebut yang masuk ke fungsi sin. Hasil perhitungan fungsi sin tersebut pun juga aproksimasi, sehingga didapatkan hasil seperti berikut, yaitu sangat kecil tetapi bukan nol:

print("Menurut NumPy, sin(pi) = " + str(np.sin(np.pi)))
Menurut NumPy, sin(pi) = 1.2246467991473532e-16

di mana “e-16” artinya “dikali 10 pangkat -16”.

Metode Bisection

Metode Bisection adalah salah satu metode yang dapat kita gunakan dalam masalah pencarian akar (root finding). Akar dari suatu persamaan didefinisikan sebagai nilai \(x\) yang memenuhi \(f(x) = 0\). Misalkan \(f\) adalah suatu fungsi kontinu terdefinisi di \([a,b]\), di mana \(f(a)\) dan \(f(b)\) berlawanan tanda (sehingga pasti ada akar pada interval tersebut, menurut Teorema Nilai Antara / Intermediate Value Theorem).

Inti sari dari metode Bisection adalah

  1. menebak bahwa akar suatu persamaan ada di dalam interval tertentu \([a, b]\);
  2. menelusuri nilai fungsi pada nilai tengah atau rata-rata dari interval tersebut;
  3. mempersempit interval dengan memanfaatkan hasil rata-rata tersebut; dan
  4. terus mencari nilai tengah dari interval yang baru, yang kemudian dipersempit lalu dicari nilai tengahnya, dan seterusnya hingga akar ditemukan, atau hingga ukuran interval sudah cukup kecil sehingga memuaskan (yaitu sudah lebih kecil dari toleransi).

Didefinisikan nilai tengah dari interval:

\(p=\frac{(a+b)}{2}\)

Akan dicari \(f(p)\) dengan syarat sebagai berikut:

  • jika \(f(p) = 0\), maka \(p\) adalah akar dari \(f\)
  • jika \(f(p)f(a) > 0\), maka \(sign(f(p)) = sign(f(a))\). Sehingga, kita dapat mempersempit interval dengan memilih batasan baru yaitu a = p dan b tidak berubah.
  • jika \(f(p)f(a) < 0\), maka \(sign(f(p)) \neq sign (f(a))\), atau \(sign(f(p)) = sign(f(b))\). Sehingga, kita dapat mempersempit interval dengan memilih batasan baru yaitu a tidak berubah dan b = p.

Metode Bisection memiliki order of convergence = 1, atau disebut memiliki kekonvergenan linier (linear convergence). Artinya, dalam proses menemukan akar persamaan (konvergen menuju jawabannya), metode Bisection tidak secepat beberapa metode lainnya yang memiliki order of convergence yang lebih tinggi.

from numpy import sin, cos, tan, log, exp, sqrt, pi

formula = input('Masukkan formula fungsi: ')

def f(x):
    return eval(formula)

def Bisection(lower, upper, tolerance):
    if f(lower)*f(upper)<0:
        p0=lower
        p=(lower+upper)/2
        if f(p)==0:
            return p
        elif f(p)*f(lower)>0:
            lower=p
        elif f(p)*f(lower)<0:
            upper=p
 
        abs_error=abs(p0-p)
        p0=p
 
        while abs_error > tolerance:
            p=(lower+upper)/2
            
            if f(p)==0:
                break
            elif f(p)*f(lower)>0:
                lower=p
            elif f(p)*f(lower)<0:
                upper=p
        
            abs_error=abs(p0-p)
            p0=p
 
        return p
 
    else:
        if f(lower)*f(upper)>0:
            return "Metode gagal mengaproksimasi akar. Silahkan ubah batas atas atau batas bawah"
        else:
            if f(lower)==0:
                return lower
            else: #f(upper)==0
                return upper

low_bound = eval(input("Masukkan batas bawah interval: "))
up_bound = eval(input("Masukkan batas atas interval: "))
tolerance = eval(input("Masukkan toleransi aproksimasi: "))

akar_bisection=Bisection(low_bound,up_bound,tolerance)

try:
    print("Akar persamaan {0} = 0 adalah x = {1:.7f}".format(formula,
    akar_bisection))
except ValueError:
    print(akar_bisection)
Masukkan formula fungsi: 2*x - 3*cos(x) + exp(-5*x) - 9
Masukkan batas bawah interval: -3
Masukkan batas atas interval: 2
Masukkan toleransi aproksimasi: 10**(-7)
Akar persamaan 2*x - 3*cos(x) + exp(-5*x) - 9 = 0 adalah x = -0.5073225

Metode Fixed-Point

Inti sari dari Metode Fixed-Point adalah mencari fixed-point (titik tetap) dari suatu fungsi (misal fungsi \(g(x)\)), yaitu suatu nilai \(p\) sehingga \(p = g(p)\), atau \(p - g(p) = 0\). Titik \(p\) disebut titik tetap, karena ketika nilai \(p\) dimasukkan ke fungsi \(g(x)\), hasilnya tetaplah \(p\). Untuk nilai \(x\) yang dekat dengan \(p\), biasanya ada kecenderungan nilai \(g(x)\) menjadi semakin mendekati \(p\).

Perhatikan bahwa, sembarang persamaan \(f(x) = 0\) bisa diubah bentuknya dengan mendefinisikan fungsi \(g(x) = x - f(x)\) (sehingga \(f(x) = x - g(x)\)). Dengan demikian, permasalahan mencari akar berubah menjadi permasalahan mencari fixed-point, yaitu mencari nilai \(p\) sehingga \(p = g(p)\) atau \(p - g(p) = 0\) (sehingga nilai \(p\) tersebut juga menyebabkan \(f(p) = 0\)).

(Tentu saja, itu bukanlah satu-satunya cara untuk mengubah permasalahan mencari akar menjadi permasalahan mencari fixed-point. Bahkan, tidak semua pilihan \(g(x)\) yang memungkinkan itu dijamin memiliki fixed-point.)

Misalkan \(g\) adalah fungsi kontinu dan memiliki fixed-point \(p\) pada interval \([a,b]\) (dan diasumsikan bahwa \(g\) memenuhi persyaratan untuk kekonvergenan metode fixed-point). Artinya, ada \(p \in [a,b]\) sehingga \(g(x) = x\). Untuk mengaproksimasi penyelesaian dari persamaan \(g(x) = x\), diperlukan suatu tebakan awal \(p_0\), kemudian iterasinya adalah:

\(p_n = g(p_{n-1})\)

Nilai tersebut terus dimasukkan ke dalam \(g\) sehingga, diharapkan, nilai \(p_n\) menjadi semakin mendekati suatu nilai \(p\) yang membuat \(g(p) = p\).

Pada umumnya, metode fixed-point memiliki kekonvergenan linier. Ketika \(g(x)\) dijamin memliki tepat satu fixed-point (atau fixed-point yang unik) pada suatu interval \([a,b]\), maka Metode Fixed-Point dengan \(p_0\) pada interval tersebut pasti memiliki kekonvergenan linier. Terkadang Metode Fixed-Point lebih cepat daripada Metode Bisection, dan terkadang Metode Bisection lebih cepat daripada Metode Fixed-Point.

Hati-hati, ada kemungkinan bahwa \(g(p_n)\) malah menjauhi \(p\), contohnya untuk \(g(x) = x^2\) dan \(p_0 > 1\) (padahal \(g(1) = 1\)). Pada kasus seperti itu, metode fixed-point tidak dijamin konvergen (artinya tidak dijamin bisa menemukan fixed-point).

Sebagai contoh penggunaan metode fixed-point, kalian bisa mencoba untuk menyelesaikan persamaan (masalah mencari akar) berikut ini,

\[f(x) = x^2 - x - 1 = 0\]

dengan sedikit manipulasi aljabar (dibagi \(x\), pindah ruas) agar mendapatkan bentuk \(x = g(x)\),

\[x = 1 + \frac{1}{x}\]

sehingga, dengan \(g(x) = 1 + \frac{1}{x}\) bisa digunakan metode fixed-point, misal dengan tebakan awal \(x = 2\) atau \(x = -3\).

(Jelas metode ini akan gagal untuk \(g(x)\) tersebut apabila dipilih tebakan awal seperti \(x=0\), \(x=-1\), atau bahkan \(x=-\frac{1}{2}\) karena akan terjadi pembagian nol. Kemungkinan terjadinya pembagian nol itu bukan hanya dari metodenya seperti metode Newton, tetapi juga dari fungsi \(f(x)\) atau \(g(x)\) yang digunakan.)

Silakan coba dengan kode di bawah ini!

Sebagai pembanding, kalian bisa menyelesaikan persamaan kuadrat \(f(x) = x^2 - x - 1 = 0\) di atas, dan mendapatkan solusi

\[x_1 = \frac{1+\sqrt{5}}{2} \approx 1.618\]

\[x_2 = \frac{1-\sqrt{5}}{2} \approx -0.618\]

Kebetulan, konstanta berikut ini yang berlambang phi kecil (\(\phi\)),

\[\phi = \frac{1+\sqrt{5}}{2}\]

adalah konstanta istimewa yang bernama golden ratio.

from numpy import cos, sin, tan, log, exp, sqrt
from tabulate import tabulate

formula = input("Masukkan Formula Fungsi: ")

def g(x):
    return eval(formula)

def FixedPoint(p0,tolerance):
    table = [["iterasi","Aproksimasi"]]
    iterasi = []
    
    i = 1
    p = g(p0)
    abs_error = abs(p-p0)
    p0 = p
    iterasi.append(i)
    iterasi.append(p)
    table.append(iterasi)

    while abs_error > tolerance:
        iterasi = []
        i += 1
        p = g(p0)
        abs_error = abs(p-p0)
        p0 = p
        iterasi.append(i)
        iterasi.append(p)
        table.append(iterasi)
    
    print(tabulate(table,headers='firstrow',tablefmt ="pretty"))
    return p0
    
    print(tabulate(table,headers = 'firstrow',tablefmt="pretty"))
    return p0 # nilai terbaru, karena sudah dipasang p0 = p

starting_point = eval(input("Masukkan titik awal iterasi: "))
tolerance = eval(input("Masukkan batas toleransi: "))

fixed_point = FixedPoint(starting_point,tolerance)

print('Ditemukan fixed point dari g(x) = {0} yaitu x = {1:.7f}'.format(formula,fixed_point))
Masukkan Formula Fungsi: 1 + 1/x
Masukkan titik awal iterasi: 2
Masukkan batas toleransi: 10**(-7)
+---------+--------------------+
| iterasi |    Aproksimasi     |
+---------+--------------------+
|    1    |        1.5         |
|    2    | 1.6666666666666665 |
|    3    |        1.6         |
|    4    |       1.625        |
|    5    | 1.6153846153846154 |
|    6    | 1.619047619047619  |
|    7    | 1.6176470588235294 |
|    8    | 1.6181818181818182 |
|    9    | 1.6179775280898876 |
|   10    | 1.6180555555555556 |
|   11    | 1.6180257510729614 |
|   12    | 1.6180371352785146 |
|   13    | 1.6180327868852458 |
|   14    | 1.618034447821682  |
|   15    | 1.618033813400125  |
|   16    | 1.6180340557275543 |
|   17    | 1.6180339631667064 |
+---------+--------------------+
Ditemukan fixed point dari g(x) = 1 + 1/x yaitu x = 1.6180340

Metode Newton biasa (dengan turunan analitik)

Misalkan \(f\) kontinu dan terturunkan (memiliki turunan) di \([a,b]\) dan ada tebakan awal \(p_0 \in\) \([a,b]\) sedemikian sehingga \(f'(p_0) \neq 0\). Iterasi pada metode Newton untuk menyelesaian \(f(x) = 0\) adalah sebagai berikut:

\(p_n = p_{n-1} - \frac{f(p_{n-1})}{f'(p_{n-1})}\)

Diharapkan bahwa, setelah banyak iterasi, nilai \(p_n\) yang diperoleh akan membuat \(f(p) = 0\) atau setidaknya sangat dekat dengan nol (lebih kecil dari batas toleransi yang kita anggap sudah memuaskan).

Metode Newton juga dapat dipandang sebagai metode fixed-point dengan \(g(x) = x - \frac{f(x)}{f'(x)}\)

Metode Newton gagal apabila, pada suatu iterasi, tiba-tiba \(f'(p_n) = 0\).

Pada umumnya, Metode Newton memiliki order of convergence = 2, atau juga disebut memiliki kekonvergenan kuadratik (quadratic convergence). Artinya, selama berhasil, Metode Newton lebih cepat daripada Metode Bisection maupun Metode Fixed-Point.

import sympy
from numpy import sin, cos, tan, log, exp, sqrt

formula = input("Masukkan fungsi: ")
def f(x):
    return eval(formula)

x = sympy.symbols("x")

df_string = str(sympy.diff(formula, x))
def df(x): # turunan f
    return eval(df_string)

def Newton(p0,tolerance):
    p = p0 - f(p0)/df(p0)
    abs_error = abs(p-p0)
    p0 = p

    while abs_error > tolerance:

        try:
            p = p0 - f(p0)/df(p0)
        except ZeroDivisionError:
            return "Metode gagal mengaproksimasi akar. Silakan pilih tebakan awal lain"
        
        abs_error = abs(p-p0)
        p0 = p
    return p

starting_point = eval(input("Masukkan tebakan awal / titik awal iterasi: "))
tolerance = eval(input("Masukkan toleransi aproksimasi: "))

akar_newton = Newton(starting_point, tolerance)

print(f"Akar dari persamaan f(x) = {formula} adalah x = {akar_newton:.7f}")
Masukkan fungsi: 2*x - 3*cos(x) + exp(-5*x) - 9
Masukkan tebakan awal / titik awal iterasi: -1
Masukkan toleransi aproksimasi: 10**(-7)
Akar dari persamaan f(x) = 2*x - 3*cos(x) + exp(-5*x) - 9 adalah x = -0.5073225

Metode Newton dengan Beda Hingga (Finite-Difference Newton’s Method)

Salah satu kekurangan Metode Newton yang biasa adalah harus mengetahui rumus turunannya secara analitik. Sebelum adanya CAS seperti SymPy, turunan analitik harus dihitung secara manual dengan kalkulus. Kalau bentuk rumus untuk \(f(x)\) sangat rumit, perhitungan turunan menjadi jauh lebih rumit. Untuk menghindari menghitung turunan secara analitik, kita dapat menggunakan definisi turunan (yang menggunakan limit):

\[f'(x) = \lim_{h \rightarrow 0} \frac{f(x+h) - f(x)}{h}\]

dengan memilih nilai \(h\) yang cukup kecil (sayangnya, kita tidak bisa membuat limit \(h\) menuju nol). Nilai \(h\) yang cukup kecil itu disebut suatu beda hingga (finite difference).

Sehingga, modifikasi metode Newton ini bisa disebut Metode Newton dengan Beda Hingga (Finite-Difference Newton’s Method). Untuk fungsi \(f\) yang kontinu, akar persamaan \(f(x) = 0\) bisa ditentukan dengan iterasi sebagai berikut:

\(p_n = p_{n-1} - \frac{f(p_{n-1})}{\left(\frac{f\left(p_{n-1}+h\right)-f(p_{n-1})}{h}\right)} = p_{n-1} - \frac{f(p_{n-1})h}{f(p_{n-1}+h)-f(p_{n-1})}\)

dengan tebakan awal \(p_0\). Perhatikan bahwa \(f'(p_{n-1})\) pada metode Newton yang biasa itu telah digantikan dengan

\[f'(p_{n-1}) \approx \frac{f(p_{n-1}+h) - f(p_{n-1})}{h}\]

Tujuan modifikasi tersebut adalah agar iterasi dapat dilakukan pada titik di mana turunannya tidak ada, atau ketika turunan analitik sulit diperoleh.

from numpy import sin, cos, tan, log, exp, sqrt

formula = input("Masukkan fungsi: ")
def f(x):
    return eval(formula)

def df(x, h=10**(-12)):
    return (f(x+h)-f(x))/h

def Finite_Difference(p0,tolerance):
    p = p0 - f(p0)/df(p0)
    abs_error = abs(p-p0)
    p0 = p

    while abs_error > tolerance:
        p = p0 - f(p0)/df(p0)
        abs_error = abs(p-p0)
        p0 = p
    return p

starting_point = eval(input("Masukkan titik awal iterasi: "))
tolerance = eval(input("Masukkan toleransi aproksimasi: "))

akar_fd = Finite_Difference(starting_point,tolerance)

print(f"Akar dari persamaan f(x) = {formula} adalah x = {akar_fd:.7f}")
Masukkan fungsi: 2*x - 3*cos(x) + exp(-5*x) - 9
Masukkan titik awal iterasi: -1
Masukkan toleransi aproksimasi: 10**(-7)
Akar dari persamaan f(x) = 2*x - 3*cos(x) + exp(-5*x) - 9 adalah x = -0.5073225

Metode Secant

Pada Metode Newton dengan Beda Hingga, nilai \(h\) konstan. Kalau kita punya dua tebakan awal yang saling dekat, misal \(p_0\) dan \(p_1\), kita bisa saja memanfaatkannya dengan memasang \(h = p_1 - p_0\). Bahkan, ketika iterasi \(p_n\) sudah semakin dekat menuju akar, jarak antara \(p_{n-1}\) dan \(p_{n-2}\) menjadi semakin kecil. Sehingga, dengan memasang nilai \(h = p_{n-2} - p_{n-1}\) atau \(h = p_{n-1} - p_{n-2}\), kita berhasil membuat limit \(h\) menuju nol.

Modifikasi ini disebut Metode Secant, dengan iterasi sebagai berikut untuk menentukan penyelesaian \(f(x) = 0\) dengan fungsi \(f\) yang kontinu:

\(p_n = p_{n-1} - \frac{f(p_{n-1})}{\left(\frac{f(p_{n-1})-f(p_{n-2})}{p_{n-1}-p_{n-2}}\right)} = p_{n-1} - \frac{f(p_{n-1})(p_{n-1} - p_{n-2})}{f(p_{n-1}) - f(p_{n-2})}\)

Dibandingkan Metode Newton yang biasa, Metode Secant menggantikan \(f'(p_{n-1})\) dengan

\[f'(p_{n-1}) \approx \frac{f(p_{n-1}) - f(p_{n-2})}{p_{n-1} - p_{n-2}}\]

sehingga, tidak seperti Metode Newton yang hanya memerlukan satu tebakan awal, Metode Secant membutuhkan dua tebakan awal, yaitu \(p_0\) dan \(p_1\). Namun, dibandingkan dengan Metode Newton dengan Beda Hingga, nilai \(h\) atau beda hingga tersebut tidak perlu ditentukan secara manual.

Menariknya, Metode Secant memiliki order of convergence = \(\phi \approx 1.618\).

from numpy import sin, cos, tan, log, exp, sqrt
formula = input("Masukkan formula fungsi: ")

def f(x):
    return eval(formula)

def Secant(p0,p1,tolerance):
    p = p1 - (f(p1)*(p1-p0))/(f(p1)-f(p0))
    abs_error = abs(p-p1)
    p0 = p1
    p1 = p

    while abs_error > tolerance:
        p = p1 - (f(p1)*(p1-p0))/(f(p1)-f(p0))
        abs_error = abs(p-p1)
        p0 = p1
        p1 = p
    return p

starting_point1 = eval(input("Masukkan titik awal pertama: "))
starting_point2 = eval(input("Masukkan titik awal kedua: "))
tolerance = eval(input("Masukkan toleransi aproksimasi: "))

akar_secant = Secant(starting_point1,starting_point2,tolerance)

print(f"Akar dari persamaan f(x) = {formula} adalah x = {akar_secant:.7f}")
Masukkan formula fungsi: 2*x - 3*cos(x) + exp(-5*x) - 9
Masukkan titik awal pertama: -1
Masukkan titik awal kedua: -2
Masukkan toleransi aproksimasi: 10**(-7)
Akar dari persamaan f(x) = 2*x - 3*cos(x) + exp(-5*x) - 9 adalah x = -0.5073225

Metode Regula Falsi (penjelasan tanpa kode)

Sejauh ini, kita sudah membahas beberapa metode root-finding atau aproksimasi akar, yaitu:

  • Metode Bisection
  • Metode Fixed-Point
  • Metode Newton biasa (dengan turunan analitik)
  • Metode Newton dengan Beda Hingga (finite-difference Newton’s method)
  • Metode Secant

Di antara semua metode tersebut, hanya Metode Bisection yang dijamin konvergen menuju akar di interval yang diberikan; semua metode lain ada kemungkinan divergen (menjauh dari akar, seperti metode fixed-point) atau gagal karena terjadi pembagian nol. Sayangnya, Metode Bisection termasuk metode yang pelan di antara metode numerik lainnya.

Untuk menjaga jaminan kekonvergenan oleh Metode Bisection tetapi memperbaiki kecepatan kekonvergenannya, kita bisa memodifikasi Metode Bisection, yaitu memodifikasi cara menentukan \(p\) yang baru yang akan mempersempit interval. Perhatikan bahwa Metode Bisection membutuhkan dua “tebakan awal” (lebih tepatnya dua batasan interval), sedangkan metode di atas yang juga membutuhkan dua tebakan awal hanyalah Metode Secant.

Apakah kita bisa menggunakan Metode Bisection, tetapi dengan modifikasi menentukan \(p\) seperti Metode Secant, agar mendapatkan order of convergence seperti Metode Secant?

Jawabannya adalah bisa, dan modifikasi tersebut dinamakan Metode Regula Falsi. Sehingga, Metode Regula Falsi bisa disebut perpaduan antara Metode Bisection dan Metode Secant.

Sebenarnya, perbedaan algoritma Metode Bisection dan Metode Regula Falsi hanya di satu baris saja, yaitu mengubah baris

\[p=\frac{a+b}{2}\]

menjadi

\[p = b - \frac{f(b)(b-a)}{f(b) - f(a)}\]

sesuai Metode Secant. Perhatikan bahwa Metode Secant biasanya membutuhkan dua tebakan awal yang tidak harus sama dengan batasan interval, sedangkan Metode Regula Falsi secara otomatis menggunakan kedua batasan interval \([a,b]\) sebagai dua tebakan awal.

Untuk pembuatan kode Metode Regula Falsi, kami serahkan ke kalian. Gampang, kok! Tinggal mengubah beberapa baris saja (baris yang menentukan nilai \(p\) yang baru) pada kode Metode Bisection, yaitu mengambil baris tersebut dari kode Metode Secant, kemudian menyesuaikan kedua tebakan awal menjadi kedua batasan interval.

Seperti Metode Secant, Metode Regula Falsi juga memiliki order of convergence = \(\phi \approx 1.618\).

Apa itu barisan? (penjelasan tanpa kode)

Suatu “barisan” (sequence) adalah sekumpulan angka yang berurut. Artinya, pada suatu barisan, ada yang bisa disebut angka pertama (atau suku pertama), angka kedua (suku kedua), angka ketiga (suku ketiga), dan sebagainya. Banyaknya suku bisa berhingga maupun tak terhingga.

Suku-suku pada suatu barisan itu bisa saja ditentukan secara manual atau sesuka hati, atau bisa juga menggunakan rumus. Intinya, suku-suku suatu barisan itu bisa diperoleh dari manapun, bahkan dari hasil iterasi metode numerik (\(p_0\), \(p_1\), \(p_2\), \(p_3\), …) juga bisa.

Oleh karena itu, contoh barisan berhingga adalah hasil iterasi fixed-point, misalnya dengan \(g(x) = 1 + \frac{1}{x}\), tebakan awal \(p_0 = 2\), dan batas toleransi \(10^{-7}\):

\((1.5, 1.6666666666666665, 1.6, 1.625, 1.6153846153846154, 1.619047619047619, \dots, 1.6180339631667064)\)

Proses tersebut berakhir setelah 17 iterasi, sehingga barisan tersebut memiliki 17 suku.

Barisan tersebut bisa diberi nama, seperti \(p_n\) dengan \(n = 1, 2, 3, \dots, 17\), yang bisa dituliskan \(\left\{p_n\right\}_{n=1}^{17}\) dengan kurung kurawal.

Contoh barisan tak berhingga adalah barisan aritmetika dan barisan geometri, seperti:

\[(-5, -2, 1, 4, 7, 10, 13, 16, 19, \dots)\]

\[\left(16, 8, 4, 2, 1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \dots\right)\]

Barisan tak berhingga dengan nama \(p_n\) yang mulai dari suku \(n=1\) bisa ditulis \(\left\{p_n\right\}_{n=1}^{\infty}\) dengan kurung kurawal, atau singkatnya \((p_n)\) saja dengan kurung biasa (dengan begitu, biasanya ada asumsi bahwa barisan tersebut tak berhingga).

Metode Aitken

Alexander Aitken menemukan bahwa, untuk sembarang barisan (termasuk sembarang metode numerik) yang memiliki kekonvergenan linier, untuk nilai \(n\) yang besar, berlaku

\[\frac{p_{n+1} - p}{p_n - p} \approx \frac{p_{n+2} - p}{p_{n+1} - p}\]

di mana \(p\) adalah nilai yang ingin dicari, sedangkan \(p_n\), \(p_{n+1}\), dan \(p_{n+2}\) adalah tiga suku barisan (atau tiga hasil aproksimasi) berturut-turut. Artinya, perbandingan error (error ratio) antar dua pasang hasil iterasi (diperoleh dari tiga hasil iterasi berturut-turut) menjadi kurang lebih sama. Dengan manipulasi aljabar, diperoleh

\[p \approx p_n - \frac{\left(p_{n+1} - p_n\right)^2}{p_{n+2} - 2p_{n+1} + p_n}\]

seolah-olah ada jalur pintas untuk langsung mendapatkan nilai yang ingin dicari.

Tentu saja, sebelum menggunakan rumus ini, kita perlu menemukan tiga hasil aproksimasi pertama, yaitu \(p_0\), \(p_1\), dan \(p_2\). Kemudian, barulah kita tentukan \(p_3\) menggunakan rumus Aitken (hasil rumus Aitken biasa disebut \(\hat{p}_n\), sehingga bisa ditulis \(p_3 = \hat{p}_0\), karena perhitungan \(p_3\) memanfaatkan \(p_0\), \(p_1\) dan \(p_2\)).

Variabel \(\hat{p}\) biasa disebut p-hat atau p-cap (kata “hat” atau “cap” artinya topi).

Apabila kita definisikan \(\Delta p_n = p_{n+1} - p_n\) dan \(\Delta^2 p_n = p_{n+2} - 2p_{n+1} + p_n\), rumus Aitken bisa ditulis

\[\hat{p}_n = p_n - \frac{(\Delta p_n)^2}{\Delta^2 p_n}\]

sehingga teknik ini biasa disebut Aitken’s delta-squared (\(\Delta^2\)) method.

Catatan: dalam pembahasan metode Aitken/Steffensen, penulisan \(\Delta^2\) BUKAN berarti \((\Delta)^2\). Itu hanya penulisan saja.

Secara umum, apabila kita punya suku-suku suatu barisan yang berturut-turut yaitu \(p_1, p_2, p_3, \dots, p_{k-3}, p_{k-2}, p_{k-1}, p_{k}\), maka rumus Aitken bisa digunakan untuk menentukan \(\hat{p}_1, \hat{p}_2, \hat{p}_3, \dots, \hat{p}_{k-3}, \hat{p}_{k-2}\), yang semuanya merupakan aproksimasi nilai yang lebih akurat untuk hasil konvergen dari barisan tersebut (dengan asumsi kekonvergenan linier).

Perhatikan: - Kita hanya bisa berhenti sampai \(\hat{p}_{k-2}\), karena perhitungannya membutuhkan \(p_{k-2}\), \(p_{k-1}\) dan \(p_k\). - Harus ada minimal 3 suku yang diketahui, artinya \(k \ge 3\).

def Aitken(p):
    k = len(p)
    if k < 3:
        return "Maaf, dibutuhkan minimal 3 suku yang diketahui."
    
    # kalau lanjut ke sini, artinya k >= 3
    list_phat = []
    for i in range(k-2):
        Delta = p[i+1] - p[i]
        DeltaSquared = p[i+2] - 2 * p[i+1] + p[i]
        phat = p[i] - (Delta)**2 / DeltaSquared
        list_phat.append(phat)
    return list_phat

try:
    # input suatu list
    p = eval(input("Masukkan list suku-suku yang diketahui: "))
except:
    print("Maaf, terjadi error. Harap masukkan list dengan benar.")
else: # kalau tidak terjadi error 
    print("Berikut hasil metode Aitken:") 
    print(Aitken(p))
Masukkan list suku-suku yang diketahui: [2, 1.5, 1.6666666666666665, 1.6, 1.625]
Berikut hasil metode Aitken:
[1.625, 1.619047619047619, 1.6181818181818182]

Metode Steffensen: Penerapan Metode Aitken pada Metode Fixed Point

Aitken hanya menemukan rumus. Johan Frederik Steffensen menemukan bahwa, karena Metode Fixed-Point memiliki kekonvergen linier, metode Aitken bisa digunakan untuk mempercepat Metode Fixed-Point.

Secara umum, apabila kita berselang-seling antara menggunakan suatu metode dan rumus Aitken (misalnya setelah memperoleh tiga hasil aproksimasi), kita dapat mempercepat kekonvergenan (accelerating convergence), seolah-olah order of convergence menjadi lebih besar dari 1. Namun, bagaimana cara selang-selingnya?

Menurut Steffensen, rumus Aitken bisa digunakan tiap tiga iterasi fixed-point, yaitu untuk \(p_3\), \(p_6\), \(p_9\), dan seterusnya.

Kita bisa memodifikasi rumus Aitken dengan menggeser indeks \(n\), yaitu menukar \(n\) dengan \(n-3\), untuk mendapatkan rumus iterasi:

\[\hat{p} = p_{n-3} - \frac{\left(p_{n-2} - p_{n-3}\right)^2}{p_{n-1} - 2p_{n-2} + p_{n-3}}\]

dan dalam hal ini, kita juga bisa mendefinisikan \(\Delta_1 = p_{n-2} - p_{n-3}\) dan \(\Delta_2 = p_{n-1} - 2p_{n-2} + p_{n-3}\) untuk mendapatkan bentuk:

\[\hat{p} = p_{n-3} - \frac{(\Delta_1)^2}{(\Delta_2)}\]

from numpy import sin, cos, tan, log, exp, sqrt

formula = input("Masukkan formula fungsi: ")

def g(x):
    return eval(formula)

def Steffensen(p0, tolerance):
    # list semua nilai p agar mudah diakses
    list_p = [p0]

    # nilai sementara
    abs_error = tolerance + 1 

    iterasi = 1 # penghitung banyaknya iterasi
    while abs_error >= tolerance:
        if iterasi % 3 == 0: # untuk kelipatan tiga, gunakan rumus Aitken
            pn_3 = list_p[iterasi - 3] # p_(n-3)
            pn_2 = list_p[iterasi - 2] # p_(n-2)
            pn_1 = list_p[iterasi - 1] # p_(n-1)
            Delta1 = pn_2 - pn_3
            Delta2 = pn_1 - 2 * pn_2 + pn_3
            pn = pn_3 - (Delta1)**2 / Delta2
        else: # selain kelipatan 3, gunakan fixed point
            pn_1 = list_p[iterasi - 1]
            pn = g(pn_1)
        
        list_p.append(pn)
        abs_error = abs( pn - pn_1 )
        iterasi += 1
    
    # return bukan hanya p, tetapi juga banyaknya iterasi
    return pn, iterasi

starting_point = eval(input("Masukkan titik awal iterasi: "))
tolerance = eval(input("Masukkan batas toleransi: "))

p_steffensen, i_steffensen = Steffensen(starting_point, tolerance)

print("Metode Steffensen")
print("Hasil: " + str(p_steffensen))
print("setelah banyaknya iterasi: " + str(i_steffensen))

print("Bandingkan banyaknya iterasi dengan hasil Metode Fixed-Point biasa:")

print(FixedPoint(starting_point, tolerance))
Masukkan formula fungsi: 1 + 1/x
Masukkan titik awal iterasi: 2
Masukkan batas toleransi: 10**(-7)
Metode Steffensen
Hasil: 1.618033988749648
setelah banyaknya iterasi: 11
Bandingkan banyaknya iterasi dengan hasil Metode Fixed-Point biasa:
+---------+--------------------+
| iterasi |    Aproksimasi     |
+---------+--------------------+
|    1    |        1.5         |
|    2    | 1.6666666666666665 |
|    3    |        1.6         |
|    4    |       1.625        |
|    5    | 1.6153846153846154 |
|    6    | 1.619047619047619  |
|    7    | 1.6176470588235294 |
|    8    | 1.6181818181818182 |
|    9    | 1.6179775280898876 |
|   10    | 1.6180555555555556 |
|   11    | 1.6180257510729614 |
|   12    | 1.6180371352785146 |
|   13    | 1.6180327868852458 |
|   14    | 1.618034447821682  |
|   15    | 1.618033813400125  |
|   16    | 1.6180340557275543 |
|   17    | 1.6180339631667064 |
+---------+--------------------+
1.6180339631667064