Search code examples
pythonpseudocodenumerical-methods

Translate pseudocode into python (secant method)


I am working on this assignment:

First, implement the f-function defined by: f(x)= exp(x)-sin(x) closest to zero.

Second, implement the Secant method on page 95 and use it to find the root of the f-function, given the input values x0 = -3.5 and x1 = -2.5

Add the following
- an absolute test: abs(f(x) ) < epsilon
- a relative test: abs(x^k - x^{k-1})/ abs(x^{k}) \leq delta
- a maximum iteration guard: k < iter_max

In each iteration print out the iteration number k, the value of current root and current f-value. Print float numbers with 20 digits.

This is the code I have to complete:

import numpy as np
from math import exp, sin
import matplotlib.pyplot as plt

def f(x: float) -> float:
    return

def secant(x0: float, x1: float, f, epsilon: float, delta: float, iter_max: int) -> float:
    return

This is the pseudocode from page 95:

input: x_0, x_1, delta, epsilon, iter_max
fx_0 <- f(x_0); fx_1 <- f(x_1)
output: 0, x_0, fx_0
output: 1, x_1, fx_1
for k = 2 to iter_max do
    if |fx_0| > |fx_1| then
        x_0 <-> x_1; fx_0 <-> fx_1
    end if
    s <- (x_1 - x_0) / (fx_1 - fx_0)
    x_1 <- x_0
    fx_1 <- fx_0
    x_0 <- x_0 - fx_0 * s
    fx_0 <- f(x_0)

    output: k, x_0, fx_0
    if |fx_0| < epsilon or |x_1 - x_0| < delta then stop
end do

And this is my own attempt:

def f(x: float) -> float:
    return exp(x) - sin(x) == 0

def secant(x0: float, x1: float, f, epsilon: float, delta: float, iter_max: int) -> float:
    fx0 = f(x0)
    fx1 = f(x1)
    return 0, x0, fx0
    return 1, x1, fx1

    for k in range(2, iter_max):
        if abs(fx0) > abs(fx1):
            x0 = x1
            x1 = x0
            fx0 = fx1
            fx1 = fx0

            s = (x1 - x0) / (fx1 - fx0)
            x1 = x0
            fx1 = fx0
            x0 = x0 - fx0 * s
            fx0 = f(x0)

        return k, x0, fx0

        if abs(fx0) < epsilon or abs(x**k - x**(k - 1))/ abs(x**(k))  <= delta:
            break

If I follow my code with

root = secant(-3.5, -2.5, f, 0.00000000001, 0.000001, 10) 
print(root) 

I get: (0, -3.5, False). So it doesn't actually do any iterations. How can I fix it?

Edit:
Photo of pseudocode enter image description here

Here: a=x_0, b=x_1 and M =iter_max.
I would like the output to be something like this: enter image description here


Solution

  • Your implementation has the followings errors:

    • The first is that the function f must return the value of the function from which you want to get the root, you should not compare it with zero.

    • The second error is caused by not reading the algorithm's objective and understanding its procedure, if it says:


    output: 0, x_0, fx_0
    output: 1, x_1, fx_1
    

    indicates that it is the result when iter_max is 0 or 1 respectively.

    • Third, in the form of exchanging values is inadequate, you must use a pivot since otherwise erase the relevant information

    • Another thing I do not understand is because you execute the following: abs(x**k - x**(k - 1))/ abs(x**(k)) instead of abs(x1 - x0).

    So correcting both errors you get the following code:

    def f(x: float) -> float:
        return exp(x) - sin(x)
    
    def secant(x0: float, x1: float, f, epsilon: float, delta: float, iter_max: int) -> float:
        fx0 = f(x0)
        fx1 = f(x1)
        if iter_max == 0:
            return 0, x0, fx0
    
        elif iter_max == 1:
            return 1, x1, fx1
    
        for k in range(2, iter_max):
            if abs(fx0) > abs(fx1):
                tmp =  x0
                x0 = x1
                x1 =tmp
                tmp = fx0
                fx0 = fx1
                fx1 = tmp
    
            s = (x1 - x0) / (fx1 - fx0)
            x1 = x0
            fx1 = fx0
            x0 = x0 - fx0 * s
            fx0 = f(x0)
            if abs(fx0) < epsilon or abs(x1 - x0)  <= delta:
                break
        return k, x0, fx0
    
    
    root = secant(-3.5, -2.5, f, 0.00000000001, 0.000001, 10) 
    
    print(root)
    

    Output:

    (5, -3.183063011933318, 4.7351012000262926e-14)