⛏️

Scientific Computing

Prerequisites

Introduction

Approximation of Real number to some procession.

Pay for precision with time.

Resources

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+3x1=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+3x1f(x)=xcos(x)-2x^2+3x-1 is continuous and

      K=0K=0,

      f(0.2)=0.280f(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 fC[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.

  1. Find intervals containing solutions to the following equations.
    • x3x=0x-3^{-x}=0

1.2 Round-off Errors and Computer arithmetic

Def. Floor function

x=floor(x)=min{nZnx}\lfloor x \rfloor=floor(x)=min\{n\in Z|n\le x\}

Def. Chopping or truncating is a function such that

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

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

Absolute error:pp\text{Absolute error}: |p-p^*|
Relative error:ppp,p0\text{Relative error}: \dfrac{|p-p^*|}{|p|},p\neq0
ppp=1pp\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

Leijen, Daan (December 3, 2001). "Division and Modulus for Computer Scientists" (PDF). Retrieved 2014-12-25.

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.

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

Is there any way to find the number of real roots of a polynomial computationally?

Bisection method

Fixed-point iteration

Newton-Rapshon method

Summary

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 (d1N tal que 1d19d_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 (diN tal que 0di92ikd_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+d1101+...+dk110k+1)10n=y,yR\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×101+0×100+1×102+5×103)×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 limx0xcosxsinx=limx0xsinx=0lim_{x\to0}xcosx-sinx=lim_{x\to0}x-sinx=0,

    entonces por regla de L’Hopital (*)

    limx0xcosxsinxxsinx=limx0ddx(xcosxsinx)ddx(xsinx)=limx0cosxxsinxcosx1cosx=limx0xsinx1cosxlim_{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 limx0xsinx=limx01cosx=0lim_{x\to0}-xsinx=lim_{x\to0}1-cosx=0,

    entonces por *

    limx0xsinx1cosx=limx0ddx(xsinx)ddx(1cosx)=limx0sinxxcosxsinx=limx01xcosxsinx=1limx0cosxsinx=1limx0cosxlimx0sinx=111=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×101)(0.1000×100)(0.9983×101)=(0.9950×101)(0.9983×101)0.1669×103=0.3300×1030.1669×103=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×101)(0.2000×101)+O(x4)==(0.1000×100)(0.5000×102)+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×102)0.6000×101+O(x5)=(0.1000×100)(0.1667×103)+O(x5)=(0.9983×101)+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×101)(0.1000×100)(0.9983×101)=(0.9950×101)(0.9983×101)(0.1000×100)(0.9983×101)=(0.3300×103)(0.1700×103)=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×1010.199899998×101=0.1113×101E^{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×1010.199899998×101=0.2988×101E^{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=2ex+x23x=\dfrac{2-e^x+x^2}{3}

con p0=1 p_0 = −1, TOL=1×105TOL = 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=2
Necesita 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=1
No 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.87
No 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

Kochanek–Bartels spline

Bezier curves

https://docs.scipy.org/doc/scipy/tutorial/interpolate.html

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 limkx[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

Los métodos más comunes para resolver este tipo de problemas incluyen el método de las diferencias finitas, el método de los elementos finitos, y el método de Runge-Kutta.

Numerical analysis, numerical methods, or symbolic calculation

Random number generation

Log

  1. Count the number of digits

By mod. Its time complexity is O(n)O(n)

def algorithmMod(n):
  count=0
  while(n>0):
    count=count+1
    n=n//10
  return count

By log. Its time complexity is O(1)O(1)

length(number)=log10(number)+1length(number)=\lfloor log_{10} (|number|)\rfloor+1
def algorithmLog(number):
  return math.floor(math.log10(math.abs(number)))+1

https://replit.com/join/zpwivezz-carlossanchez14

Modulo operation

Remainder after division.

a=nq+r,qZr=mod(dividend,divisor)=remainder=anana=nq+r,q\in\mathbb{Z}\\ r=mod(dividend,divisor)=remainder=a-n\lfloor \dfrac{a}{n} \rfloor

Worked examples

  1. Extracting individual digits.
    get(position,number)=trunc(number10position) mod 10get(position,number)=trunc(\dfrac{number}{10^{position}})\text{ }mod\text{ }10

Truncate

Floor and ceil

Rounding

Round half to even, a variant of the round-to-nearest method.

This method is called the round-to-even method. Other names include the round-half-to-even method, the round-ties-to-even method, and "bankers' rounding." This variant of the round-to-nearest method is also called convergent rounding, statistician's rounding, Dutch rounding, Gaussian rounding, odd–even rounding, or bankers' rounding.

Banker's rounding: the value is rounded to the nearest even number. Also known as "Gaussian rounding", and, in German, "mathematische Rundung".

Standard rounding: the value is rounded to the nearest number (be it odd or even). In German it is known as "kaufmännische Rundung".

Banker's rounding
Banker's rounding VBspeed © 2000-10, updated: 30-May-2002 Banker's rounding Banker's rounding: the value is rounded to the nearest even number. Also known as "Gaussian rounding", and, in German, "mathematische Rundung". Standard rounding: the value is rounded to the nearest number (be it odd or even). In German it is known as "kaufmnnische Rundung".
http://xbeat.net/vbspeed/i_BankersRounding.htm

754-2019 - IEEE Standard for Floating-Point Arithmetic since 1985.

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

If the fractional part of x is 0.5, then y is the even integer nearest to x.

⚠️
0 is even number.
round(1.7)=2,round(2.7)=3round(1.5)=2,round(2.5)=2round(1.3)=1,round(2.3)=2round(1.7)=2,round(2.7)=3round(1.5)=2,round(2.5)=2round(1.3)=1,round(2.3)=2round( 1.7 )=2,round( 2.7 )=3\\ round( 1.5 )=2,round( 2.5 )=2\\ round( 1.3 )=1, round( 2.3 )= 2\\ round( -1.7 )=-2,round( -2.7 )=-3\\ round( -1.5 )=-2,round( -2.5 )=-2\\ round( -1.3 )=-1, round( -2.3 )=-2
round(x)={xif xx<0.5 or (isEven(x) and xx=0.5)xif xx>0.5 or (isOdd(x) and xx=0.5)round(x) = \begin{cases} \lfloor x \rfloor &\text{if } |\lfloor x \rfloor -x|<0.5 \text{ or (}isEven(x)\text{ and }|\lfloor x \rfloor -x|=0.5 )\\ \lceil x \rceil &\text{if } |\lceil x \rceil -x|>0.5 \text{ or (}isOdd(x)\text{ and }|\lceil x \rceil -x|=0.5 ) \end{cases}
JSFiddle
Test your JavaScript, CSS, HTML or CoffeeScript online with JSFiddle code editor.
https://jsfiddle.net/2hs9nqbk/2/
function roundIt(n, d = 0) {
            var m = Math.pow(10, d); 
            var n = +(d ? n * m : n).toFixed(8);
            var i = Math.floor(n), 
                diff = n - i; // getting the difference 
            var e = 1e-8; // Rounding errors in var(diff) 
            // Checking if the difference is less than or 
            // greater than, based on that adding the 1 to it. 
            var r = (diff > 0.5 - e && diff < 0.5 + e) ? 
                ((i % 2 == 0) ? i : i + 1) : Math.round(n);
            return d ? r / m : r; // if d != 0 then returning r/m else r 
}

https://www.geeksforgeeks.org/gaussian-bankers-rounding-in-javascript/

Experimental

import time					
import numpy as np

def timer(f):
  x=np.random.rand(1,100000)[0]
  times = []
  for i in range(10):
    tic = time.perf_counter()
    f(x)
    toc = time.perf_counter()
    times.append(toc - tic)
  print(f"Build finished in {np.mean(times):0.4f} +- {np.std(times):0.4f} seconds")

Seminumerical Algorithms. Arithmetic.

Numerical analysis.

Matrix operations.

The Fast Fourier Transform

Integer and Polynomial Arithmetic

Number-Theoretic Algorithms

References

Burden, R. L. y Faires, J. D. (2015). Análisis numérico. (9a ed.).
Cengage Learning.


Chapra, S. C. (2017). Applied numerical methods with MATLAB for engineers and scientists. (4a ed.). Mcgraw-Hill Education.


Demanet, L., (2012). Introduction to numerical analysis. MIT
opencourseware.
https://ocw.mit.edu/courses/mathematics/18-330introduction-to-numerical-analysis-spring-2012/


Epperson, J. F., (2021). An introduction to numerical methods and analysis. (3a ed). Wiley.


Gezerlis, A., (2020). Numerical methods in physics with python. Cambridge University Press.


Howard II, J. P., (2020). Computational methods for numerical analysis with R. CRC Press

Fausett, L.V. (2002). Numerical methods: algorithms and applications. Pearson.


Gerald, C.F. y Wheatley, P.O. (2001). Análisis numérico con aplicaciones (6a ed.). Prentice Hall.


Stoer, J. & Bulirsch, R. (1993). Introduction to numerical analysis.Springer-Verlag.