Logistics and Agenda

Homework 3 due Monday Sept 21 at 11:59pm

Agenda

• This week: Root-finding, numerical integration
• Next week: Random numbers and distributions

Today

• Root finding methods

Midterm 1:

Do it at home part

• You will be given a task which would require to splice couple covered concepts

• Organize yourself in groups of 2 or less

  • give us group composition by Friday Sep 25

• You have almost two weeks to work on assignment

  • prepare code
  • report

• Midterm write up and code is due Sunday October 4th 11:59pm

• You may use AI freely to write code, debug, and test.

• You are responsible for the code you submit.

• We expect you to understand what the code is doing.

• You may not use AI to write the report.

  • Report should be done by you and only one per group is needed

In class part

• on Monday October 5th, with computers but no AI
  • we will give you modified assignment which would be similar
  • it will require modification or reuse of prepared code

Definition and canonical form

We want to find root of equation ($x_0$) which satisfies

$$f(x_0) = 0$$

Often it would work a bit different

$$h(x_0) = g(x_0)$$

Which can be moved to canonical form

$$f(x_0) = h(x_0) - g(x_0) = 0$$

Example

$$2\sin(x_0) = 1 \to 2\sin(x_0) - 1 = 0$$

Trial and error method

One can try to get the solution by just guessing with a hope to hit the solution. This is not highly scientific.

However, each guess can provide some clues.

A general algorithm

• decide on required precision (stop conditions)
• make a guess ($x_i$)
• make intelligent new guess ($x_{i+1}$) which uses available information
• repeat as long as $|f(x_{i+1})| > \varepsilon_f$ and $|x_{i+1} - x_i| > \varepsilon_x$

Bisection method pseudocode

Bisection python implementation

AI uses prompt in bisection_ai_assignment.txt

#!/usr/bin/env python3
"""
Bisection method implementation.

The function bis(f, x_m, x_p) finds a root of f in the interval [x_m, x_p]
assuming f(x_m) < 0 and f(x_p) > 0.

Parameters:
    f: callable, function whose root is to be found.
    x_m: float, lower bound of the bracket (f(x_m) < 0).
    x_p: float, upper bound of the bracket (f(x_p) > 0).
    eps_f: float, tolerance for function value (default 1e-6).
    eps_x: float, tolerance for interval width (default 1e-6).

Returns:
    x0: float, approximation of the root.
    f(x0): float, function value at the root approximation.
"""

def bis(f, x_m, x_p, eps_f=1e-6, eps_x=1e-6):
    """Bisection method to find a root of f in the bracket [x_m, x_p].

    The algorithm iteratively narrows the bracket until both the function
    value and the bracket width are within the specified tolerances.

    Args:
        f: callable, function whose root is to be found.
        x_m: float, lower bound of the bracket (f(x_m) < 0).
        x_p: float, upper bound of the bracket (f(x_p) > 0).
        eps_f: float, tolerance for function value (default 1e-6).
        eps_x: float, tolerance for interval width (default 1e-6).

    Returns:
        x0: float, approximation of the root.
        f(x0): float, function value at the root approximation.
    """
    x_m= x_m
    x_p= x_p
    while True:
        x_g = (x_p+ x_m) / 2.0
        f_g = f(x_g)
        if abs(f_g) <= eps_f and abs(x_p- x_g) <= eps_x:
            return x_g, f_g
        if f_g > 0:
            x_p= x_g
        else:
            x_m= x_g

if __name__ == "__main__":

# Simple test: find sqrt(2) as root of f(x) = x^2 - 2

def f(x):
        return x * x - 2.0
    root, val = bis(f, 1.0, 2.0)
    print(f"Root approximation: {root}")
    print(f"f(root) = {val}")
    print(f"x mismatch = {abs(root - 2**0.5)}")

Solution convergence expression

We say that algorithm has defined convergence if it is possible to express

$$\lim_{k\to\infty} (x_{k+1}-x_0)=c (x_k-x_0)^m$$

Where $x_0$ is true root of the equation, $c$ is some constant, and $m$ is the order of convergence.

The best algorithms have quadratic convergence, i.e. $m=2$

• The bisection algorithm has the linear rate of convergence: ($m=1$) and $c=1/2$

Generally the speed of the algorithm is related to its convergence order. However, other factors may affect the speed.

• For some algorithm it is hard to define convergence parameters

Newton-Raphson method

Numerical derivative of a function

Mathematical definition

$$f'(x)=\lim_{h \to 0}\frac{f(x+h)-f(x)}{h}$$

The initial intent is to calculate it at very small $h$.

Remember about round off errors

For computers with $h$ small enough $f(x+h)-f(x)=0$.

Let’s be smarter. Recall Taylor series expansion

$$f(x+h)=f(x)+\frac{f'(x)}{1!}h+\frac{f''(x)}{2!}h^2+\cdots$$

So we can see

$f_c'(x)=$$\frac{f(x+h)-f(x)}{h}$$= f'(x)$$+\frac{f''(x)}{2}h +\cdots$

Here computed approximation and algorithm error.

There is a range of optimal $h$ when both the round off and the algorithm errors are small.

Derivative via Forward Difference (fd)

$$f_{fd}'(x) = \frac{f(x+h)-f(x)}{h}$$

Algorithm error for small $h$

$$\varepsilon_{fd} \approx \frac{f''(x)}{2} h$$

This is quite bad since error is proportional to $h$.

Derivative via Central Difference (cd)

This is much better way (though it use an extra computation for a function)

$$f_c'(x) = \frac{f(x+h)-f(x-h)}{2h}$$

Algorithm error $$\varepsilon_{cd} \approx \frac{f'''(x)}{6}h^2$$