численные методы

нахождение корня функции f : f (x) = 0

метод бисекции

требует априорного знания значений функции на границах промежутка [a,b], в котором ищется корень. предполагается, что этот промежуток содержит корень:

    f(a) ⋅ f(b) < 0

метод использует итерации для последовательного разделения промежутка на половины, одна из которых содержит корень. отбрасывается та часть [a,b]/2 на которой значение функции не меняет знака. процесс останавливается при достижении заданной точности приближения

Bisection Method

псевдокод:

алгоритм задает последовательность cₙ , которая сходится к значению корня

    lim cₙ = r
   n → ∞ 

listing Haskell

метод фиксированной точки

для того, чтобы найти корни функции f(x) можно использовать итерацию фиксированной точки

фиксированная точка x для данной функции g(x) есть точка такая, что g(x) = x

имея монотонную функцию g(x) : x ∈ [a,b] & g(x) ∈ [a,b] (т.е. функция отображает свою область определения в нее же

    g : [a,b] → [a,b]  

доказательство:
положим, что мы имеем две фиксированные точки p и q функции g(x) : g(p) = p , g(q) = q
по теореме о среднем значении производной

    |p - q| = |g(p) - g(q)| = |g'(χ) ⋅ (p - q)| ≤ |g'(χ)| ⋅ |p - q| ≤ k ⋅ |p - q|

    где χ ∈ [p,q] ,  0 < k < 1 

    |p - q| ≤ k ⋅ |p - q| < |p - q| , 0 < k < 1 

как найти фиксированную точку? необходимо выбрать начальную точку p₀

    p₁=g(p₀),
    p₂=g(p₁)
    ..
    pn+1=g(pn) 

    пусть p=g(p)
    p₀, p₁, p₂ .. pn есть последовательность предположений

    |pn - p| = |g(pn) - g(p)| ⇒ |g'(χ) * (pn - p)| = |g'(χ)| * |pn - p| ≤ k * |pn - p|

    ( g(pn) - g(p) ) / ( pn - p ) = g'(χ)  
    
    где χ лежит между pn and p

итерация сгенерит последовательность p_n , n = 0..∞ такую что

  limn → ∞ g(pn) = g(p) = p

например, найдем корень для функции f(x) = x^3 - 2*x - 5

разрешим уравнение f(x) относительно x в точке, где f(x) = 0 :

    0 = x^3 - 2*x - 5
    x^3 = 2 * x + 5
    x = (2 * x + 5) ^ (1/3)
  

listing Haskell

запустим в repl:

λ> main
2
2.079
2.091
2.092
2.093
λ> f 2.093
-1.73016429999997e-2

метод Ньютона

если есть функция f(x) чья производная f'(x) непрерывна возле значения корня r, тогда для нахождения корня может быть использован метод Ньютона, который сходится быстрее чем метод биекций. он тоже требует первоначального предположения о величине r (initial guess)

метод использует уравнение для касательной в первоначальной точке:

f'(x₀) = f(x₀) / (x - x₀)

полагая, что f(x) определена на [a,b] и ∃ r ∈ [a,b] : f(r) = 0, метод Ньютона определяет следующую последовательность итераций:

    x_n+1 = x_n - f(x_n) / f'(x_n) для n ≥ 0

точка x_n+1 это точка пересечения касательной и оси абсцисс, которая для следующей итерации берется уже в качестве первоначальной. этот процесс повторяется, пока не будет достигнута желаемая точность в определении корня

попробуем найти ⎷2 : listing Haskell

запустим в repl:

λ> main
2.0000000000000  
1.5000000000000
1.4166666666667
1.4142156862745
1.4142135623747
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
1.4142135623731
  

"я Катя, я дура, но я вот нашла корень из двух"

полином Лагранжа и матрица Вандермонда

          n
  Pn(x) = Σ ak * xk  
          0

  f(x_k) = Pn(x_k), k=0,1,2,..,n

получаем систему линейных уравнений для коэффициентов a0, a1, a2, .., an:

    f(x₀)           1    x₀    x₀2    ...    x₀n        a0 
    f(x₁)           1    x₁    x₁2    ...    x₁n        a1
    ...            ...
    f(xn)           1    xn    xn2    ...    xnn        an 

матрица называется матрицей Вандермонда: octave listing

можно увеличить число точек для того, чтобы результат стал точнее, но для больших значений N численное решение уравнений становится затратным и приводит к проблемам со стабильностью. если вы используете octave/matlab, то очень скоро начнете получать ворнинги :

    Warning: Matrix is close to singular or badly scaled.
    Results may be inaccurate. RCOND = 9.377592e-48.

для функции f(x) = cos(x) на [0 , 2 * pi] (красный график) ее апроксимация (синий график) будут выглядеть как-то так:

ошибка интерполяции показана на правом графике. по горизонтальной оси откладывется количество точек апроксимации. видно, что с увеличением точек ошибка резко увеличивается. а вы ожидали, что увеличивая количество точек, вы ее уменьшаете

графики построены для 50 точек - если вы увеличите количество точек, то ошибка будет более заметной

поэтому требуется какой-то другой путь для полиномиальной интерполяции

Say you’re trying to find an interpolation polynomial that agrees with the function f(x) at two points, a P_1(x). The function can then be approximated by taking linear combination of two polynomial of degree 1:

L_{1,0}=\frac{x₁ -x}{x₁ -x₀} L_{1,0}=\frac{x -x₀}{x₁ -x₀}$$

Hence, whenever L_{1,0}(x₀)=1 then L_{1,1}(x₁)=0, and whenever L_{1,1}(x₀)=0 then L_{1,1}(x₁)=1.This way, the following equation will creates secant lines connecting to each other: P_1(x)=L_{1,0}(x)+L_{1,1}f(x₁).

The above equation can be extended for n+1 points as the following:

The index n denotes the degree of the polynomial and k indicates f(x_k) being the coefficient in P_n(x).

The Lagrange polynomials satisfies:

L_{n,k}(x_j) = \begin{cases} 1 & j = k\\ 0 & j \neq k \end{cases}

The same thing happens with Lagrange Interpolation that happened in Vandermonde Matrix. This problem, or misbehavoir, is known as Runge’s phenomenom.

Runge’s Phenomenon

Runge’s Phenomenon:

Its a problem that occurs with large data during interpolation of polynomial of higher degrees. The more data points used, the higher the polynomial degree. To investigate this phenomena, I will use f(x)=\frac{1}{1+25x^2} on [-1,1].

Equally Spaced Nodes:

error plot of f(x) with 10 interp point

f(x) with 10 interpolation points

The error plot looks great. Anyone would be happy to get an error like that, almost machine precision(10e-16). However, the reconstruction of the function(in blue) in poor. With more points, we see the runge’s phenomena, the oscillation toward the end:

The error plot of this interpolation shows the impact of Runge’s Phenomena:

The error toward the end point, where the oscillations occur in the interpolation, is very poor.

Chebyshev Points:

x = \cos \frac{k\pi}{n}\text{ for }0 \le k \le n

Chebyshev point is better in polynomial interpolation because it reduces the ocilation as n increases. Below is a plot of f(x) with just 16 points.

Interpolation Using Splines

This method has advantages over the use of high degree polynomial, which have tendency to oscillate between data values. This method is a series of third degree polynomials joining the points together. Only four unknown coefficients need to computed to construct the cubic spline.

s_0 (x) = a_0 + b_{0}x + c_0 x^2 + d_0 x^3 \\ s_1 (x) = a_1 + b_{1}x + c_1 x^2 + d_1 x^3 \\ s_2 (x) = a_2 + b_{2}x + c_2 x^2 + d_2 x^3 \\ \vdots \\ s_n (x) = a_n + b_n x + c_n x^2 + d_nx^3$$

Each of these cubic polynomial must pass through the two points it joins.

Smoothness conditions is most important in the spline method. Therefore, continuity of the slope f'(x) and the second derivative f''(x), which determines the concavity of the function f(x), must also be agree between the adjacent data points. These conditions controls the oscillations that usually happens in the high order polynomials from happening in Cubic Spline.

Properties of Cubic Spline:

1. S(x)=s_k(x)=s_{k,0} + s_{k,1}(x-x_k) + s_{k,2}(x-x_k)^2 + s_{k,3}(x-x_k)^3 for k=0,1,2, \cdots n
2. s(x_k)= f(x_k) for k=0,1,2, \cdots n. The splines passes through each data points.
3. s_k(x_{k+1})=s_{k+1}(x_{k+1}) for k=0,1,2, \cdots n. The spline forms a continuous functions on the interval.
4. s'_k(x_{k+1})=s'_{k+1}(x_{k+1}) for k=0,1,2,\cdots n. The spline forms a smooth function.
   
5. s''_k(x_{k+1})=s''_{k+1}(x_{k+1}) for k=0,1,2,\cdots n. The second derivative is continuous.

The properties of cubic splines makes sure that Runge phenomenon does not occur in this interpolation method.

Matlab Code:

#!/usr/local/bin/octave -q

close all; clc; clear all

x = linspace(-1,1,15);
x=-cos(pi*(0:20)/20); % the quadrature points [-1,1]
y = 1./(1+25*x.^2);

xx = linspace(-1,1,100);
yy= 1./(1+25*xx.^2);
y = spline(x,y,xx);

plot(xx,yy,’b’)
hold on
plot(xx,y,’r’)

pause ()

The result of this code is the following:

10 equally spaced points

10 chebyshev points

15 points

Quadrature

Quadrature is the numerical evaluation of definit inegrals in the following form:

                    b
             F(x) = ∫ f(x) dx
                    a

This can approximated by in this form:

            F(x) ≅ Σ  w_k f(x_k) 
                   k

There are many forms of quadratures available. Each of these methods are in the same form above. The only difference is the choice of x_k called the quadrature nodes, and weights w_k.

Interpolatory Quadrature:

Previous subject area had been about polynomial interpolation. The integral of the polynomial intepolant p_n(x) of the f(x) maybe be used to derive the quadrature formula of F(x). Lets this polynomial integral be denoted as \hat{F(x)} from this point on.

\hat{F(x)}=\hat{F(x)}=\int_a^b {p_n(x)}dx$$

Example using Lagrange form:

The polynomial interpolation in the Lagrange form is: p_n(x) = \sum_{k=0}^n {L_{n,k}(x) f(x_k)}. The integral of a given function f(x) which is interpolated in Lagrange form can be approximated by:

Interchanging the order of summation and integration, we obtain:

The quadrature formula is then:

the weight is now w_k = \int_a^b {L_{n,k}(x) dx}. The most commonly used and the simplest quadrature formula are those with equally spaced nodes and its knows as the (n+1) Newton-Cotes formulas. Here is an example of 2-point closed Newton-Cotes formula using the Lagrange form above:

set x₀=a , x₁=b and \Delta{x}=h=b-a. The interpolatory quadrature is equivalent to \hat{F_1 (x)} = w_0 f(x₀) + w_1 f(x₁), where

From my previous on "Lagrange Interpolation",

w_0 = \int_{x₀}^{x₁} {\frac{x-x₁}{x₀-x₁}dx} and w_1 = \int_{x₀}^{x₁} {\frac{x-x₀}{x₁-x₀}dx}

Here is the Interpolatory quadrature:

function [sums]=sidrewquad(B,lend,rend)
  syms x;
    for i=1:length(B)
    c(i)=B(i)*x^(i-1);
    ss(i)=int(c(i),x,lend,rend);
  end

  sums=sum(ss);
  sums=double(sums);
end

This function file take the Interpolation polynomial and integrate it. Using this code to integrate f(x)= sin(x) on [- \pi, \pi], the numerical result of the code is:

sidrewquad(B,lend,rend)
ans =
    4.118443451920250e-11

sidrewquad(B,lend,rend)
ans =

This form can be extended for n+1 point of any Polynomial base interpolation:

The above equation can be used for any of the Rienmann Sum rules.

The above function file will work for any of the polynomial interpolation code, but not cubic spline. So we wrote one to Integrate cubic spline:

function [sums]=sidrewquadspline(a,b,c,d,lend,rend)
  syms x;
  for i = 1 : length(b)
    C(i)=(a(i)*x^(i-1));
   ss(i)=int(C(i),x,lend,rend);
  end

  sums=sum(ss);
  sums=double(sums);

Problem with this code is that we forgot cubic spline has a number of cubic polynomials, not just one. So this code is currently not working.

Gauss-Chebyshev Quadrature

The Newton-Cotes formulas are based on interpolating polynomials. However, interpolation at equidistant points experiences the Runge Phenomenon even for well-behaved function. For explicitly given f(x), other choices of quadrature should be considered, such as the Gaussian quadrature.

where

Here is a Matlab code implementing Gauss-Chebyshev Quadrature:

function int = chebquad(f, n, a, b) ;
  x = cos (((2 * (1 : n) – 1) / (2 * n)) * pi) ;
  w = pi / n ;
  fx = ((b - a) / 2) * f(((b - a) / 2) * x + (b + a) / 2) .* sqrt (1 - x .^ 2) ;
  int = sum (w .* fx) ;
end