MFDS-NM-03: Solusi Numerik Persamaan Non-Linear

Outline:

• Referensi/Materi
• Referensi Video
• Code/Module
• Forum Diskusi

Referensi/Materi

• https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/9781118673515.app8
• http://www.math.pitt.edu/~trenchea/math1070/MATH1070_5_Rootfinding.pdf
• Materi Bahasa Indonesia

Referensi Video

• https://www.youtube.com/watch?v=pStQf71WwB4 [Bisection]
• https://www.youtube.com/watch?v=6vh8QP1GliY [Regula Falsi]
• https://www.youtube.com/watch?v=szQUIRPrAgQ [Newton-Rhapson]
• https://www.youtube.com/watch?v=hZ4l5wXRV4g [Newton-Rhapson]
• https://www.youtube.com/watch?v=-JY0oavOhfw [Secant 01]
• https://www.youtube.com/watch?v=DbZ45Ej8YMw [Secant 02]

Open Code in Google Colaboratory

Metode Numerik
Solusi Numerik Persamaan nonlinear
¶

(C) Taufik Sutanto - 2020

tau-data Indonesia ~ https://tau-data.id/mfds-nm-03/

III. Solusi Numerik Persamaan nonlinear:

1. Fixed Points Problem: x = g(x)
2. Bisection : c = (a+b)/2
3. A little fix to Bisection: Regula Falsi
4. Approximate Location of Roots
5. Newton Rhapson Method
6. Secant Method

* All with Convergence Rate, Error Analysis, Application, & Practice..

Permasalahan mencari Akar (Roots) dari suatu persamaan

• Pada bidang ilmu Matematika, Statistika, Fisika, Teknik, dsb sering dijumpai permasalahan seperti berikut:
• _Diberikan sebuah fungsi kontinu $f(x)$, tentukan nilai $r$ sedemikian sehingga $f(r)=0$.
• Permasalahan seperti diatas disebut sebagai "root finding problems"

Diskusi ¶

• Apa contoh aplikasi penting mencari akar?

Akar Persamaan (Roots of Equations): Contoh¶

• $r=-2, 3, -1$ adalah akar-akar dari persamaan: $x^4-3x^3-7x^2+15x=-18$
• Solusi tersebut dapat divalidasi secara analytic dengan menuliskan persamaannya sebagai: $(x+2)(x-3)^2(x+1)=0$
• Kita bahkan memiliki solusi akar persamaan untuk sembarang polinomial derajat dua.
• $ax^2+bx+c = 0 $ memiliki akar-akar $x_{(1,2)} = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$
• Namun sayangnya hanya sedikit sekali persamaan yang bisa kita selesaikan secara analytic.
• Bahkan kita tidak memiliki solusi analytic untuk persamaan yang sederhana seperti $x=cos(x)$

• Pada kesempatan ini kita akan mencoba mencari cara (numerik) yang secara iteratif akan menemukan solusi pendekatan untuk "sembarang" fungsi Non-linear.

• What about linear? ... :)

In [1]:

# Python code untuk Fungsi
def P(x=3):
    return x**2 - x -2
P(x=1)

Out[1]:

-2

In [2]:

def pungsi(x):
    return x**2

pungsi(2)

Out[2]:

4

In [3]:

#inline pylab
# Python code untuk plot suatu fungsi
import pylab, numpy

X = numpy.linspace(-3,3,100) # 100 titik antara -3 sampai 2
y = [P(x) for x in X] # nilai fungsi di setiap titik
y2 = [pungsi(x) for x in X]
pylab.plot(X,y,'blue') # "red", atau r* ro , dsb
pylab.plot(X,y2,'red') # "red", atau r* ro , dsb
pylab.plot(X,[0]*100,'k') # sumbu x, 'k'=hitam

Out[3]:

[<matplotlib.lines.Line2D at 0x2d06cb2e438>]

Secara umum terdapat 3 cara dalam menyelesaikan persamaan tidak linear:

1. Solusi Analytic (seperti contoh sebelumnya)
2. Solusi Grafik (biasanya untuk perkiraan/initial guess metode lain)
3. Solusi Numerik (yang akan kita bahas di Mata Kuliah ini)

Kasus I: Fixed Points

• Nilai-nilai x yang memenuhi x=f(x) untuk x∈ℜ (lihat gambar) disebut sebagai titik-titik tetap (fixed points).

Teorema 1:¶

• Misal $f\in C[a,b]$ dan $y=f(x)$, jika $a\leq y \leq b$ ketika $a<x<b$, maka $f$ memiliki fixed point pada interval $[a,b]$
• Jika $f$ memiliki turunan di $[a,b]$ dan terdapat $K$ sedemikian sehingga $|f'(x)|\leq K<1$ untuk setiap $x\in (a,b)$, maka $f$ memiliki tepat satu fixed point di $[a,b]$

** Ponder this ... Diskusi cakupan teorema

Teorema 2 (fixed point theorem):

• Misal f dan f′ kontinu di lingkungan δ dari P = (P−δ,P+δ)=(a,b) dan memuat satu fixed point unik P. Jika suatu iterasi dimulai dari Po∈(a,b), maka:
• jika |f′(x)|≤K<1∀a≤x≤b, maka iterasi xn=f(xn-1) akan konvergen ke P. P disebut sebagai "Attractive Fixed Point"
• jika |f′(x)|>1∀a≤x≤b, maka iterasi xn=f(xn-1) tidak akan konvergen ke P. P disebut sebagai "Repulsive Fixed Point"

image Source: https://www.mathworks.com/matlabcentral/fileexchange/69171-the-general-iteration-method-fixed-point-iteration-method

Corolarry¶

• Misal kondisi pertama teorema 2 terpenuhi, maka $|P-P_n|\leq \frac{K^n|P_1-P_0|}{1-K}$

• Mengapa corolarry ini penting?

• Bagaimana menangani tentang nilai $K$ yang tidak diketahui?

Contoh:¶

• $f(x) = 1+x-x^2/4$ ==> Solusi P = {-2, 2} --Verify this!
• $f'(x) = 1-x/2$

In [4]:

# Gambarnya
def f(x):
    return 1-(x**2/4) # Kenapa fungsinya jadi seperti ini?

X = numpy.linspace(-2,2,100) # 100 titik antara -3 sampai 2
y = [f(x) for x in X] # nilai fungsi di setiap titik
pylab.plot(X,y,'r') # "red", atau r* ro , dsb
pylab.plot(X,[0]*100,'k') # sumbu x, 'k'=hitam

Out[4]:

[<matplotlib.lines.Line2D at 0x2d06ebe8f98>]

In [5]:

# Contoh Iterasi Titik Tetap
def f(x):
    return 1+x-(x**2/4)

n = 7
xo = 1
for i in range(n):
    x = f(xo)
    print(x, end =', ')
    xo = x
# modify for function in page 51

1.75, 1.984375, 1.99993896484375, 1.9999999990686774, 2.0, 2.0, 2.0, 

In [6]:

xo = 1
x = f(xo)
toleransi = 10**-9
maxIter = 1000
iterasi =0
while (abs(x-xo)>toleransi) and iterasi<maxIter:
    print(x, end =', ')
    xo = x
    x = f(xo)
    iterasi += 1

1.75, 1.984375, 1.99993896484375, 1.9999999990686774, 

Diskusi :¶

• Apa sebaiknya "stopper" iterasinya?
• Bagaimana mendeteksi apakah iterasinya konvergen/divergen secara numerik?
• Bagaimana menduga error jika kita menggunakan n iterasi? (asumsikan konvergen)

Metode BiSection (Bolzano)¶

• Teorema: Asumsikan $f\in C[a,b]$ dan terdapat $r\in [a,b]$ sedemikian sehingga $f(r)=0$. jika $f(a)*f(b)<0$ dan $c$ adalah titik tengah $a$ dan $b$, maka $|r-c_n|\leq \frac{b-a}{2^{n+1}}$ untuk $n=1,2,3,...$ dan $\lim_{n\rightarrow \inf}c_n=r$

Source: https://commons.wikimedia.org/wiki/File:Bisection_anime.gif

In [7]:

# Contoh Aplikasi Bisection di Python
# di Buku Hal 58
import math

def f(x):
    return x*math.sin(x) - 1

a = 0
b = 2
n = 9 # iterasi

if f(a)*f(b)<0:
    for i in range(n):
        c = (a+b)/2 # THE BiSection
        print(a,b,c,f(c))
        if f(a)*f(c)<0:
            b = c
        else:
            a = c
else:
    print('Error, a dan b tidak mengapit akar')

0 2 1.0 -0.1585290151921035
1.0 2 1.5 0.4962424799060816
1.0 1.5 1.25 0.18623077419448286
1.0 1.25 1.125 0.015051043361482108
1.0 1.125 1.0625 -0.0718266313849869
1.0625 1.125 1.09375 -0.028361722663138855
1.09375 1.125 1.109375 -0.0066427748602908565
1.109375 1.125 1.1171875 0.004208034010642736
1.109375 1.1171875 1.11328125 -0.0012164904180159697

In [8]:

a = 0
b = 2
toleransi = 10**-6
iterasi = 0
if f(a)*f(b)<0:
    c = (a+b)/2
    while (abs(f(c))> toleransi ):
        iterasi = iterasi +1 # iterasi =+`
        c = (a+b)/2 # THE BiSection
        print('iterasi ',iterasi, a,b,c,f(c))
        if f(a)*f(c)<0:
            b = c
        else:
            a = c
else:
    print('Error, a dan b tidak mengapit akar')

iterasi  1 0 2 1.0 -0.1585290151921035
iterasi  2 1.0 2 1.5 0.4962424799060816
iterasi  3 1.0 1.5 1.25 0.18623077419448286
iterasi  4 1.0 1.25 1.125 0.015051043361482108
iterasi  5 1.0 1.125 1.0625 -0.0718266313849869
iterasi  6 1.0625 1.125 1.09375 -0.028361722663138855
iterasi  7 1.09375 1.125 1.109375 -0.0066427748602908565
iterasi  8 1.109375 1.125 1.1171875 0.004208034010642736
iterasi  9 1.109375 1.1171875 1.11328125 -0.0012164904180159697
iterasi  10 1.11328125 1.1171875 1.115234375 0.0014960036578040015
iterasi  11 1.11328125 1.115234375 1.1142578125 0.00013981310235822164
iterasi  12 1.11328125 1.1142578125 1.11376953125 -0.0005383247225386745
iterasi  13 1.11376953125 1.1142578125 1.114013671875 -0.0001992523031036919
iterasi  14 1.114013671875 1.1142578125 1.1141357421875 -2.971872073043169e-05
iterasi  15 1.1141357421875 1.1142578125 1.11419677734375 5.50474110863064e-05
iterasi  16 1.1141357421875 1.11419677734375 1.114166259765625 1.2664400200756987e-05
iterasi  17 1.1141357421875 1.114166259765625 1.1141510009765625 -8.527146514669681e-06
iterasi  18 1.1141510009765625 1.114166259765625 1.1141586303710938 2.0686302812933377e-06
iterasi  19 1.1141510009765625 1.1141586303710938 1.1141548156738281 -3.2292572573755507e-06
iterasi  20 1.1141548156738281 1.1141586303710938 1.114156723022461 -5.803132732129512e-07

Diskusi:¶

• Apa kelebihan dan kekurangan Bisection.
• Apa stopper-nya?
• Jika ada >1 akar diapit apakah tetap akan konvergen?

Regula Falsi

• Kita tau bahwa $m = \frac{f(b)-f(a)}{b-a}$, jika kita ambil 2 pasang titik : $(a,f(a))$ ~ $(b,f(b))$ dan $(b,f(b))$ ~ $(c,0)$ maka $\frac{f(b)-f(a)}{b-a}=\frac{0-f(b)}{c-b}$.
• sehingga kita mendapatkan cara menghitung $c$ dengan cara lain: $c= b - \frac{f(b)(b-a)}{f(b)-f(a)}$

Source: https://commons.wikimedia.org/wiki/File:Regula_falsi_method.png

In [9]:

a = 0
b = 2
toleransi = 10**-6
iterasi = 0
if f(a)*f(b)<0:
    c = (a+b)/2
    while (abs(f(c))> toleransi ):
        iterasi = iterasi +1 # iterasi =+`
        c = b - (f(b)*(b-a)) / (f(b)-f(a)) # THE RegulaFalsi
        print('iterasi ',iterasi, a,b,c,f(c))
        if f(a)*f(c)<0:
            b = c
        else:
            a = c
else:
    print('Error, a dan b tidak mengapit akar')

iterasi  1 0 2 1.0997501702946164 -0.020019210242675722
iterasi  2 1.0997501702946164 2 1.1212407359645027 0.00983461086237658
iterasi  3 1.0997501702946164 1.1212407359645027 1.1141611949626335 5.630358231867305e-06
iterasi  4 1.0997501702946164 1.1141611949626335 1.1141571430336825 3.0022619945668794e-09

In [10]:

# Regula Falsi di Python
# Buku Halaman 61

def f(x):
    return 2*x**2-5*x #x*math.sin(x) - 1

a = 1 #0
b = 5 #2
n = 4 # iterasi

if f(a)*f(b)<0:
    for i in range(n):
        c = b - (f(b)*(b-a)) / (f(b)-f(a)) # THE RegulaFalsi
        print(a,b,c,f(c))
        if f(a)*f(c)<0:
            b = c
        else:
            a = c
else:
    print('Error, a dan b tidak mengapit akar')
    
(2.105-2.5)/2.5

1 5 1.4285714285714284 -3.0612244897959187
1.4285714285714284 5 1.818181818181818 -2.479338842975208
1.818181818181818 5 2.105263157894737 -1.6620498614958432
2.105263157894737 5 2.285714285714286 -0.9795918367346932

Out[10]:

-0.158

In [11]:

import numpy as np


def f(x):
    return np.exp(x) -2 - x

a = 0.0
b = 4.0
n = 4 # iterasi

if f(a)*f(b)<0:
    for i in range(n):
        c = b - (f(b)*(b-a)) / (f(b)-f(a)) # THE RegulaFalsi
        print(a,b,c,f(c))
        if f(a)*f(c)<0:
            b = c
        else:
            a = c
else:
    print('Error, a dan b tidak mengapit akar')

0.0 4.0 0.08064816928306762 -0.996658720599791
0.08064816928306762 4.0 0.15941157777968584 -0.9865910235413711
0.15941157777968584 4.0 0.23582803134445873 -0.9698714447500276
0.23582803134445873 4.0 0.30947960743481406 -0.9467638250137362

Diskusi:¶

• Apa kelebihan dan kekurangan Regula Falsi?
• Apakah stopper Regula Falsi berbeda? Kenapa?

Not Bracketing, tapi hanya butuh satu titik: Newton-Rhapson¶

• Teorema: Misal $f\in C^2[a,b]$ dan terdapat $p\in [a,b]$ dimana $f(p)=0$. Jika $f'(p)\neq 0$, maka terdapat $\delta >0$ sedemikian sehingga barisan $f(p)$ berikut akan konvergen ke $p$:
• $p_n = p_{n-1}-\frac{f(p_{n-1})}{f'(p_{n-1})}$ ... (Newton-Rhapson Iteration)

Source: https://commons.wikimedia.org/wiki/File:NewtonIteration_Ani.gif

Newton-Rhapson untuk menghitung akar¶

• Corollary: Misal $A>0 \in \Re$ dan $p>0$ adalah nilai titik awal untuk pendekatan ke $\sqrt{A}$, maka
• $p_n = \frac{p_{n-1}+\frac{A}{p_{n-1}}}{2}$ konvergen ke $\sqrt{A}$ ... (Coba turunkan formula ini).

In [12]:

# Contoh Buku hal 74
A = 5
p = 2 # Po titik awal
n = 4 # Jumlah iterasi

print(p, end = ', ')
for i in range(n):
    p = 0.5 * (p + A/p)
    print(p, end = ', ')
    
print('\n Akar 5 = ',math.sqrt(5))
print('Error = ', abs(math.sqrt(5)-p))

2, 2.25, 2.236111111111111, 2.2360679779158037, 2.23606797749979, 
 Akar 5 =  2.23606797749979
Error =  0.0

Order of Convergence Newton:¶

• Definisi Order of Convergence $R$
• Misal $p_n \rightarrow p$, dan $e_n = p-p_n$ untuk $n\geq 0$. Jika terdapat $A\neq 0$ dan $R>0$ dan $\lim_{n\rightarrow\inf} \frac{|p-p_{n+1}|}{|p-p_n|^R} = \lim_{n\rightarrow\inf} \frac{|e_{n+1}|}{|e_n|^R}=A$, maka barisannya dikatakan konvergen ke $p$.
• $A$ disebut sebagai konstanta error asimtotik.
• Jika $R=1$ convergence rate-nya Linear, R=2 kuadratik.

Teorema Convergence Rate Newton-Rhapson¶

• Misal barisan newton $p_n$ konvergen ke $p$. jika p adalah (simple) root, convergence rate metode Newton adalah kuadratik dan
• $|e_{n+1}|\approx \frac{|f"(p)|}{2|f'(p)|}|e_n|^2$ untuk suatu nilai $n$ yang cukup besar.
• Jika $p$ adalah akar order $M \rightarrow$ konvergensinya linear dan
• $|e_{n+1}|\approx \frac{M-1}{M}|e_n|$

** "order of root":

• $p$ order $M$ jika: $f(p)=0, f'(p)=0, f"(p)=0, ..., f^{(M)}\neq 0$

Diskusi:¶

• Apa kelebihan dan kekurangan metode Newton?
• Apa kekurangan terbesar metode Newton?
• Bagaimana cara terbaik menghitung "A" pada Convergence rate Newton?

In [13]:

# Contoh Aplikasi Newton-Rhapson di Python
# Contoh di Buku Hal 77
def f(x):
    return x**3 + x**2 -3*x -3 # x**3-3*x+2

def df(x): # Turunan pertama f
    return 3*x**2 + 2*x - 3 #3*x**2-3

x = 1 # -2.4 # Initial point/guess
n = 9 # Jumlah iterasi

print(x, end = ', ')
for i in range(n):
    x = x - f(x)/(df(x))
    print(x, end = ', ')
    
abs(2.2-1.732)/1.732

1, 3.0, 2.2, 1.8301507537688444, 1.7377954531428215, 1.7320722915449542, 1.7320508078710555, 1.7320508075688774, 1.7320508075688772, 1.7320508075688774, 

Out[13]:

0.2702078521939955

Metode Secant¶

• Kelemahan utama metode Newton adalah ia memerlukan turunan (pertama) analitik fungsinya.
• Dalam aplikasinya bahkan sebenarnya yg dibutuhkan turunan kedua (Why?)
• Metode Secant mengganti turunan pertama dengan pendekatannya.
• Kenapa boleh? Karena turunan pertama di Metode Newton sebenarnya hanya digunakan sebagai "Arah". Apa Maksudnya?
• Turunan pertama fungsi dapat di aproksimasi dengan berbagai cara, cara yang paling mudah adalah menggunakan limit berikut:
• $f'(x) = \lim_{h\rightarrow 0} \frac{f(x+h)-f(x)}{h}$ ... sering disebut sebagai forward differencing
• Aproksimasi dilakukan dengan menggunakan suatu nilai $h$ yang relatif cukup kecil (misal $h=10^{-5}$)

In [14]:

# Contoh Aplikasi Newton-Rhapson di Python
# Contoh di Buku Hal 77
def dfs(x, h=10**-2): 
    return (f(x+h)-f(x))/h # Aproksimasi Turunan pertama f

x = 1 # -2.4 # Initial point/guess
n = 5 # Jumlah iterasi

print(x, end = ', ')
for i in range(n):
    x = x - f(x)/(dfs(x))
    print(x, end = ', ')

abs(2.2-1.732)/1.732

1, 2.9606882015587788, 2.1817616576450876, 1.8255016711899639, 1.7378198680514445, 1.7321097411084792, 

Out[14]:

0.2702078521939955

In [15]:

def f(x):
    return np.exp(x) -2 - x


def dfs(x, h=10**-2): 
    return (f(x+h)-f(x))/h # Aproksimasi Turunan pertama f

x = 4.0
n = 3 # Jumlah iterasi

print(x, end = ', ')
for i in range(n):
    x = x - f(x)/(dfs(x))
    print(x, end = ', ')

4.0, 3.097896825942087, 2.2958526267774717, 1.668309817179185, 

In [16]:

np.exp(-2)
# np.math.factorial(4)

Out[16]:

0.1353352832366127

In [17]:

def P(x):
    return np.exp(5)*(1+7*(x-1)+ (7**2*(x-1)**2)/np.math.factorial(2) + 
               (7**3*(x-1)**3)/np.math.factorial(3) + (7**4*(x-1)**4)/np.math.factorial(4) )
P(0)

Out[17]:

9108.85763992064

In [18]:

Error_Relatif = (np.exp(-2) - P(0))/np.exp(-2)
Error_Relatif

Out[18]:

-67304.86009854663

End of Week Module "Solusi Numerik Persamaan nonlinear "

---

Jangan lupa kolom komentar bukan untuk berdiskusi/bertanya, untuk keperluan tersebut silahkan klik tautan untuk menuju Forum terkait modul ini.

What's Next?

[Previous Lesson]

[Diskusi]

[Next Lesson]

[Kembali ke Kurikulum]

Referensi:

1. Johansson, R. (2018). Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib. Apress.
2. John H. Mathews, Numerical Methods for Mathematics, Prentice Hall, 1992 [Refernsi Utama]
3. Heath, M. T. (2018). Scientific computing: an introductory survey (Vol. 80). SIAM.
4. Conte, S. D., u0026amp; De Boor, C. (2017). Elementary numerical analysis: an algorithmic approach (Vol. 78). SIAM.