⛏️

Numerical analysis

TagsComputer scienceMath
Created
Updated

Prerequisites

Resources

Resources

NameTags
Burden, Richard L., and J. Douglas Faires. Numerical Analysis. 7th ed. Belmont, CA: Brooks Cole, 2000. ISBN: 0534382169.numerical analysis
http://www.sam.math.ethz.ch/~hiptmair/tmp/NumCSE/NumCSE15.pdf
https://ocw.mit.edu/courses/18-335j-introduction-to-numerical-methods-spring-2019/

Key questions

Numerical analysis, numerical methods, or symbolic calculation

Differences

NameDifferencesSimilarities
Numerical analysisProofRigorousnumerical methodsAlgorithms
Numerical methodsAlgorithms
Symbolic calculationDeductive

Octave, Matlab, Mathematica, Maxima, Python, or Julia

Story

William Kahan

[1] https://archive.siam.org/meetings/la09/talks/benzi.pdf

https://mathshistory.st-andrews.ac.uk/Biographies/Seidel/

Octave

Glosario

Modelo matemático. Descripción de un sistema (fenómeno, procedimiento, procedimiento …) usando la teoría (gramatica, principios, sus reglas) de alguna de las áreas de las matemáticas (análisis numérico, ecuaciones diferenciales, teoría de automatas, probabilidades, topología, …) .

Exactitud. Cualidad de conocimiento sobre la magnitud real.

Precision.  Dispersión del conjunto de valores obtenidos o medidos de un instrumento.

Sesgo o inexacitud. Diferencia entre la magnitud real y los valores obtenidos.

Incertidumbre o imprecision. Dispersión tendiendo a cero de los valores obtenidos de un instrumento.

Accuracy and precision

Bias.

ISO 5725-1 ISO/DIS 5725-1

Mathematical preliminaries and Error Analysis

1.1 Review of Calculus

Worked examples

  1. Show that the following equations have at least one solution in the given intervals.
    • xcos(x)−2x2+3x−1=0,[0.2,0.3] and [1.2,1.3]xcos(x)-2x^2+3x-1=0,[0.2,0.3]\text{ and }[1.2,1.3]

      As f(x)=xcos(x)−2x2+3x−1f(x)=xcos(x)-2x^2+3x-1 is continuous and

      K=0K=0,

      f(0.2)=−0.28≤0≤f(0.3)=0.006009f(0.2) =-0.28\le 0 \le f(0.3)=0.006009 \\

      The Intermediate Value Theorem implies that a number x exists, with 0.2<c1<0.30.2<c_1<0.3.

      If f∈C[1.2,1.3]f\in C[1.2,1.3] and 0 is a number between f(1.2)=0.1548f(1.2)=0.1548 and f(1.3)=−0.1322f(1.3)=-0.1322, then there exists a number cc in (1.2,1.3)(1.2,1.3) for which f(c)=0f(c)=0.

https://gist.github.com/sanchezcarlosjr/b7d9adee444ea4bbd1864fa8c0386dd1
  1. Find intervals containing solutions to the following equations.
    • x−3−x=0x-3^{-x}=0

1.2 Round-off Errors and Computer arithmetic

Def. Floor function

⌊x⌋=floor(x)=min{n∈Z∣n≤x}\lfloor x \rfloor=floor(x)=min\{n\in Z|n\le x\}

Def. Chopping or truncating is a function such that

truncate(x,n)=⌊10nx⌋10ntruncate(x,n)=\dfrac{\left \lfloor {10^nx}\right \rfloor}{10^n}

Def. Suppose that p∗p^* is an approximation to pp.

Absolute error:∣p−p∗∣\text{Absolute error}: |p-p^*|
Relative error:∣p−p∗∣∣p∣,p≠0\text{Relative error}: \dfrac{|p-p^*|}{|p|},p\neq0
∣p−p∗∣∣p∣=∣1−p∗p∣\dfrac{|p-p^*|}{|p|}=|1-\dfrac{p^*}{p}|

Ceil and floor

https://stackoverflow.com/questions/6208488/implementation-of-ceil-and-floor

Cast

https://stackoverflow.com/questions/22330575/how-are-typecasts-in-c-implemented

https://en.wikipedia.org/wiki/Loop_unrolling

How Many Decimals Do We Really Need?

Floating-Point Computation by Pat Sterbenz

How calculating more precise than double or long double? How to increase precision by hundreds of decimals?

Use C++, and then use https://gmplib.org/#DOWNLOAD, other languages https://en.wikipedia.org/wiki/GNU_Multiple_Precision_Arithmetic_Library

How can minimize the error upon arithmetic operations?

Worked examples

  1. Determine the six-digit a) chopping and b) rounding values of the ϕ,e\phi,e.
https://gist.github.com/sanchezcarlosjr/de302811a791735664a65b150534571a

Chapter notes

https://nptel.ac.in/content/storage2/courses/122104019/numerical-analysis/Rathish-kumar/numerr/new4.htm

https://www.youtube.com/watch?v=JTinepDn1dI&ab_channel=OscarVeliz

Solutions of Equations in One Variable

Bisection method

Fixed-point iteration

https://gist.github.com/sanchezcarlosjr/c814b73dfcf4e415315fc170578d4807

Newton-Rapshon method

https://www.geogebra.org/classic/dpbxakmc

Summary

https://gist.github.com/sanchezcarlosjr/5ef537a436941173a5bc0cfd4c28cff1

Worked examples

  1. Da la definicion de numero de maquina en la forma de punto flotante decimal normalizado.

    Sea el punto flotante decimal normalizado analogo al punto flotante binario normalizado para representar las operaciones de los numeros de maquina segun el IEEE754 en base 10 y siendo inspirado en la normalizacion de la notacion cientifica (aqui tomamos d0=0d_0=0), a saber, el punto flotante decimal normalizado tiene la forma:

    ±0.d1d2d2d3...di...dk×10n\pm0.d_1d_2d_2d_3...d_i...d_k \times 10^n

    donde d1d_1 es distinto de 0 (d1∈N tal que 1≤d1≤9d_1\in N \text{ tal que }1\le d_1\le9) y los siguientes terminos del significando o mantisa did_i son cualquier numero menor a la base 10 (di∈N tal que 0≤di≤9∨2≤i≤kd_i \in N \text{ tal que } 0\le d_i \le 9 \lor2\le i\le k), 10 es nuestra base bb, nn el expontente y kk digitos decimales del numero de maquina tal que

    ±(d0+d110−1+...+dk−110−k+1)10n=y,y∈R\pm(d_0+d_110^{-1}+...+d_{k-1}10^{-k+1})10^{n}=y,y\in R

    Ejemplo:

    y=10.15y=10.15

    normalize(y)=0.1015×102normalize(y)=0.1015\times10^2

    porque

    (0+1×10−1+0×100+1×10−2+5×10−3)×102=(0+1\times10^{-1}+0\times{10^0}+1\times10^{-2}+5\times10^{-3})\times10^{2}=

    10+0+0.1+0.5=10.1510+0+0.1+0.5=10.15

    2.a

    Si limx→0xcosx−sinx=limx→0x−sinx=0lim_{x\to0}xcosx-sinx=lim_{x\to0}x-sinx=0,

    entonces por regla de L’Hopital (*)

    limx→0xcosx−sinxx−sinx=limx→0ddx(xcosx−sinx)ddx(x−sinx)=limx→0cosx−xsinx−cosx1−cosx=limx→0−xsinx1−cosxlim_{x\to0}\dfrac{xcosx-sinx}{x-sinx}=\\ lim_{x\to0}\dfrac{\dfrac{d}{dx}(xcosx-sinx)}{\dfrac{d}{dx}(x-sinx)}=\\ lim_{x\to0}\dfrac{cosx-xsinx-cosx}{1-cosx}=\\ lim_{x\to0}\dfrac{-xsinx}{1-cosx}

    Si limx→0−xsinx=limx→01−cosx=0lim_{x\to0}-xsinx=lim_{x\to0}1-cosx=0,

    entonces por *

    limx→0−xsinx1−cosx=limx→0ddx(−xsinx)ddx(1−cosx)=limx→0−sinx−xcosxsinx=limx→0−1−xcosxsinx=−1−limx→0cosxsinx=−1−limx→0cosxlimx→0sinx=−1−11=−2lim_{x\to0}\dfrac{-xsinx}{1-cosx}=lim_{x\to0}\dfrac{\dfrac{d}{dx}(-xsinx)}{\dfrac{d}{dx}(1-cosx)}\\ =lim_{x\to0}\dfrac{-sinx-xcosx}{sinx}\\ =lim_{x\to0}-1-\dfrac{xcosx}{sinx}\\ =-1-lim_{x\to0}\dfrac{cosx}{\dfrac{sin}{x}}\\ =-1-\dfrac{lim_{x\to0}cosx}{lim_{x\to0}\dfrac{sin}{x}}\\= -1-\dfrac{1}{1}\\=-2

    2. b

Redondeo a 4 digitos.

f(0.1000×100)=(0.1000×100)cos(0.1000×100)−sin(0.1000×100)(0.1000×100)−sin(0.1000×100)=(0.1000×100)(0.9950×100)−(0.9983×10−1)(0.1000×100)−(0.9983×10−1)=(0.9950×10−1)−(0.9983×10−1)0.1669×10−3=−0.3300×10−30.1669×10−3=−0.1977×101f(0.1000\times 10^0)\\=\dfrac{(0.1000\times 10^0)cos(0.1000\times 10^0)-sin(0.1000\times 10^0)}{(0.1000\times 10^0)-sin(0.1000\times 10^0)}\\ =\dfrac{(0.1000\times 10^0)(0.9950\times10^0)-(0.9983\times10^{-1})}{(0.1000\times 10^0)-(0.9983\times10^{-1})}\\ =\dfrac{(0.9950\times 10^{-1})-(0.9983\times10^{-1})}{0.1669\times10^{-3}}\\ =\dfrac{-0.3300\times10^{-3}}{0.1669\times10^{-3}}\\ =-0.1977\times10^{1}

2.c

cos(0.1000×100)=(0.1000×100)−(0.1000×100)2(0.2000×101)+O(x4)=(0.1000×100)−(0.1000×10−1)(0.2000×101)+O(x4)==(0.1000×100)−(0.5000×10−2)+O(x4)=0.9950×100+O(x4)cos(0.1000\times 10^0)=(0.1000\times10^0)-\dfrac{(0.1000\times10^0)^2}{(0.2000\times10^1)}+O(x^4)=\\ (0.1000\times10^0)-\dfrac{(0.1000\times10^{-1})}{(0.2000\times10^1)}+O(x^4)=\\ =(0.1000\times10^0)-(0.5000\times10^{-2})+O(x^4)\\ =0.9950\times10^{0}+O(x^4)

sin(0.1000×100)=(0.1000×100)−(0.1000×100)30.6000×101+O(x5)=(0.1000×100)−(0.1000×10−2)0.6000×101+O(x5)=(0.1000×100)−(0.1667×10−3)+O(x5)=(0.9983×10−1)+O(x5)sin(0.1000\times10^0)=(0.1000\times10^0)-\dfrac{(0.1000\times10^0)^3}{0.6000\times 10^1}+O(x^5)\\ = (0.1000\times10^0)-\dfrac{(0.1000\times10^{-2})}{0.6000\times 10^1}+O(x^5)=\\ (0.1000\times10^0)-(0.1667\times10^{-3})+O(x^5)= (0.9983\times10^{-1})+O(x^5)

f(0.1000×100)=(0.1000×100)(0.9950×100)−(0.9983×10−1)(0.1000×100)−(0.9983×10−1)=(0.9950×10−1)−(0.9983×10−1)(0.1000×100)−(0.9983×10−1)=(−0.3300×10−3)(0.1700×10−3)=−0.1941×101f(0.1000\times 10^0)=\dfrac{(0.1000\times 10^0)(0.9950\times10^{0})-(0.9983\times10^{-1})}{(0.1000\times 10^0)-(0.9983\times10^{-1})} =\\ \dfrac{(0.9950\times10^{-1})-(0.9983\times10^{-1})}{(0.1000\times 10^0)-(0.9983\times10^{-1})}=\\ \dfrac{(-0.3300\times 10^{-3})}{(0.1700\times10^{-3})}=\\ -0.1941\times10^{1}

  1. d
ERf(0.1)=∣−0.199899998×101+0.1977×101∣∣−0.199899998×101∣=0.1113×10−1E^{f (0.1)} _R=\dfrac{|-0.199899998\times 10^1+0.1977\times10^{1}|}{|-0.199899998\times 10^1|}\\ =0.1113\times10^{-1}

2.e

ERf(0.1)=∣−0.199899998×101+0.1941×101∣∣−0.199899998×101∣=0.2988×10−1E^{f (0.1)} _R=\dfrac{|-0.199899998\times 10^1+0.1941\times10^{1}|}{|-0.199899998\times 10^1|}\\ =0.2988\times10^{-1}
  1. Usa el metodo de Steffensen para encontrar una aproximacion a la raiz de
x=2−ex+x23x=\dfrac{2-e^x+x^2}{3}

con p0=−1 p_0 = −1, TOL=1×10−5TOL = 1 \times 10^{−5} y usa al menos 10 digitos en los calculos.

Vas a poner las formulas o ecuaciones que usaste en cada uno de los calculos aritmeticos que realizaste.

  1. Find root of a number using Newton’s method

https://opensource.apple.com/source/Libm/Libm-47.1/ppc.subproj/sqrt.c.auto.html

https://math.mit.edu/~stevenj/18.335/newton-sqrt.pdf

https://math.stackexchange.com/questions/296102/fastest-square-root-algorithm

https://stackoverflow.com/questions/63450864/ieee-754-conformant-sqrt-implementation-for-double-type

Tabla comparativa

MetodoDiferenciasVentajasDesventajas
BiseccionBasado en el teorema del valor intermedio, determinando cualquier intervalo lo satisface hasta lograr la tolerancia deseada. Por tanto, necesita un intervalo [a,b]. a=1a=1El mas simple de implementar. No necesita calculo de derivada. Siempre converge a la raiz.El mas lento por su radio de convergencia. No encuentra raices complejas. Buenas aproximaciones pueden ser rechazadas.
Netwon-RapshonMetodo de punto fijo, que necesita del calculo de la derivada. Necesitando un punto p0p_0 como entrada. a=2a=2 siempre que M=El mas rapido, su radio de convergencia es a=2a=2Necesita calculo previo de derivada. No encuentra raices complejas. f’(x)=0f’(x)=0  el metodo diverge. Depende en gran medida de p0p0 para converger (no esta en el intervalo[p−δ,p+δ] [p-\delta, p+\delta])
SecanteMetodo de punto fijo, que NO necesita del calculo de la derivada, usando su aproximacion mediante la def. Necesitando dos puntos de la secante p0,p1p_0,p_1 como entrada. a=ϕa=\phiNo necesita calculo de derivada.No encuentra raices complejas. Puede diverger.
Regula falsiAdaptacion del metodo de la secante y biseccion que asegura que la raiz este en el intervalo [p0,p1][p_0,p_1]. a=1a=1No necesita calculo de derivada. No encuentra raices complejas. No es recomendado por (Burden, 73).
MüllerBasado en metodo de la secante, pero en vez de una linea usamos una parabola, por tanto, necesitamos 3 puntos. p0,p1,p2p_0,p_1,p_2. a=1.87a=1.87No necesita calculo de derivada. Calcula raices complejas.Puede diverger.
Halley
Householder-n

http://www.math.usm.edu/lambers/mat460/fall09/lecture7.pdf

https://math.stackexchange.com/questions/2473794/asymptotic-order-versus-order-of-convergence

https://math.stackexchange.com/questions/976942/whats-the-formal-difference-between-analytical-and-numerical?rq=1

Interpolation

https://github.com/sanchezcarlosjr/uabc/blob/main/essays/numerical-analysis/uabc_MA_Interpolaci_n_por_el_m_todo_de_diferencias_divididas_de_Newton.pdf

https://gist.github.com/sanchezcarlosjr/6cd5f027740e964c4a45cb0880dc18ae

https://gist.github.com/sanchezcarlosjr/74ab156989087338a6d27227817a7a22

Numerical linear algebra

Solving Linear Systems

https://gist.github.com/sanchezcarlosjr/9edad3c66853ca24de3e97104cf68261

Gauss-Jordan

LU factorization

Approximate solutions by iterations

Exact solutions tend to be inefficiency O(n3)O(n^3). Our computers save us by their power processing in complex calculations. Thus the methods were born.

We say that a method approximate to the solution xx from an arbitrary initial vector x[0]x^{[0]} where each k iteration tend to x[0],x[1],x[2],...,x[k],...→xx^{[0]},x^{[1]},x^{[2]},...,x^{[k]},...\rarr x such that Ax=bAx=b.

The method converges when limk→∞x[k]=xlim_{k\rarr \infty} x^{[k]}=x, otherwise it diverges.

Jacobi method

This method was discovered by Carl Gustav Jacob Jacobi, which use triangular matrices.

https://sci-hub.do/https://doi.org/10.1080/0020739980290204

https://arxiv.org/pdf/1304.2097.pdf

https://www.coe.pku.edu.cn/attachments/dc028aad92c84efdafe6b9f65f1f0d48.pdf

https://core.ac.uk/download/pdf/82202439.pdf

Seidel method

This method was discovered by Philipp Ludwig von Seidel. We now (inappropriately) call the Gauss-Seidel method [1], which improves the Jacobi method. Honor to whom honor is due.

Methods of Conjugate Gradients

SVD

https://cims.nyu.edu/~donev/Teaching/NMI-Fall2014/Lecture5.handout.pdf

Elementary function computations

For computations with very high precision, the algorithms based on the elliptic integrals are among the fastest known today for the logarithm, the number π, and the exponential function.

Logarithms, pow, gamma,

Notes

CORDIC

https://www.jjj.de/fxt/fxtbook.pdf

Numerical differential equations

https://gist.github.com/sanchezcarlosjr/80aa8c039d71bea3decc86fb744f719f

Numerical analysis, numerical methods, or symbolic calculation