---
title: 數學軟體實作 - Root finding
tags: 2020 Fall - 數學軟體實作
GA: G-77TT93X4N1
---

# Lecture - Root finding

## 數列極限
已知 $\{a_n\}$ 是個收斂數列, 欲判斷數列極限
$$
\lim_{n\to\infty} a_n = L.
$$
想法就是看 $n$ 很大的時候 $a_n$ 的值是多少, 如果沒有變動太多我們就覺得他收斂了!

也就是說, 我們會設定一個容許值(tolerance), 如果 $|a_{n+1}-a_{n}|$ 小於這個數字, 我們就覺得差不多了, 並且覺得 $L \approx a_{n+1}$.

以下我們求 $a_n = \sin(\frac{1}{n})$的極限試試:

% 求數列極限範例
tol = 10^(-5);                       % 設定容許值

difference = 100;                    % 前後項的差(初始值)
ii=1;                                % 第一步
a0 = sin(1/ii);                      % a0 = 前一項(a_n)
while difference > tol               % 如果前後項差大於容許值就繼續做
    ii = ii+1;                       % 下一步
    a1 = sin(1/ii);                  % a1 = 下一項(a_{n+1})
    difference = abs(a1-a0);         % 計算前後項的差
    a0 = a1;                         % 設定下個迴圈的初始值為 a_1
    
    if(ii>10^5)                     % 設定跳脫次數
        disp('****WARNING****')      % 顯示示警及跳脫次數
        disp(['Stop!! Too many iterations. ', 'Have not found limit yet.'])
        break
    end
end
disp(['# of iterations= ', num2str(ii)])
disp(['L= ', num2str(a1)])           % 顯示估計出的極限值

---

#### Remark
容許值(tolerance) 越小, 估計出的極限值越準. 不過**容許值不等於誤差值**.

---

### Local function

如果我們常會需要算 $a_{n+1} = f(n)$ 這種樣子數列的極限, 可以用 local function 的方式定義 $f(x)$, 這樣以後若需要找數列極限只需要改 $f$ 的定義就好. 不用在函數裡翻來翻去看哪些地方需要改.

#### Remark
script 中的 local function 需要在整個 script 的最後, 並且每個 local function 都要有 `end` 做結尾

%
% 求數列極限範例
%
% 數列 a_{n+1} = f(a_n), 函數 f 定義在此檔案最後的 function y=f(x)
%
tol = 10^(-5);                       % 設定容許值

difference = 100;                    % 前後項的差(初始值)
ii=1;                                % 第一步
a0 = f(ii);                        % a0 = 前一項(a_n)
while difference > tol
    ii = ii+1;                       % 下一步
    a1 = f(ii);                    % a1 = 下一項(a_{n+1})
    difference = abs(a1-a0);         % 計算前後項的差
    a0 = a1;                         % 設定下個迴圈的初始值為 a_1
    
    if(ii>10^5)                     % 設定跳脫次數
        disp('****WARNING****')      % 顯示示警及跳脫次數
        disp(['Stop!! Too many iterations. ', 'Have not found limit yet.'])
        break
    end
end
disp(['# of iterations= ', num2str(ii)])
disp(['L= ', num2str(a1)])           % 顯示估計出的極限值

function y = f(x)                    % 定義函數 f(x)
y = sin(1./x);
end

---

### Exercise

考慮數列 $a_{n+1} = 10\sin(a_n)$, 試用不同初始值 $a_1$ 求其極限.

---

## 函數實根個數

給定一個函數 $f(x)$, $x\in[a, b]$, 我們想要把它所有的實根求出來. 也就是找到所有 $x_r$ 使得 $f(x_r)=0$.

不過首先我們需要先知道有幾個根. 我們用暴力法解決這個問題.

先將 $[a, b]$ 區間分成 $N$ 等分, 這樣我們會得到 $a=x_0<x_1<\cdots<x_{N+1}=b$ 這些區間. 接著我們用中間值定理來看每個區間內有沒有函數的實根. 如果以下式子成立:
$$
f(x_i)\cdot f(x_{i+1})<0,
$$
那我們就知道在 $[x_i, x_{i+1}]$ 這區間內至少有一個實根.

---

### Exercise

找出以下函數有多少個實根:
$$
f(x) = \frac{\sin(x)}{x} - \frac{1}{10}, \quad x\in [-10, 10]
$$

#### code structure - example

x = linspace(a, b, N+1);
num = 0;
for ii=1:N
    if f(x(ii))*f(x(ii+1))<0
        num = num+1;

return num

---

## 函數求根 - 固定點迭代

假設我想要求 $f(x) = \frac{\sin(x)}{x} - \frac{1}{10}$ 這個函數的實根. 我可以定義一個數列如下:
$$
\frac{\sin(a_n)}{a_{n+1}} -\frac{1}{10}=0, \quad \rightarrow \quad a_{n+1} = 10\sin(a_n).
$$
我們知道如果這數列收斂, 那收斂的值就會是 $f(x)$ 的實根.
* **不過它不一定會收斂.**

一般而言, 我們可以定義所謂的固定點迭代
$$
a_{n+1} = f(a_n).
$$
如果 ${a_n}$ 收斂了, 那他就是函數 $f$ 的固定點, $a = f(a)$.

### Fixed point theorem
如果有個函數 $f$ 滿足
$$
|f'(x)|\le \alpha<1, \quad x\in[a, b],
$$
那這個函數就有唯一的固定點, 並且固定點迭代一定會收斂.

### Newton's iteration
假設我們想解 $g(x)=0$, 那我們可以定義另一個函數 $F$:
$$
F(x) = x - \frac{g(x)}{g'(x)},
$$
* **$F$ 的固定點 = $g$ 的實根**

那這樣 $g$ 的實根會是 $F$ 的固定點, 並且在每個 $g$ 的實根附近一定有個區間使得以下固定點迭代會收斂, 並且收斂的很快:
$$
a_{n+1} = a_n - \frac{g(a_n)}{g'(a_n)}.
$$
* **不過我們不知道收斂區間的大小**

延伸閱讀 - [Te-Sheng Lin - Fixed point iteration](https://teshenglin.github.io/post/2019_fixed_point/)

* 牛頓法最大迭代次數設定為**10**

---

## Exercise - Root finding

試寫一函數求出 $\frac{\sin(x)}{x}=c$ 所有的實數解.

此函數 input 為 $c$, output 為所有實數解.

* $\left|\frac{\sin(x)}{x}\right| \le \frac{1}{|x|}$, 所以實數解只會落在 $\left[-\frac{1}{|c|}, \frac{1}{|c|}\right]$ 這個區間
* $|x|<\epsilon$ 需要單獨定義, 否則會出錯
* 若用牛頓法求根, 需驗證求出的解是否落在區間內
* 若用牛頓法求根10步不收斂, 就再將區間切細