数值计算详细笔记（二）：非线性方程解法

2. Nonlinear Equation f(x) = 0

2.1 Bisection Method

2.1.1 Root-Finding Problem

Root-Finding or Zero-Finding Problem

Given a function $f(x)$f(x) in one variable $x$x, finding a root $x$x of an equation of the form $f(x)=0$f(x)=0.

Solution $x$x is called a root of equation $f(x)=0 $, or zero of function $f(x)$f(x).

2.1.2 Bisection Algorithm

Prejudgment

By the Intermediate Value Theorem, if

$f\in C[a,b], and \ f(a)*f(b)<0,$f∈C[a,b],and f(a)∗f(b)<0,

then there exists at least a point $x^*\in (a,b)$x∗∈(a,b), such that $f(x^*)=0$f(x∗)=0.

Algorithm Description

• INPUT: endpoints $a,b$a,b; tolerance $TOL$TOL; maximum number of iterations $N$N.

• OUTPUT: approximate solution $c$c or message of failure.

• Step $1$1: Set $k = 1,FA=f(a)$k=1,FA=f(a);

• Step $2$2: While $k\leq N$k≤N, do Steps $3～6$3～6

  • Step $3$3: Set $c= \displaystyle\frac{a+b}{2}$c=2a+b​; and compute $FC=f(c)$FC=f(c).
  • Step $4$4: If $FC=0$FC=0 or $\displaystyle\frac{|b-a|}{2}<TOL$2∣b−a∣​<TOL, then output $c$c, (Procedure complete successfully.) Stop!
  • Step $5$5: If $FA*FC<0$FA∗FC<0, then set $b=c$b=c; else set $a=c$a=c
  • Step $6$6: Set $k=k+1$k=k+1.

• Step $7$7: OUTPUT “Method failed after $N$N iterations.” STOP!

Other Stopping Criteria

Other Stopping Criteria for Iteration procedures with a given tolerance $\varepsilon >0:$ε>0:

$|p_n-p_{n-1}|<\varepsilon\\ \displaystyle\frac{|p_n-p_{n-1}|}{|p_n|} < \varepsilon \\ |f(p_n)|< \varepsilon$∣pn​−pn−1​∣<ε∣pn​∣∣pn​−pn−1​∣​<ε∣f(pn​)∣<ε

2.1.3 Convergence Analysis

Theorem

Suppose that $f\in C[a,b]$f∈C[a,b], and $f(a)*f(b)<0$f(a)∗f(b)<0. The Bisection method generates a sequence $\{{p_n}\}_1^\infty${pn​}1∞​ approximating a zero point $p$p of $f$f with

$|p_n-p|\leq \displaystyle\frac{b-a}{2^n},n\geq1.$∣pn​−p∣≤2nb−a​,n≥1.

Proof

By the procedure, we know that

$|b_1-a_1|=|b-a|,\\ |b_2-a_2|=\displaystyle\frac{|b_1-a_1|}{2}=\displaystyle\frac{|b-a|}{2},\\ ...\\ |b_n-a_n|=\displaystyle\frac{|b_{n-1}-a_{n-1}|}{2}=\displaystyle\frac{|b-a|}{2^{n-1}},$∣b1​−a1​∣=∣b−a∣,∣b2​−a2​∣=2∣b1​−a1​∣​=2∣b−a∣​,...∣bn​−an​∣=2∣bn−1​−an−1​∣​=2n−1∣b−a∣​,

Since $p_n=\displaystyle\frac{a_n+b_n}{2}$pn​=2an​+bn​​ and $p\in (a_n,p_n]$p∈(an​,pn​] or $p\in [p_n,b_n)$p∈[pn​,bn​) for all $n\geq 1$n≥1, it follows that

$|p_n-p|\leq \displaystyle\frac{|b_n-a_n| }{2}=\displaystyle\frac{b-a}{2^n}.$∣pn​−p∣≤2∣bn​−an​∣​=2nb−a​.

Convergence Rate

Since

$|p_n-p|\leq \displaystyle\frac{|b_n-a_n|}{2}=\displaystyle\frac{|b-a|}{2^n}$∣pn​−p∣≤2∣bn​−an​∣​=2n∣b−a∣​

The Sequence $\{p_n\}_{n=1}^\infty${pn​}n=1∞​ converges to $p$p with rate of convergence $O(\displaystyle\frac{1}{2^n})$O(2n1​), that is

$p_n=p+O(\displaystyle\frac{1}{2^n})$pn​=p+O(2n1​)

Other Property

Bisection is certain to converge, but does so slowly.

Given starting interval $[a,b]$[a,b], length of interval after $k$k iterations is $\displaystyle\frac{b-a}{2^k}$2kb−a​, so achieving error tolerance of $\varepsilon(\displaystyle\frac{(b-a)}{2^k}<\varepsilon)$ε(2k(b−a)​<ε) requires $k\approx[log_2^{\frac{b-a}{\varepsilon}}]$k≈[log2εb−a​​] iterations, regardless of function $f$f involved.

2.2 Fixed Point Method

2.2.1 Introduction

• Fixed point of given function $g:\mathbb{R}\rightarrow \mathbb{R}$g:R→R is value $x^*$x∗ such that $x^*=g(x^*)$x∗=g(x∗)
• Many iterative methods for solving nonlinear equations use fixed-point iteration scheme of form
  $x_{k+1}=g(x_k)$xk+1​=g(xk​)
  where fixed points for $g$g are solutions for $f(x)=0$f(x)=0.

Example

If $f(x)=x^2-x-2$f(x)=x2−x−2, it has two roots $x^*=2$x∗=2 and $x^*=-1$x∗=−1. Then fixed points of each of functions
$g(x)=x^2-2 \\ g(x) = \sqrt{x+2}\\ g(x) = 1+\displaystyle\frac{2}{x} \\ g(x) = \displaystyle\frac{x^2+2}{2*x-1}$g(x)=x2−2g(x)=x+2​g(x)=1+x2​g(x)=2∗x−1x2+2​
are solutions to equation $f(x)=0$f(x)=0.

2.2.2 Algorithm

Definition

• To approximate the fixed point of a function $g(x)$g(x), we choose an initial approximation $p_0$p0​, and generate the sequence $\{p_n\}^\infty_{n=0}${pn​}n=0∞​ by letting
  $\left\{ \begin{aligned} Given & \ \ \ \ p_0 \\ p_n = & g(p_{n-1}),n=0,1,..., \end{aligned} \right.${Givenpn​=​    p0​g(pn−1​),n=0,1,...,​
  for each $n\geq 1$n≥1.
• If the sequence $\{p_n\}^\infty_{n=0}${pn​}n=0∞​ converges to $p$p and $g(x)$g(x) is continuous, then we have
  $p = \lim_{n\rightarrow \infty}p_n=\lim_{n\rightarrow \infty}g(p_n)=g(\lim_{n\rightarrow \infty}p_n)=g(p).$p=n→∞lim​pn​=n→∞lim​g(pn​)=g(n→∞lim​pn​)=g(p).
  and a solution to $x=g(x)$x=g(x) is obtained.
• This technique described above is called fixed point iteration (or functional iteration).

Example

Pseudo-Code

• INPUT: Initial approximation $p_0$p0​, tolerance $TOL$TOL, Maximum number of iteration $N$N.

• OUTPUT: approximate solution $p$p or message of failure.

• Step $1$1: Set $n = 1$n=1.

• Step $2$2: While $n\leq N$n≤N, do Steps $3～5$3～5.

  • Step $3$3: Set $p= g(p_0)$p=g(p0​).
  • Step $4$4: If $|p-p_0|<TOL$∣p−p0​∣<TOL, then output $p$p; (Procedure complete successfully.) Stop!
  • Step $5$5: Set $n=n+1, p_0=p$n=n+1,p0​=p.

• Step $6$6: OUTPUT “Method failed after $N$N iterations.” STOP!

2.2.3 Existence and Uniqueness

Theorem: Sufficient Conditions for the Existence and Uniqueness of a Fixed Point

1. 【Existence】If $g(x)\in C[a,b]$g(x)∈C[a,b] and $g(x)\in [a,b]$g(x)∈[a,b] for all $x\in [a,b]$x∈[a,b], then $g(x)$g(x) has a fixed point in $[a,b]$[a,b]. ($g(x)\in [a,b]$g(x)∈[a,b] means $max\{g(x)\}\leq b,min\{g(x)\}\geq a$max{g(x)}≤b,min{g(x)}≥a)
2. 【Uniqueness】If, in addition, $g'(x)$g′(x) exists on $(a,b)$(a,b), and a positive constant $k<1$k<1 exists with $|g'(x)|\leq k$∣g′(x)∣≤k, for all $x\in (a,b)$x∈(a,b).

Meet the two conditions above means the fixed point in $[a,b]$[a,b] is unique.

Proof for Existence

1. If $g(a)=a$g(a)=a or $g(b)=b$g(b)=b, then $g(x)$g(x) has a fixed point at an endpoint.
2. Suppose not, then it must be true that $g(a)>a$g(a)>a and $g(b)<b$g(b)<b.
3. Thus the function $h(x)=g(x)-x$h(x)=g(x)−x is continuous on $[a,b]$[a,b], and we have
   $h(a)=g(a)-a>0,h(b)=g(b)-b<0.$h(a)=g(a)−a>0,h(b)=g(b)−b<0.
4. The Intermediate Value Theorem implies that there exists $p\in (a,b)$p∈(a,b) for $h(x)=g(x)-x$h(x)=g(x)−x which $h(p)=0$h(p)=0.
5. Thus $g(p)-p=0$g(p)−p=0, and $p$p is a fixed point of $g(x)$g(x).

The proving process above changes $g(x)=x$g(x)=x to $f(x)=g(x)-x$f(x)=g(x)−x, changing fixed point problem into zero point existing problem.

Proof for Uniqueness

Using contradiction method to prove this characteristic.

1. Suppose, in addition, that $|g'(x)\leq k\leq 1|$∣g′(x)≤k≤1∣ and that $p$p and $q$q are both fixed points in $[a,b]$[a,b] with $p\not=q.$p​=q.
2. Then by the Mean Value Theorem, a number $\xi$ξ exists between $p$p and $q$q. Hence, in $[a,b]$[a,b],
   $\displaystyle\frac{g(p)-g(q)}{p-q}=g'(\xi)$p−qg(p)−g(q)​=g′(ξ)
   exists.
3. Then
   $|p-q|=|g(p)-g(q)|=|g'(\xi)||p-q|\leq k*|p-q|<|p-q|,$∣p−q∣=∣g(p)−g(q)∣=∣g′(ξ)∣∣p−q∣≤k∗∣p−q∣<∣p−q∣,
   which is a contradiction.
4. So $p=q$p=q, and the fixed point in $[a,b]$[a,b] is unique.

2.2.4 Convergence Analysis

Theorem

Meet the two conditions above, we can find that, for any number $p_0\in[a,b]$p0​∈[a,b], the sequence $\{p_n\}^\infty_0${pn​}0∞​ defined by
$p_n=g(p_{n-1}), n\geq 1$pn​=g(pn−1​),n≥1
converges to the unique fixed point $p$p in $[a,b]$[a,b].

Proof

Using the Intermediate Value Theorem, we can prove the existence.

Using the fact that $|g'(x)|\leq k$∣g′(x)∣≤k and the Mean Value Theorem, we have
$|p_n-p|=|g(p_{n-1}-g(p))|=|g'(\xi)||p_{n-1}-p|\leq k*|p_{n-1}-p|,$∣pn​−p∣=∣g(pn−1​−g(p))∣=∣g′(ξ)∣∣pn−1​−p∣≤k∗∣pn−1​−p∣,
where $\xi \in (a,b)$ξ∈(a,b).

Applying this inequality inductively shows
$|p_n-p|\leq k*|p_{n-1}-p|\leq k^2*|p_{n-2}-p|\leq ...\leq k^n*|p_{0}-p|.$∣pn​−p∣≤k∗∣pn−1​−p∣≤k2∗∣pn−2​−p∣≤...≤kn∗∣p0​−p∣.

Since $k<1$k<1,
$\lim_{n\rightarrow\infty}|p_n-p|\leq \lim_{n\rightarrow\infty}k^n|p_0-p|=0,$n→∞lim​∣pn​−p∣≤n→∞lim​kn∣p0​−p∣=0,
and $\{p_n\}_0^\infty${pn​}0∞​ converges to $p$p.

Convergence Rate

Since
$|p_n-p|\leq k^n*|p_0-p|\leq k^n*|b-a|$∣pn​−p∣≤kn∗∣p0​−p∣≤kn∗∣b−a∣,
the Sequence $\{p_n\}_{n=0}^\infty${pn​}n=0∞​ converges to $p$p with rate of convergence $O(k^n)$O(kn) with $k< 1$k<1, that is
$p_n=p+O(k^n),k<1$pn​=p+O(kn),k<1.

2.2.5 Bounds for the Error

Corollary

If the solution for $g(x)=x$g(x)=x exists and is unique, then the bounds for the error involved in using $p_n$pn​ to approximate $p$p are given by
$|p_n-p|\leq k^n*max\{p_0-a,b-p_0\}$∣pn​−p∣≤kn∗max{p0​−a,b−p0​}
and
$|p_n-p|\leq \displaystyle\frac{k^n}{1-k}*|p_1-p_0|,$∣pn​−p∣≤1−kkn​∗∣p1​−p0​∣,
for all $n\geq 1$n≥1.

Proof

$|p_n-p_{n-1}|\leq |g(p_{n-1})-g(p_{n-2})|\leq k*|p_{n-1}-p_{n-2}|\leq ...\leq k^{n-1}|p_1-p_0|.$∣pn​−pn−1​∣≤∣g(pn−1​)−g(pn−2​)∣≤k∗∣pn−1​−pn−2​∣≤...≤kn−1∣p1​−p0​∣.

Let $m>n$m>n, using the theorem $|a+b|\leq |a|+|b|$∣a+b∣≤∣a∣+∣b∣, then we have
$|p_m-p_n|\leq |p_m-p_{m-1}|+|p_{m-1}-p_{m-2}|+...+|p_{n+1}-p_n|\leq (k^{m-1}+k^{m-2}+...+k^n)|p_1-p_0|\\ |p_m-p_m|\leq k^n*(1+k+...+k^{m-n-1})|p_1-p_0|\leq k^n*\displaystyle\frac{1-k^{m-n-1}}{1-k}*|p_1-p_0|$∣pm​−pn​∣≤∣pm​−pm−1​∣+∣pm−1​−pm−2​∣+...+∣pn+1​−pn​∣≤(km−1+km−2+...+kn)∣p1​−p0​∣∣pm​−pm​∣≤kn∗(1+k+...+km−n−1)∣p1​−p0​∣≤kn∗1−k1−km−n−1​∗∣p1​−p0​∣.

Let $m\rightarrow\infty$m→∞, we have
$\lim_{m\rightarrow\infty}|p_m-p_n|=|p-p_n|\leq \displaystyle\frac{k^n}{1-k}*|p_1-p_0|.$m→∞lim​∣pm​−pn​∣=∣p−pn​∣≤1−kkn​∗∣p1​−p0​∣.

2.3 Newton’s Method

2.3.1 Introduction

Status

The Newton-Raphson (or simply Newton’s) method is one of the most powerful and well-known numerical methods for solving a root-finding problem
$f(x)=0.$f(x)=0.

Rough Description

(1) Suppose that $f\in C^2[a,b]$f∈C2[a,b], and $x^*$x∗ is a solution of $f(x)=0$f(x)=0.

(2) Let $\hat{x}\in [a,b]$x^∈[a,b] be an approximation to $x^*$x∗ such that $f'(\hat{x})\not=0$f′(x^)​=0 and $|\hat{x}-x^*|$∣x^−x∗∣ is “small”.

Consider the first Taylor polynomial for $f(x)$f(x) expanded about $\hat{x}$x^,
$f(x)=f(\hat{x})+(x-\hat{x})f'(\hat{x})+\displaystyle\frac{(x-\hat{x})^2}{2}f''(\xi(x)).$f(x)=f(x^)+(x−x^)f′(x^)+2(x−x^)2​f′′(ξ(x)).
where $\xi(x)$ξ(x) lies between $x$x and $\hat{x}$x^.

Thus, consider $f(x^*)=0$f(x∗)=0, and gives
$0=f(x^*)\approx f(\hat{x})+(x^*-\hat{x})f'(\hat{x}).$0=f(x∗)≈f(x^)+(x∗−x^)f′(x^).

Solution for finding $x^*$x∗ is
$x^*\approx \hat{x}-\displaystyle\frac{f(\hat{x})}{f'(\hat{x})}.$x∗≈x^−f′(x^)f(x^)​.

2.3.2 Algorithm

Definition

• Starts with an initial approximation $x_0$x0​
• Defined iteration scheme by
  $x_n=x_{n-1}-\displaystyle\frac{f(x_{n-1})}{f'(x_{n-1})},\forall n\geq 1$xn​=xn−1​−f′(xn−1​)f(xn−1​)​,∀n≥1
• This scheme generates the sequence $\{x_n\}_0^\infty${xn​}0∞​

Example

Pseudo-Code

The function $f$f is differentiable and $p_0$p0​ is an initial approximation.

• INPUT: Initial approximation $p_0$p0​, tolerance $TOL$TOL, Maximum number of iteration $N$N.

• OUTPUT: approximate solution $p$p or message of failure.

• Step $1$1: Set $n = 1$n=1.

• Step $2$2: While $n\leq N$n≤N, do Steps $3～5$3～5.

  • Step $3$3: Set $p= p_0-\displaystyle\frac{f(p_0)}{f'(p_0)}$p=p0​−f′(p0​)f(p0​)​.
  • Step $4$4: If $|p-p_0|<TOL$∣p−p0​∣<TOL, then output $p$p; (Procedure complete successfully.) Stop!
  • Step $5$5: Set $n=n+1, p_0=p$n=n+1,p0​=p.

• Step $6$6: OUTPUT “Method failed after $N$N iterations.” STOP!

2.3.3 Convergence

Theorem

(1) $f\in C^2[a,b].$f∈C2[a,b].

(2) $p\in [a,b]$p∈[a,b] is such that $f(p)=0$f(p)=0 and $f'(p)\not=0$f′(p)​=0.

Then there exists a $\delta>0$δ>0 such that Newton’s method generates a sequence $\{p_n\}_1^\infty${pn​}1∞​ converging to $p$p for any initial approximation $p_0\in[p-\delta,p+\delta]$p0​∈[p−δ,p+δ].

Proof

$p_n=g(p_{n-1})\\ g(x)=x-\displaystyle\frac{f(x)}{f'(x)}$pn​=g(pn−1​)g(x)=x−f′(x)f(x)​

Things need to prove: According to the proofing process of the fixed point method, we need to find an interval $[p-\delta,p+\delta]$[p−δ,p+δ] that $g$g maps into itself, and $|g'(x)|\leq k<1$∣g′(x)∣≤k<1 for all $x\in[p-\delta,p+\delta].$x∈[p−δ,p+δ]. (Existence and Uniqueness)

Proving Process:

1. $\exists \delta_1>0,g(x) \in C[p-\delta_1,p+\delta_1]$∃δ1​>0,g(x)∈C[p−δ1​,p+δ1​]

   Since $f'(p)\not=0$f′(p)​=0 and $f'$f′ is continuous, there exists $\delta_1>0$δ1​>0 such that $f'(x)\not= 0$f′(x)​=0 and $f’(x)\in C[a,b] $ with $\forall x\in [p-\delta_1,p+\delta_1]$∀x∈[p−δ1​,p+δ1​].

2. $g'(x) \in C[p-\delta_1,p+\delta_1]$g′(x)∈C[p−δ1​,p+δ1​]

   $g'(x)=\displaystyle\frac{f(x)f''(x)}{(f'(x))^2}$g′(x)=(f′(x))2f(x)f′′(x)​ for all $x\in [p-\delta_1,p+\delta_1]$x∈[p−δ1​,p+δ1​]. Since $f\in C^2[a,b]$f∈C2[a,b], $g'(x)\in [p-\delta_1,p+\delta_2]$g′(x)∈[p−δ1​,p+δ2​].

3. $|g'(x)|\leq k< 1$∣g′(x)∣≤k<1

   $f(p)=0,g'(p)=0$f(p)=0,g′(p)=0

   Since $g'\in C[p-\delta_1,p+\delta_1]$g′∈C[p−δ1​,p+δ1​], there exists a $\delta$δ with $0<\delta<\delta_1$0<δ<δ1​ and

   $|g'(x)|\leq k<1, \forall x\in [p-\delta,p+\delta].$∣g′(x)∣≤k<1,∀x∈[p−δ,p+δ].

4. $g\in[p-\delta,p+\delta]\mapsto [p-\delta,p+\delta]$g∈[p−δ,p+δ]↦[p−δ,p+δ]

   According to the Mean Value Theorem, if $x\in[p-\delta,p+\delta]$x∈[p−δ,p+δ], there exists
   $\xi \in [x,p],|g(x)-g(p)|=|g'(\xi)||x-p|\\ |g(x)-p|=|g(x)-g(p)|=|g'(\xi)||x-p|\leq k|x-p|<|x-p|<\delta$ξ∈[x,p],∣g(x)−g(p)∣=∣g′(ξ)∣∣x−p∣∣g(x)−p∣=∣g(x)−g(p)∣=∣g′(ξ)∣∣x−p∣≤k∣x−p∣<∣x−p∣<δ
   Thus, $g\in[p-\delta,p+\delta]\mapsto [p-\delta,p+\delta]$g∈[p−δ,p+δ]↦[p−δ,p+δ].

According to the proving process above, all the hypotheses of the Fixed-Point Theorem are satisfied for $g(x)=x-\displaystyle\frac{f(x)}{f'(x)}$g(x)=x−f′(x)f(x)​. Therefore, the sequence ${p_n}_{n=1}^\infty $ defined by
$p_n=g(p_n-1),\forall n\geq 1$pn​=g(pn​−1),∀n≥1
converges to $p$p for any $p_0\in[p-\delta,p+\delta].$p0​∈[p−δ,p+δ].

2.3.4 Secant Method

Background

For Newton’s method, each iteration requires evaluation of both function $f(x_k)$f(xk​) and its derivative $f'(x_k)$f′(xk​), which may be inconvenient or expensive.

Improvement

Derivative is approximated by finite difference using two successive iterates, so iteration becomes
$x_{k+1}=x_k-f(x_k)*\displaystyle\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}.$xk+1​=xk​−f(xk​)∗f(xk​)−f(xk−1​)xk​−xk−1​​.
This method is known as secant method.

Example

Pseudo-Code

• INPUT: Initial approximation $p_0,p_1$p0​,p1​, tolerance $TOL$TOL, Maximum number of iteration $N$N.

• OUTPUT: approximate solution $p$p or message of failure.

• Step $1$1: Set $n = 2,q_0=f(p_0),q_1=f(p_1)$n=2,q0​=f(p0​),q1​=f(p1​).

• Step $2$2: While $n\leq N$n≤N, do Steps $3～5$3～5.

  • Step $3$3: Set $p= p_1-q_1*\displaystyle\frac{p_1-p_0}{q_1-q_0}$p=p1​−q1​∗q1​−q0​p1​−p0​​.(Compute $p_i$pi​)
  • Step $4$4: If $|p-p_1|<TOL$∣p−p1​∣<TOL, then output $p$p; (Procedure complete successfully.) Stop!
  • Step $5$5: Set $n=n+1, p_0=p_1,p_1=p, q_0=q_1,q_1=f(p)$n=n+1,p0​=p1​,p1​=p,q0​=q1​,q1​=f(p).

• Step $6$6: OUTPUT “Method failed after $N$N iterations.” STOP!

2.3.5 False Position Method

Definition

To find a solution of $f(x)=0$f(x)=0 for a given continuous function $f$f on the interval $[p_0,p_1]$[p0​,p1​], where $f(p_0)$f(p0​) and $f(p_1)$f(p1​) have opposite signs
$f(p_0)*f(p_1)<0.$f(p0​)∗f(p1​)<0.

The approximation $p_2$p2​ is chosen in same manner as in Secant Method, as the $x$x-intercept ($x$x轴截距) of the line joining $(p_0,f(p_0))$(p0​,f(p0​)) and $(p_1,f(p_1))$(p1​,f(p1​)).

To decide the Secant Line to compute $p_3$p3​, we need to check $f(p_2)*f(p_1)$f(p2​)∗f(p1​) or $f(p_2)*f(p_0)$f(p2​)∗f(p0​).

If $f(p_2)*f(p_1)$f(p2​)∗f(p1​) is negative, we choose $p_3$p3​ as the $x$x-intercept for the line joining $(p_1,f(p_1)$(p1​,f(p1​) and $p_2,f(p_2)$p2​,f(p2​).

In a similar manner, we can get a sequence $\{p_n\}_2^\infty${pn​}2∞​ which approximates to the root.

Example

Pseudo-Code

• INPUT: Initial approximation $p_0,p_1$p0​,p1​, tolerance $TOL$TOL, Maximum number of iteration $N$N.

• OUTPUT: approximate solution $p$p or message of failure.

• Step $1$1: Set $n = 2,q_0=f(p_0),q_1=f(p_1)$n=2,q0​=f(p0​),q1​=f(p1​).

• Step $2$2: While $n\leq N$n≤N, do Steps $3～6$3～6.

  • Step $3$3: Set $p= p_1-q_1*\displaystyle\frac{p_1-p_0}{q_1-q_0}$p=p1​−q1​∗q1​−q0​p1​−p0​​.(Compute $p_i$pi​)
  • Step $4$4: If $|p-p_1|<TOL$∣p−p1​∣<TOL, then output $p$p; (Procedure complete successfully.) Stop!
  • Step $5$5: Set $n=n+1, q=f(p)$n=n+1,q=f(p).
  • Step $6$6: If $q*q_1<0$q∗q1​<0, set $p_0=p, q_0=q$p0​=p,q0​=q; else set $p_1=p, q_1=q.$p1​=p,q1​=q.

• Step $6$6: OUTPUT “Method failed after $N$N iterations.” STOP!

2.4 Error Analysis for Iteration Methods

2.4.1 The Rate of Sequence Convergence

Definition

Suppose $\{p_n\}_{n=0}^\infty${pn​}n=0∞​ is a sequence that converges to $p$p, with $p_n\not=p$pn​​=p for all n.

If positive constants $\lambda$λ and $\alpha$α exist with
$\lim_{n\rightarrow\infty} \displaystyle\frac{|p_{n+1}-p|}{|p_n-p|^\alpha}=\lambda,$n→∞lim​∣pn​−p∣α∣pn+1​−p∣​=λ,

then $\{p_n\}_{n=0}^\infty${pn​}n=0∞​ converges to $p$p of order $\alpha$α, with asymptotic error constant $\lambda$λ.

Properties

1. A sequence with a high order of convergence converges more rapidly than a sequence with a lower order.
2. The asymptotic constant affects the speed of convergence but is not as important as the order.

Example

1. If $\alpha=1$α=1, the sequence is linearly convergent.
2. If $\alpha =2$α=2, the sequence is quadratically convergent.

Summary

Using the Mean Value Theorem to prove Linear Convergence and the Taylor’s Theorem to prove Quadratic Convergence with $g'(p)=0.$g′(p)=0.

2.4.2 Convergent Order of Fixed-Point Iteration (Improved)

Convergent Oder of Fixed-Point Iteration

(1) $g\in C[a,b]$g∈C[a,b] for all $x\in[a,b]$x∈[a,b]
(2) $g'(x)$g′(x) is continuous on $(a,b)$(a,b) and a positive constant $0<k<1$0<k<1 exists with $|g'(x)|\leq k$∣g′(x)∣≤k, for all $x\in(a,b)$x∈(a,b).

If $g'(p)\not=0$g′(p)​=0, then for any number $p_0$p0​ in $[a,b]$[a,b], the sequence $p_n=g(p_{n-1})$pn​=g(pn−1​), for $n\geq 1$n≥1, converges only linearly to the unique fixed point $p$p in $[a,b]$[a,b].

Proof

$p_{n+1}-p=g(p_n)-g(p)=g'(\xi_n)(p_n-p),$pn+1​−p=g(pn​)−g(p)=g′(ξn​)(pn​−p),
where $\xi_n$ξn​ is between $p_n$pn​ and $p$p.

Since $\{p_n\}_{n=0}^\infty${pn​}n=0∞​ converges to $p$p, and $\xi_n$ξn​ is between $p_n$pn​ and $p$p, thus $\{\xi_n\}_{n=0}^\infty${ξn​}n=0∞​ also converges to $p$p.

Thus,
$\lim_{n\rightarrow\infty}\displaystyle\frac{|p_{n+1}-p|}{|p_n-p|}=\lim_{n\rightarrow\infty}|g'(\xi_n)|=|g'(p)|,$n→∞lim​∣pn​−p∣∣pn+1​−p∣​=n→∞lim​∣g′(ξn​)∣=∣g′(p)∣,
fixed-point iteration exhibits linear convergence with asymptotic error constant $|g'(p)|$∣g′(p)∣ whenever $g'(p)\not=0$g′(p)​=0, which also implies that higher-order convergence for fixed-point methods can occur only when $g'(p)=0$g′(p)=0.

Quadratical Convergence

Let $p$p be a solution of the equation $x=g(x)$x=g(x).

(1) $g'(p)=0$g′(p)=0

(2) $g''$g′′ is continuous and strictly bounded by $M$M on an open interval $I$I containing $p$p.

Then there exists a $\delta >0$δ>0 such that, for $p_0\in [p-\delta, p+\delta]$p0​∈[p−δ,p+δ], the sequence defined by $p_n=g(p_{n-1})$pn​=g(pn−1​), when $n\geq 1$n≥1, converges at least quadratically to $p$p.

Moreover, for sufficiently large values of $n$n,
$|p_{n+1}-p|<\displaystyle\frac{M}{2}|p_n-p|^2.$∣pn+1​−p∣<2M​∣pn​−p∣2.

Proof

Due to the two conditions described above,
$g(x)=g(p)+g'(p)(x-p)+\displaystyle\frac{g''(\xi)}{2}(x-p)^2=p+\displaystyle\frac{g''(\xi)}{2}*(x-p)^2$g(x)=g(p)+g′(p)(x−p)+2g′′(ξ)​(x−p)2=p+2g′′(ξ)​∗(x−p)2
is derived, that $\xi$ξ lies between $x$x and $p$p.

Thus,
$p_{n+1}=g(p_n)=p+\displaystyle\frac{g''(\xi)}{2}*(p_n-p)^2\\ p_{n+1}-p=\displaystyle\frac{g''(\xi)}{2}*(p_n-p)^2\\ \lim_{n\rightarrow\infty}g''(\xi)=g''(p)\\ \lim_{n\rightarrow\infty}\displaystyle\frac{|p_{n+1}-p_n|}{|p_n-p|^2}=\lim_{n\rightarrow\infty}\displaystyle\frac{g''(\xi)}{2}=\displaystyle\frac{g''(p)}{2}$pn+1​=g(pn​)=p+2g′′(ξ)​∗(pn​−p)2pn+1​−p=2g′′(ξ)​∗(pn​−p)2n→∞lim​g′′(ξ)=g′′(p)n→∞lim​∣pn​−p∣2∣pn+1​−pn​∣​=n→∞lim​2g′′(ξ)​=2g′′(p)​

Since $g''$g′′ is strictly bounded by $M$M on the interval $|p-\delta,p+\delta|$∣p−δ,p+δ∣, for sufficiently large values of $n$n,
$|p_{n+1}-p|<\displaystyle\frac{M}{2}|p_n-p|^2$∣pn+1​−p∣<2M​∣pn​−p∣2
is also derived.

Construct a quadratically convergent fixed-point problem

Let
$g(x)=x-\phi(x)f(x)\\ g'(x)=1-\phi'(x)f(x)-\phi(x)f'(x)$g(x)=x−ϕ(x)f(x)g′(x)=1−ϕ′(x)f(x)−ϕ(x)f′(x)

And the condition is $g'(p)=0$g′(p)=0, thus $\phi(p)=\displaystyle\frac{1}{f'(p)}$ϕ(p)=f′(p)1​.

A reasonable approach is to let $\phi(x)=\displaystyle\frac{1}{f'(x)}$ϕ(x)=f′(x)1​, which is the Newton’s method.

Remarks

1. the convergence rate of Fixed-Point iteration is usually linear, with constant $C=|g'(p)|$C=∣g′(p)∣.
2. But if $g'(p)=0$g′(p)=0, then the convergence rate is at least $quadratic$quadratic.

2.4.3 Zero of Multiplicity

Definition

A solution $p$p of $f(x)=0$f(x)=0 is a zero of multiplicity $m$m of $f(x)$f(x) if for $x\not=p$x​=p, we can write
$f(x)=(x-p)^mq(x),$f(x)=(x−p)mq(x),
where
$\lim_{x\rightarrow p}q(x)\not=0.$x→plim​q(x)​=0.

Theorem

1. $f\in C^1[a,b]$f∈C1[a,b] has a simple zero at $p$p in $(a,b)$(a,b) if and only if $f(p)=0$f(p)=0, but $f'(p)\not=0$f′(p)​=0.
2. The function $f\in C^m[a,b]$f∈Cm[a,b] has a zero of multiplicity $m$m at $p$p if and only if
   $0=f(p)=f'(p)=f''(p)=...=f^{(m-1)}(p).$0=f(p)=f′(p)=f′′(p)=...=f(m−1)(p).
   but $f^m(p)\not=0$fm(p)​=0.

Proof

If $f$f has a simple zero at $p$p, then
$f(p)=0\\ f(x)=(x-p)*q(x)\\ \lim_{x\rightarrow p}q(x)\not=0.$f(p)=0f(x)=(x−p)∗q(x)x→plim​q(x)​=0.

Since $f\in C^1[a,b]$f∈C1[a,b],
$f'(p)=\lim_{x\rightarrow p}f'(x)=\lim_{x\rightarrow p}[q(x)+(x-p)*q'(x)]=\lim_{x\rightarrow p}q(x)\not=0.$f′(p)=x→plim​f′(x)=x→plim​[q(x)+(x−p)∗q′(x)]=x→plim​q(x)​=0.

2.4.4 Convergence of Newton’s Method

Property

Newton’s method transforms nonlinear equation $f(x)=0$f(x)=0 into fixed-point problem $x=g(x)$x=g(x) with $g(x)=x-\displaystyle\frac{f(x)}{f'(x)}$g(x)=x−f′(x)f(x)​.

1. If $p$p is a simple root, $f(p)=0,f'(p)\not=0,g'(p)=0$f(p)=0,f′(p)​=0,g′(p)=0, thus the convergence rate is quadratic. (Iterations must start close enough to root.)

2. If $p$p is a root of multiplicity,
   $f(x)=(x-p)^mq(x)\\ g'(p)=1-\displaystyle\frac{1}{m}\not=0,$f(x)=(x−p)mq(x)g′(p)=1−m1​​=0,
   thus the convergence rate is linear.

the Method of avoiding multiple root

$f(x)=(x-p)^mq(x)\\ u(x)=\displaystyle\frac{f(x)}{f'(x)}\\ u(x)=(x-p)*\displaystyle\frac{q(x)}{mq(x)+(x-p)q'(x)}\\ u'(x)=\displaystyle\frac{1}{m}\not=0.$f(x)=(x−p)mq(x)u(x)=f′(x)f(x)​u(x)=(x−p)∗mq(x)+(x−p)q′(x)q(x)​u′(x)=m1​​=0.

Thus $p$p is a simple root of $u(x)$u(x). Then we change the Newton’s method into
$g(x)=x-\displaystyle\frac{u(x)}{u'(x)}=x-\displaystyle\frac{f(x)f'(x)}{f'(x)^2-f(x)f''(x)}$g(x)=x−u′(x)u(x)​=x−f′(x)2−f(x)f′′(x)f(x)f′(x)​
whose convergence rate is also quadratic.

2.4.5 Convergence rate of Secant Method

1. Convergence rate of secant method is normally superlinear, with $r\approx1.618$r≈1.618, which is lower than Newton’s method.
2. Secant method need to evaluate two previous functions per iteration, there is no requirement to evaluate the derivative.
3. Its disadvantage is that it needs two starting guesses which close enough to the solution in order to converge.

2.5 Accelerating Convergence

2.5.1 Aitken’s method

Background

Accelerating the convergence of a sequence that is linearly convergent, regardless of its origin or application.

$\lim_{n\rightarrow \infty} \displaystyle\frac{p_{n+1}-p}{p_n-p}=\lambda,\lambda\not=0.$n→∞lim​pn​−ppn+1​−p​=λ,λ​=0.
Thus, when $n$n is sufficiently large,
$\displaystyle\frac{p_{n+1}-p}{p_n-p}\approx \displaystyle\frac{p_{n+2}-p}{p_n+1-p}\\ p\approx\displaystyle\frac{p_n*p_{n+2}-p_{n+1}^2}{p_{n+2}-2*p_{n+1}+p_n}\\ p\approx p_n-\displaystyle\frac{(p_{n+1}-p_{n})^2}{p_{n+2}-2*p_{n+1}+p_n}$pn​−ppn+1​−p​≈pn​+1−ppn+2​−p​p≈pn+2​−2∗pn+1​+pn​pn​∗pn+2​−pn+12​​p≈pn​−pn+2​−2∗pn+1​+pn​(pn+1​−pn​)2​

Aitken’s $\Delta$Δ method is to define a new sequence ${\hat{p}}_{n=0}^\infty:$p^​n=0∞​:
$\hat{p}=p_n-\displaystyle\frac{(p_{n+1}-p_{n})^2}{p_{n+2}-2*p_{n+1}+p_n},$p^​=pn​−pn+2​−2∗pn+1​+pn​(pn+1​−pn​)2​,
whose convergence rate is faster than the original sequence $\{p_n\}_{n=0}^\infty${pn​}n=0∞​.

Definition

Given the sequence $\{p_n\}_{n=0}^\infty${pn​}n=0∞​, the forward difference $\Delta p_n$Δpn​ is defined by
$\Delta p_n=p_{n+1}-p_n,n\geq 0.$Δpn​=pn+1​−pn​,n≥0.
Higher powers $\Delta^kp_n$Δkpn​ are defined recursively by
$\Delta^k p_n=\Delta(\Delta^{k-1}p_{n}),k\geq 2.$Δkpn​=Δ(Δk−1pn​),k≥2.
For example,
$\Delta^2p_n=\Delta(\Delta p_{n})=\Delta(p_{n+1}-p_n )=\Delta(p_{n+1})-\Delta(p_n)=p_{n+2}-2p_{n+1}+p_n.$Δ2pn​=Δ(Δpn​)=Δ(pn+1​−pn​)=Δ(pn+1​)−Δ(pn​)=pn+2​−2pn+1​+pn​.
Thus, $\hat{p_n}=p_n-\displaystyle\frac{(\Delta p_n)^2}{\Delta^2p_n}$pn​^​=pn​−Δ2pn​(Δpn​)2​.

Theorem

Condition:
$\lim_{n\rightarrow \infty} \displaystyle\frac{p_{n+1}-p}{p_n-p}=\lambda,\lambda\not=0.\\ (p_n-p)(p_{n+1}-p)>0$n→∞lim​pn​−ppn+1​−p​=λ,λ​=0.(pn​−p)(pn+1​−p)>0

Result: the sequence $\{\hat{p_n}\}_{n=0}^\infty${pn​^​}n=0∞​ converges to $p$p faster than $\{p_n\}_{n=0}^\infty${pn​}n=0∞​ in the sense that
$\lim_{n\rightarrow\infty}\displaystyle\frac{\hat{p}_n-p}{p_n-p}=0.$n→∞lim​pn​−pp^​n​−p​=0.

Proof:

2.5.2 Steffensen’s method

Definition

The function is $p=g(p)$p=g(p), and the initial approximation is $p_0$p0​, $\hat{p_0}=p_0-\displaystyle\frac{(\Delta p_0)^2}{\Delta^2p_0}$p0​^​=p0​−Δ2p0​(Δp0​)2​.

Assume that $\hat{p_0}$p0​^​ is a better approximation than $p_2$p2​, so applying fixed-point iteration to $\hat{p_0}$p0​^​ instead of $p_2$p2​, and the computing process shows below.

Theorem

Suppose that $x=g(x)$x=g(x) has the solution $p$p with $g'(p)\not=1$g′(p)​=1.

If there exists a $\delta>0$δ>0 such that $g\in C^3[p-\delta,p+\delta]$g∈C3[p−δ,p+δ],

then Steffensen’s method gives quadratic convergence for any $p_0\in [p-\delta,p+\delta].$p0​∈[p−δ,p+δ].

Pseudo-Code

• INPUT: Initial approximation $p_0$p0​, tolerance $TOL$TOL, Maximum number of iteration $N$N.

• OUTPUT: approximate solution $p$p or message of failure.

• Step $1$1: Set $n = 1$n=1.

• Step $2$2: While $n\leq N$n≤N, do Steps $3～5$3～5.

  • Step $3$3: Set $p_1=g(p_0),p_2=g(p_1),p= p_0-\displaystyle\frac{(p_1-p_0)^2}{p_2-2p_1+p_0}$p1​=g(p0​),p2​=g(p1​),p=p0​−p2​−2p1​+p0​(p1​−p0​)2​.
  • Step $4$4: If $|p-p_0|<TOL$∣p−p0​∣<TOL, then output $p$p; (Procedure complete successfully.) Stop!
  • Step $5$5: Set $n=n+1, p_0=p$n=n+1,p0​=p.

• Step $6$6: OUTPUT “Method failed after $N$N iterations.” STOP!

2.6 Zeros of Polynomials and Muller’s Method

2.6.1 Polynomial Theorem

2.6.2 Horner’s Method

Background

A more efficient method to calculate the $P(x_0)$P(x0​) and $P'(x_0)$P′(x0​) for a polynomial $P(x)$P(x).

Theorem

Let
$P(x)=\sum\limits_{i=0}^{i=n}a_ix^i.$P(x)=i=0∑i=n​ai​xi.

• Construction process for P(x_0) (Substitute formulas one by one to verify)

if $b_n=a_n$bn​=an​ and
$b_k=a_k+b_{k+1}x_0,k\in [0,n-1],$bk​=ak​+bk+1​x0​,k∈[0,n−1],
then $b_0=P(x_0)$b0​=P(x0​).

Moreover, if
$Q(x)=\sum\limits_{i=1}^{n}b_ix^{i-1}$Q(x)=i=1∑n​bi​xi−1
then
$P(x)=(x-x_0)Q(x)+b_0.$P(x)=(x−x0​)Q(x)+b0​.

• Construction process for P’(x_0) (Substitute formulas one by one to verify)

$P(x)=(x-x_0)Q(x)+b_0\\ P'(x)=Q(x)+(x-x_0)Q'(x)\\ P'(x_0)=Q(x_0)$P(x)=(x−x0​)Q(x)+b0​P′(x)=Q(x)+(x−x0​)Q′(x)P′(x0​)=Q(x0​)

Let $Q(x)=\sum\limits_{i=1}^{n}b_ix^{i-1}=(x-x_0)R(x)+c_1$Q(x)=i=1∑n​bi​xi−1=(x−x0​)R(x)+c1​, where $R(x)=\sum\limits_{i=2}^{n}c_ix^{i-2}$R(x)=i=2∑n​ci​xi−2. Thus
$c_n=b_n,\\ c_k=b_k+c_{k+1}x_0,k\in[1,n-1],\\ Q(x_0)=c_1=P'(x_0).$cn​=bn​,ck​=bk​+ck+1​x0​,k∈[1,n−1],Q(x0​)=c1​=P′(x0​).

Pseudo-Code

To compute the value $P(x_0)$P(x0​) and $P'(x_0)$P′(x0​) for a function $P(x)=\sum\limits_{i=0}^{n}a_ix^i.$P(x)=i=0∑n​ai​xi.

• INPUT: Degree $n$n, coefficients $a_0,a_1,...,a_n$a0​,a1​,...,an​ of polynomial $P(x)$P(x), point $x_0$x0​.

• OUTPUT: Values of $P(x_0)$P(x0​) and $P'(x_0)$P′(x0​).

• Step $1$1: Set $y = a_n$y=an​ ($b_n$bn​ for $Q$Q), $z=0$z=0 ($c_{n+1}$cn+1​ for $R$R).

• Step $2$2: For $j=n-1,n-2,...,0$j=n−1,n−2,...,0, set

  • $z=y+z*x_0$z=y+z∗x0​ ($c_{j+1}$cj+1​ for $R$R),
  • $y =a_j+y*x_0$y=aj​+y∗x0​ ($b_j$bj​ for $Q$Q).

• Step $3$3: OUTPUT $y:P(x_0)$y:P(x0​) and $z:P'(x_0)$z:P′(x0​).

2.6.3 Deflation Method

Newton’s method combined with Honor’s method

Deflation Method (压缩技术)

Suppose that $x_N$xN​ in the Nth iteration of the Newton-Raphson procedure, is an approximation zero of $P(x)$P(x), then
$P(x)=(x-x_N)Q(x)+b_0=(x-x_N)Q(x)+P(x_N)\approx (x-x_N)Q(x).$P(x)=(x−xN​)Q(x)+b0​=(x−xN​)Q(x)+P(xN​)≈(x−xN​)Q(x).

Let $\hat{x_1}=x_N$x1​^​=xN​ be the approximate zero of $P$P, and $Q_1(x)=Q(x)$Q1​(x)=Q(x) be the approximate factor, then we have
$P(x)\approx (x-\hat{x_1})Q_1(x).$P(x)≈(x−x1​^​)Q1​(x).

To find the second approximate zero of $P(x)$P(x), we can use the same procedure to $Q_1(x)$Q1​(x), give $Q_1(x)\approx(x-\hat{x_2})Q_2(x)$Q1​(x)≈(x−x2​^​)Q2​(x), where $Q_2(x)$Q2​(x) is a polynomial of degree $n-2$n−2. Thus $P(x)\approx (x-\hat{x_1})Q_1(x)\approx (x-\hat{x_1})(x-\hat{x_2})Q_2(x)$P(x)≈(x−x1​^​)Q1​(x)≈(x−x1​^​)(x−x2​^​)Q2​(x).

Repeat this procedure, till $Q_{n-2}(x)$Qn−2​(x) which is an quadratic polynomial and can be solved by quadratic formula. We can get all approximate zeros of $P(x)$P(x). This method is called deflation method.

2.6.4 Muller’s Algorithm