«

Apunte: Elementos de cálculo numérico

$\def\lim{\mathop{\text{lím}}\limits} \def\sup{\mathop{\mathrm{sup}}\limits} \def\inf{\mathop{\text{ínf}}\limits} \def\max{\mathop{\text{máx}}\limits} \def\min{\mathop{\text{mín}}\limits} \def\limsup{\mathop{\text{lím}\,\text{sup}}\limits} \def\liminf{\mathop{\text{lím}\,\text{ínf}}\limits} \def\varlimsup{\mathop{\overline{\text{lím}}}\limits} \def\varliminf{\mathop{\underline{\text{lím}}}\limits} \def\sen{\mathop{\mathrm{sen}}} \def\ord{\mathop{\mbox{\rm ord}}} \def\sig{\mathop{\text{sig}}} \def\adj{\mathop{\mbox{\rm adj}}} \def\Im{\mathop{\mathrm{Im}}} \def\im{\mathop{\mathrm{im}}} \def\coker{\mathop{\mathrm{coker}}} \def\codim{\mathop{\mathrm{codim}}} \def\Hom{\mathop{\mathrm{Hom}}} \def\Nu{\mathop{\mathrm{Nu}}} \def\Dom{\mathop{\mathrm{Dom}}} \def\mcm{\mathop{\mathrm{mcm}}} \def\mcd{\mathop{\mathrm{mcd}}} \def\Obj{\mathop{\mathrm{Obj}}} \def\End{\mathop{\mathrm{End}}} \def\Aut{\mathop{\mathrm{Aut}}} \def\Gal{\mathop{\mathrm{Gal}}} \def\car{\mathop{\mathrm{car}}} \def\Spec{\mathop{\mathrm{Spec}}} \def\nil{\mathop{\mathrm{nil}}} \def\ann{\mathop{\mathrm{ann}}} \def\Tr{\mathop{\mathrm{Tr}}} \def\Res{\mathop{\mathrm{Res}}} \def\dive{\mathrel{||}} \def\ldot{\mathrel{.}} \def\symdiff{\mathrel{\triangle}} \def\vacio{\varnothing} \def\blacksquare{\hbox{\vrule width 5pt height 10pt depth 0pt}} \def\qed{\ \ \ \hbox{}\nolinebreak $\blacksquare$ \par{}\medskip} \def\multint{\int\!\!\cdot\!\cdot\!\cdot\!\!\int} \newcommand{\pila}[2]{\genfrac{}{}{0pt}{}{#1}{#2}} \newcommand{\comb}[2]{\binom{#1}{#2}} \newcommand{\parti}[1]{\frac{\partial}{\partial{#1}}} \newcommand{\partii}[2]{\frac{\partial{#1}}{\partial{#2}}} \def\syss{\Leftrightarrow} \def\then{\Rightarrow} \def\tri{\bigtriangleup} \def\medio{\textstyle{\frac12}} \def\ms{\medskip} \def\mm{\mathrm} \def\ni{\noindent} \def\la{\langle} \def\ra{\rangle} \def\flecha{\mathop{\to}\limits} \def\longflecha{\mathop{\longrightarrow}\limits} \def\longinvflecha{\mathop{\longleftarrow}\limits} \def\1{\mathbbm{1}} \def\Re{\mathop{\mm{Re}}} \def\ext{\mathop{\mm{ext}}} \def\co{\mathop{\mm{co}}} \def\pr{\mathop{\mm{pr}}} \def\id{\mm{id}} \def\ind{\mathop{\mm{ind}}} \def\rg{\mathop{\mm{rg}}} \def\GL{\mm{GL}} \def\cte{\mm{cte}} \def\ev{\mm{ev}} %Cosos que uso: % \smallsetminus es el \ de conjuntos % \leqq, \geqq, \lneqq, \gneqq, \subset, \supset, \subsetneqq, \supsetneqq % \overline, \hat, \underbrace, \overbrace % \rtimes, \ltimes producto semidirecto % \coprod \amalg % \hookrightarrow \twoheadrightarrow \dashrightarrow \leadsto \def\menos{\smallsetminus} \def\cdtos{\cdots} \def\bl{\as\begin{leftbar}} \def\el{\end{leftbar}} \def\N{\mathbb{N}} \def\R{\mathbb{R}} \def\C{\mathbb{C}} \def\Z{\mathbb{Z}} \def\HH{\mathcal{H}} \def\FF{\mathcal{F}} \def\PP{\mathcal{P}} \def\E{\mm{E}} \def\Var{\mm{Var}} \def\Cond{\mm{Cond}}$

Punto flotante

Los números de máquina son de la forma $0,a_1a_2\ldots a_m\cdot B^l=q\cdot B^l$, donde $a_1\neq0$, $a_i\in\{0,\ldots,B-1\}$, $B$ es la base y $-M_1\leqq l\leqq M_2$. Dado $x=0,a_1a_2\ldots\cdot B^l\in\R$ definimos el truncado $x^*=0,a_1a_2\ldots a_m \cdot B^l$; vale $\frac{|x-x^*|}{|x|}\leqq B^{-m+1}$. El redondeo cumple $\frac{|x-x^*|}{|x|}\leqq \frac12B^{-m+1}$. Tenemos pues $x^*=(1+\delta)x$, con $|\delta|\leqq\epsilon=\frac12B^{-m+1}$, donde $\epsilon$ se llama $\epsilon$ de máquina; es el menor número tal que $(1+\epsilon)^*\neq 1$. Sumar números muy distintos en magnitud es malo; restar números muy parecidos también; la mejor manera de sumar muchos números es de menor a mayor.

Normas matriciales

Equivalencia de normas. Las normas en $\R^n$ son equivalentes, es decir, dadas $\|{\cdot}\|$, $\|{\cdot}\|'$ hay $K_1,K_2>0$ con $K_1\|{\cdot}\|'\leqq\|{\cdot}\|\leqq K_2\|{\cdot}\|'$.

Prueba. Basta probarlo para $\|x\|'=\|x\|_2=\sqrt{\sum_{i=1}^nx_i^2}$; tenemos $\|\sum_{i=1}^nx_ie_i\|\leqq\sum_{i=1}^n|x_i|\|e_i\|\leqq \sqrt{\sum_{i=1}^n\|e_i\|}\|x\|_2$ por Cauchy-Schwarz; la otra por el absurdo: si no, para todo $n$ hay $x_n$ con $\|x_n\|<\frac{1}{n}\|x_n\|_2$; podemos tomar $\|x_n\|_2=1$ y $\|x_n\|<\frac{1}{n}$; luego hay una subsecuencia con $\|x_{n_k}-x\|_2\to0$, luego $\|x_{n_k}-x\|\to0$, $\|x\|=0$, $x=0$, absurdo.∎

Normas matriciales. Dada una norma $\|{\cdot}\|$ de $\R^n$ obtenemos una norma para $\R^{n\times n}$ dada por $\|A\|=\max_{\|x\|=1}\|Ax\|$, y tenemos $\|Ax\|\leqq \|A\|\|x\|$ y $\|AB\|\leqq\|A\|\|B\|$. Vale $\|A\|_\infty=\max_i\left\{\sum_{j=1}^n|A_{ij}|\right\}$, $\|A\|_1=\max_j\left\{\sum_{i=1\vphantom{j}}^n|A_{ij}|\right\}$, $\|A\|_2=\sqrt{\rho(A^tA)}$, donde $\rho(A)=|\lambda_{\max}|$ es el radio espectral, con $\lambda_{\max}\in\C$ el autovalor de norma máxima.

Prueba. Sea $v_1,\ldots,v_n$ una base ortonormal de autovectores de $A^tA$ con autovalores $\lambda_1,\ldots,\lambda_n$. Dada $x=\sum_{i=1}^na_iv_i$ tenemos $\|Ax\|_2^2=\la Ax,Ax\ra=\la A^tAx,x\ra=$ $\bigl\langle\sum_{i=1}^n\lambda_ia_iv_i,\sum_{i=1}^na_iv_i\bigl\rangle=\sum_{i=1}^n\lambda_ia_i^2\leqq |\lambda_{\max}|\|x\|_2^2$, con igualdad si $x=v_{\max}$, luego $\|A\|_2=\sqrt{\rho(A^tA)}$.∎

Número de condición. Si $A\in\R^{n\times n}$ es inversible definimos $\Cond(A)=\|A\|\|A^{-1}\|$. Vale $\Cond(A)\geqq 1$ y $\frac1{\Cond(A)}=\inf_{\text{$B$ singular}}\frac{\|A-B\|}{\|A\|}$. Aplicación: si $Ax=b$ y $A\tilde x=\tilde b$, $\frac{\|\tilde x-x\|}{\|x\|}\leqq\Cond(A)\frac{\|\tilde b-b\|}{\|b\|}$.

Prueba. Si $B$ es singular hay $x\neq0$ con $Bx=0$; vale $\|x\|=\|A^{-1}(A-B)x\|$ luego $1\leqq\|A^{-1}\|\|A-B\|$ y $\frac1{\Cond(A)}\leqq \frac{\|A-B\|}{\|A\|}$. Ahora tomo $y$ con $\|A^{-1}y\|=\|A^{-1}\|\|y\|$ y $x$ con $Ax=y$, $\|y\|=\frac{1}{\|A^{-1}\|}$ y $\|x\|=1$; tomo $z$ tal que $z^tu\leqq 1$ si $\|u\|=1$ pero $z^tx=1$ (hiperplano que pasa por $x$ y deja a la bola unitaria de un lado); $B=A-yz^t$ es singular porque $Bx=0$; $\|yz^tu\|\leqq\|y\|=\frac{1}{\|A^{-1}\|}$ si $\|u\|=1$ con igualdad si $u=x$, luego $\|A-B\|=\|yz^t\|=\frac1{\|A^{-1}\|}$.∎

Radio espectral. Dado $\epsilon>0$ hay una norma matricial $\|{\cdot}\|$ tal que $\rho(A)\leqq \|A\|\leqq \rho(A)+\epsilon$. Equivalen $A^n\to 0$ y $\rho(A)<1$. Además, $\|A^n\|^{\frac1n}\to\rho(A)$.

Prueba. Primero, sea $v_1,\ldots,v_n\in\C^n$ una base que pone a $A$ en su forma de Jordan; multiplicando los $v_i$ por potencias de $\epsilon$ logramos que los números fuera de la diagonal sean $\epsilon$. Defino la norma $\|\sum_{i=1}^n a_iv_i\|=\max_i|a_i|$; vale $\|A\|\leqq \rho(A)+\epsilon$ en esa norma. Segundo, si $\rho(A)<1$ con $\epsilon$ chico obtengo $\|A\|<1$, luego $\|A^n\|\leqq\|A\|^n\to0$; si $A^n\to 0$, sea $v_i$ con $Av_i=\lambda_{\max}v_i$; si $v_i=x+yi$, con $x,y\in\R^n$, $A^nv_i=A^nx+A^nyi\to0$, pero $\|A^nv_i\|=\rho(A)^n$, luego $\rho(A)^n\to0$ y $\rho(A)<1$. Tercero, de $\|A\|\leqq \rho(A)+\epsilon$ obtengo $\|A^n\|^{\frac1n}\leqq \rho(A)+\epsilon$; por otro lado $\rho(A)^n\leqq \rho(A^n)\leqq\|A^n\|$ luego $\rho(A)\leqq \|A^n\|^{\frac1n}\leqq \rho(A)+\epsilon$ si $n$ grande.∎

Matrices

Descomposición LU. Dada $A\in\R^{n\times n}$ hay una matriz de permutación $P$ con $PA=LU$, $L$ diagonal inferior y con 1 en la diagonal y $U$ diagonal superior.

Descomposición de Cholesky. Dada $A\in\R^{n\times n}$ simétrica y definida positiva (o sea $\la x,Ax\ra>0$ si $x\neq0$) existe una única $L$ diagonal inferior con elementos de la diagonal positivos tal que $A=LL^t$.

Prueba. Vale $A_{ij}=\sum_{k=1}^{\min\{i,j\}}L_{ik}L_{jk}$. Si tenemos los $L_{ij}$ con $j< k$, tenemos $A_{kk}=\sum_{j=1}^k L_{kj}^2$, luego $L_{kk}^2=A_{kk}-\sum_{i=1}^{k-1}L_{kj}^2$. Veamos que podemos tomar una raíz positiva: si no, tomamos raíz $L_{kk}\in\C$; tenemos que $A_k,L_k\in\R^{k\times k}$, con $A_{kij}=A_{ij}$ y $L_{kij}=L_{ij}$ cumplen $A_k=L_kL_k^t$, luego $0<|A_k|=|L_k|^2=L_{11}^2\ldots L_{kk}^2$ (por Sylvester), con lo que $L_{kk}^2>0$. Faltarían los $L_{ik}$ con $i>k$; ahora $A_{ik}=\sum_{j=1}^{k}L_{ij}L_{kj}$, luego $L_{ik}=\frac1{L_{kk}}(A_{ik}-\sum_{j=1}^{k-1}L_{ij}L_{kj})$.∎

Métodos de Jacobi y Gauss-Seidel. Dada $A\in\R^{n\times n}$ escribimos $A=L+D+U$, con $L_{ij}=0$ si $i\leqq j$, $D$ diagonal y $U_{ij}=0$ si $i\geqq j$. Queremos resolver $Ax=b$. Métodos iterativos: dado $x^0$ ponemos $x^{k+1}=Bx^k+c$, asumiendo $x=Bx+c$. Jacobi: $x^{k+1}=-D^{-1}(L+U)x^k+D^{-1}b$; Gauss-Seidel: $x^{k+1}=-(D+L)^{-1}Ux^k+(D+L)^{-1}b$. O sea Jacobi es $x^{k+1}_i=\frac1{a_{ii}}(b_i-\sum_{j\neq i}a_{ij}x^k_j)$ y Gauss-Seidel, $x^{k+1}_i=\frac1{a_{ii}}(b_i-\sum_{j<i}a_{ij}x^{k+1}_j-\sum_{j>i}a_{ij}x^k_j)$. El error $e_k=x^k-x$ cumple $e_{k+1}=Be_k=B^ke_1$, luego el método converge para todo punto inicial si y sólo si $B^k\to0$, es decir, sii $\rho(B)<1$. Si $A$ es estrictamente diagonal dominante, es decir, si $|a_{ii}|>\sum_{j\neq i}|a_{ij}|$ para $i=1,\ldots,n$, entonces tanto Jacobi como Gauss-Seidel convergen. Si $A$ es simétrica y definida positiva entonces Gauss-Seidel converge. Si $A$ es tridiagonal ($A_{ij}=0$ si $|i-j|>1$) entonces Jacobi converge si y sólo si Gauss-Seidel converge, y $\rho(B_{GS})=\rho(B_J)^2$.

Prueba. Primero, que GS converge si $A$ cumple $|a_{ii}|>\sum_{j\neq i}|a_{ij}|$: $B=-(D+L)^{-1}U$ cumple $B_{ij}=-\frac1{A_{ii}}(A_{ij}-\sum_{k=1}^{i-1}A_{ik}B_{kj})$ si $i< j$ y 0 si no; queremos ver que $\|B\|_\infty=\max_{i}\{\sum_{j>i}|B_{ij}|\}<1$; supongamos que sabemos que $\sum_{j>k}|B_{kj}|<1$ para $k< i$; $\sum_{j>i}|B_{ki}|\leqq\frac1{|A_{ii}|}\sum_{j>i}(|A_{ij}|+\sum_{k< i}|A_{ik}||B_{kj}|)=\frac1{|A_{ii}|}(\sum_{j>i}|A_{ij}|+\sum_{k<i}|A_{ik}|\sum_{j>i}|B_{kj}|)\leqq$ $\frac1{|A_{ii}|}(\sum_{j>i}|A_{ij}|+\sum_{k<i}|A_{ik}|)<1$ y vale para $i$, luego $\|B\|_\infty<1$. Segundo, que GS converge si $A$ es simétrica y definida positiva: sea $\lambda\in\C$, $x\in\C^n$, $x\neq0$, con $Bx=\lambda x$; es $-L^tx=\lambda(D+L)x$, $Ax=(1-\lambda)(D+L)x$; tenemos $\frac{1}{1-\lambda}=\frac{x^*(D+L)x}{x^*Ax}$ y conjugando $\frac{1}{1-\overline\lambda}=\frac{x^*(D+L^*)x}{x^*Ax}$; sumo: $2\Re(\frac{1}{1-\lambda})=1+\frac{x^*Dx}{x^*Ax}>1$, lo que da $|\lambda|<1$. Tercero, si $A$ es tridiagonal vale $|A|=|D+\alpha^{-1}L+\alpha U|$ si $\alpha\neq0$: $A+\alpha^{-1} L+\alpha U$ es $A$ en la base $v_1,\alpha v_2,\ldots,\alpha^{n-1}v_n$; con eso se ve que si $\lambda\neq0$, $\lambda^2$ es autovalor de $B_{GS}$ sii $\lambda$ es autovalor de $B_J$.∎

Raíces

Punto fijo. Dada $f:I\to I$, con $I=[a,b]$ intervalo, se busca $r$ con $f(r)=r$. El algoritmo: tomar $x_0$ y luego $x_n=f(x_{n-1})$. Si $|f'|\leqq\lambda<1$ converge y se tiene $|x_n-r|\leqq \frac{\lambda^n}{1-\lambda}|x_1-x_0|$.

Prueba. $|x_n-r|=|f(x_{n-1})-f(r)|\leqq \lambda|x_{n-1}-r|\leqq \cdtos\leqq \lambda^n|x_0-r|$. Por otro lado $|x_1-x_0|=|f(x_0)-x_0|\geqq |x_0-f(r)|-|f(r)-f(x_0)|\geqq |x_0-r|-\lambda|x_0-r|=(1-\lambda)|x_0-r|$.∎

Newton-Raphson. Dada $f:I\to\R$ consiste en buscar una raíz $r$ como punto fijo de $g(x)=x-\frac{f(x)}{f’(x)}$. Si $f''$ está acotada y $f'(r)\neq0$, hay un entorno de $r$ en el que el algoritmo converge con orden 2, es decir $\frac{|e_{n+1}|}{|e_n|^2}$ está acotado, con $e_n=x_n-r$. Si $f''>0$ converge partiendo de cualquier punto.

Prueba. Por Taylor $\left|\frac{e_{n+1}}{e_n^2}\right|=\frac12\frac{|f''(\xi)|}{|f'(x_n)|}$, que está acotado, lo que prueba lo primero. La convexidad $f''>0$ da por Taylor que si $f$ decrece en $I_1=(-\infty,v)$ y crece en $I_2=(v,\infty)$, si $x_0\in I_i$, $x_n\in I_i$ para todo $n$, y $x_n$ es creciente en $I_1$ y decreciente en $I_2$, por lo que converge en ambos casos.∎

Secante. Algoritmo: $x_{n+1}=x_n-f(x_n)\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}$. Si $r$ es raíz, $f''$ acotado y $f'(r)\neq0$, hay un entorno de $r$ en el que converge con orden $p>0$, $p^2=p+1$, es decir, $\frac{|e_{n+1}|}{|e_n|^p}$ está acotado.

Prueba. Se ve que $e_{n+1}=e_ne_{n-1}\frac{f[x_n,r,x_{n-1}]}{f[x_n,x_{n-1}]}=\frac12\frac{f''(\eta_n)}{f'(\xi_n)}e_ne_{n-1}$. Sea $y_n=\frac{|e_{n+1}|}{|e_n|^p}$, $c_n=\frac12\frac{f''(\eta_n)}{f'(\xi_n)}$ y $c=\frac12\frac{f''(r)}{f'(r)}$. Tenemos $y_n=|c_n|y_{n-1}^{1-p}$; si $c_n\to c\neq0$ se ve que $|y_n|\to |c|^{1/p}$.∎

Interpolación polinómica

Forma de Newton. Defino $f[x_0]=f(x_0)$ y $f[x_0,\ldots,x_n]=\frac{f[x_0,\ldots,x_{n-1}]-f[x_1,\ldots,x_n]}{x_0-x_n}$. Defino $\PP_n$ como el conjunto de polinomios de grado a lo sumo $n$. El $p_n\in\PP_n$ que interpola a $f$ en $x_0,\ldots,x_n$ es

$$p_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+\cdots+f[x_0,\ldots,x_n](x-x_0)\ldots(x-x_{n-1})$$

y vale $f(x)-p_n(x)=f[x_0,\ldots,x_n,x](x-x_0)\ldots(x-x_n)$ para todo $x$. Si $f\in C^n[a,b]$ hay $\xi$ con $f[x_0,\ldots,x_n]=\frac{f^{(n)}(\xi)}{n!}$.

Prueba. Lo primero es obvio por inducción, probando además que $f[x_0,\ldots,x_n]$ es simétrica. Para lo segundo $F(x)=f(x)-p_{n-1}(x)-f[x_0,\ldots,x_n](x-x_0)\ldots(x-x_{n-1})$ tiene $n+1$ raíces $x_0,\ldots,x_n$, luego hay $\xi$ con $F^{(n)}(\xi)=0$, es decir, $f^{(n)}(\xi)-f[x_0,\ldots,x_n]n!=0$.∎

Interpolación de Hermite. Sea $f\in C^n(I)$. La función $f[x_0,\ldots,x_n]:I^n\menos\bigcup_{i\neq j}\{x_i=x_j\}\to\R$ se puede extender a una función continua sobre $I^n$ calculada ordenando $(x_0,\ldots,x_n)$ y usando la fórmula recursiva con base $f[\underbrace{x,\ldots,x}_{\text{$k+1$ veces}}]=\frac{f^{(k)}(x)}{k!}$. Con eso dados puntos $x_0,\ldots,x_n$ obtenemos $p_n\in\PP_n$ por la fórmula de Newton y es el único que satisface $p_n^{(j)}(x)=f^{(j)}(x)$ si $x$ aparece $j+1$ veces. Sigue valiendo $f(x)-p_n(x)=f[x_0,\ldots,x_n,x](x-x_0)\ldots(x-x_n)$ y $f[x_0,\ldots,x_n]=\frac{f^{(n)}(\xi)}{n!}$.

Prueba. Que $p_n$ cumple la condición se ve reordenando $(x_i)$ para que los $x$ estén al principio (se ve que el $p_n$ obtenido es el mismo perturbando los $(x_i)$ para que sean distintos). Que es único: si $q_n$ cumple $q_n(x)-p_n(x)=q_n[x_0,\ldots,x_n,x](x-x_0)\ldots(x-x_n)$ pero $q_n[x_0,\ldots,x_n,x]=\frac{q_n^{(n+1)}(\xi)}{(n+1)!}=0$ (porque $q_n\in\PP_n$), entonces $q_n=p_n$.∎

Polinomios de Chebyshev. Sea $T_n(x)=\cos(n\cos^{-1}(x))$, donde $\cos^{-1}:[-1,1]\to[0,\pi]$. Restringido a $[-1,1]$ está en $\PP_n$, sigue la recursión $T_{n+1}=2xT_n-T_{n-1}$, $T_1=x$, $T_0=1$, tiene coeficiente principal $2^{n-1}$, tiene $n$ raíces distintas en $[-1,1]$, $|T_n|\leqq1$ en $[-1,1]$ y lo alcanza en $n+1$ puntos. Si $P$ es un polinomio mónico de grado $n$ entonces $\|P\|_{L^\infty[-1,1]}\geqq \frac1{2^{n-1}}$. Por lo tanto la mejor interpolación sobre $n$ puntos en $[-1,1]$ es sobre las raíces de $T_{n+1}$. Esto se lleva a $[a,b]$ tomando $\tilde T_n(x)=T_n(\frac{2x-a-b}{b-a})$. Se obtiene $\|f-p_n\|_{L^\infty[a,b]}\leqq \frac{\|f^{(n+1)}\|_{L^\infty[a,b]}}{2^{2n+1}(n+1)!}(b-a)^{n+1}$.

Prueba. Lo primero es claro. Supongamos que hay $P$ mónico de grado $n$ con $|P(x)|<\frac1{2^{n-1}}$ para todo $x\in[-1,1]$. Sea $W_n=\frac1{2^{n-1}}T_n$ y $-1=y_0<y_1<\cdots<y_n=1$ con $|W_n(y_i)|=\frac1{2^{n-1}}$, alternando de signo. Entonces $(P(y_i)-W_n(y_i))(P(y_{i+1})-W_n(y_{i+1}))<0$, luego $P-W_n$ tiene $n$ raíces, pero grado a lo sumo $n-1$; luego es nulo y $P=W_n$, absurdo.∎

Splines cúbicos. Dada $f\in C[a,b]$ y $a=x_0<x_1<\cdots<x_n=b$ hay una única $S\in C^2[a,b]$ tal que $S|_{[x_i,x_{i+1}]}$ es cúbica, $S(x_i)=f(x_i)$ y $S''(x_0)=S''(x_n)=0$.

Cuadrados mínimos

Productos internos. Sea $V$ un $\R$-espacio vectorial; un producto interno es una función bilineal $\langle{\cdot},{\cdot}\rangle:V^2\to\R$ simétrica tal que $\langle x,x\rangle\geqq0$, con igualdad sii $x=0$. Define una norma $\|x\|=\sqrt{\langle x,x\rangle}$. Vale Cauchy-Schwarz: $\langle x,y\rangle\leqq\|x\|\|y\|$, desigualdad triangular $\|x+y\|\leqq\|x\|+\|y\|$ y Pitágoras: $\|x+y\|^2=\|x\|^2+\|y\|^2$ si $\la x,y\ra=0$. Sea $S\subset V$ un subespacio, y $x\in V$; $y\in S$ minimiza $\|x-y\|$ sii para todo $s\in S$ vale $\langle x-y,s\rangle$; en ese caso $y$ es único. Si $S$ es finitamente generado se puede tomar una base ortonormal $v_1,\ldots,v_n$, con $\|v_i\|=1$ y $\la v_i,v_j\ra=0$ si $i\neq j$. En ese caso $y=\sum_{i=1}^n \la x,v_i\ra v_i$.

Prueba. Si $y\in S$ hace a $\|x-y\|$ mínimo, $\|x-y+ts\|\geqq \|x-y\|$ para $s\in S$, $t\in\R$, luego $t^2\langle s,s\rangle+2t\langle x-y,s\rangle\geqq0$, imposible si $\langle x-y,s\rangle\neq0$; además hay igualdad sii $s=0$, lo que dice que $y$ es único. Si $\langle x-y,s\rangle=0$ para todo $s\in S$, $\|x-y+s\|^2=\|x-y\|^2+\|s\|^2\geqq\|x-y\|^2$ y $\|x-y\|$ es mínimo. Dado $S$ con base $\{u_1,\ldots,u_n\}$ hacemos Gram-Schmidt: ponemos $v_k=u_k-\sum_{i=1}^{k-1}\frac{\la u_k,v_i\ra}{\la v_i,v_i\ra}v_i$ y queda una base ortonormal $e_n=\frac1{\|v_n\|}v_n$.∎

Regresión lineal. Dado $y\in\R^n$ y $x\in\R^{n\times k}$ queremos $\beta\in\R^k$ tal que $\sum_{i=1}^n \left(y_i-\sum_{j=1}^k\beta_jx_{ij}\right)^2$ sea mínimo. Se da sii $x^tx\beta=x^ty$. Si las columnas de $x$ son independientes la solución es única: $\beta=(x^tx)^{-1}x^ty$.

Lema. Si $A\in\R^{n\times k}, x\in\R^k, y\in\R^n$ entonces $\la Ax,y\ra=\la x,A^ty\ra$.

Prueba. $S=\{x\beta\mid \beta\in\R^k\}$ es un subespacio; queremos $\beta$ con $\|y-x\beta\|$ mínimo; luego $\la y-x\beta,xs\ra=0$ para todo $s\in\R^k$, pero por el lema es $\la x^t(y-x\beta),s\ra=0$, o sea $x^tx\beta=x^ty$. Si las columnas de $x$ son independientes, $\beta\mapsto x\beta$ es inyectiva, luego $\beta$ es única; además $\la x^tx\beta,\beta\ra=\la x\beta,x\beta\ra$, luego $x^tx\beta=0\then \beta=0$ y $x^tx$ es inversible.∎

Descomposición QR. Si $A\in\R^{n\times k}$, sean $a_1,\ldots,a_k$ sus columnas; tomo $u_i=a_i-\sum_{j=1}^{i-1}\frac{\la a_i,u_j\ra}{\la u_j,u_j\ra}u_j$ y $e_i=\frac1{\|u_i\|}u_i$; se ve $a_i=\sum_{j=1}^i \la a_i,e_j\ra e_j$; luego $A=QR$ con $Q_{ij}=e_{ji}$ y $R=\la a_i,e_j\ra$, $Q$ ortogonal (o sea $Q^tQ=1$) y $R$ triangular superior.

Aproximación por polinomios. Wierstrass: si $f\in C[0,1]$ hay una sucesión $p_n\in\PP_n$ que converge uniformemente: $p_n=\sum_{k=0}^n f\left(\frac kn\right)\comb{n}{k}x^k(1-x)^k$. Si $w\in C[0,1]$ es positiva, $C[0,1]$ tiene un producto interno $\la f,g\ra=\int_0^1 wfg$ y si $p_n^*$ es la proyección de $f$ en $\PP_n$ tenemos $\|f-p_n^*\|\to0$, y por lo tanto (Parseval) $\|f\|^2=\sum_{n=0}^\infty \la f, e_n\ra^2$ si $\{e_n\}$ es la base ortonormal de polinomios.

Prueba. Sea $f\in C[0,1]$. Sea $x\in[0,1]$; dado $n$ sean $X_n$ variables aleatorias iid con $P(X_n=1)=x$, $P(X_n=0)=1-x$, y sea $X=\frac1n\sum_{i=1}^nX_i$. Vale $\E(X)=x$, $\Var(X)=\frac1nx(1-x)\leqq \frac1{4n}$, $P(X=\frac kn)=\comb{n}{k}x^k(1-x)^k$ para $k=0,\ldots,n$. Dado $\delta>0$ tenemos $P(|X-x|\geqq\delta)\leqq\frac{\Var(X)}{\delta^2}\leqq \frac1{4n\delta^2}$ por Chebyshev. Dado $\epsilon>0$ hay $\delta>0$ con $|x-y|<\delta\then|f(x)-f(y)|<\epsilon$ si $x,y\in[0,1]$, luego $\E(|f(X)-f(x)|)=$$\E_{|X-x|<\delta}(|f(X)-f(x)|)+\E_{|X-x|\geqq\delta}(|f(X)-f(x)|)\leqq\epsilon+\frac{2\|f\|_\infty}{4n\delta^2}$. Luego $|\E(f(X))-\E(f(x))|\leqq \E(|f(X)-f(x)|)\to0$ uniformemente cuando $n\to\infty$; ahora $\E(f(X))=\sum_{k=0}^n f\left(\frac kn\right)\comb{n}{k}x^k(1-x)^k$, luego $\|p_n-f\|_\infty\to0$, y probamos Wierstrass. Ahora dada $w$, $\la f,g\ra=\int_0^1 wfg$, $f$ y $p_n^*$ tenemos $\|f-p_n^*\|\leqq\|f-p_n\|=\sqrt{\int_0^1w(f-p_n)^2}\leqq$$\sqrt{\int_0^1w\epsilon^2}=\epsilon\|1\|\to0$ si $n\to\infty$.∎

Integración numérica

Dado $[a,b]$ se buscan puntos $x_0,\ldots,x_n$ y pesos $A_0,\ldots,A_n$ para aproximar $I(f)=\int_a^bfw$, con $w>0$, por $Q(f)=\sum_{i=0}^n A_if(x_i)$. El grado de exactitud es el máximo $k$ tal que $I(p)=Q(p)$ para todo $p\in\PP_k$; por linealidad basta testear $x^k$. Si hay $M$ con $|Q(f)|\leqq M(b-a)\|f\|_\infty$ entonces si $f\in C^{k+1}[a,b]$ se tiene $|R(f)|=|I(f)-Q(f)|\leqq \frac{(1+M)(b-a)^{k+2}}{(k+1)!}\|f^{(k+1)}\|_{\infty}$. Si partimos $[a,b]$ en $a=x_0<\cdots<x_n=b$ y tenemos $|Q_i(f|_{[x_{i-1},x_i]})|\leqq M(x_i-x_{i-1})\|f\|_\infty$, las combinamos en $Q=\sum_{i=1}^nQ_i$ y tenemos $|R(f)|\leqq \frac{(1+M)(b-a)\max\{x_i-x_{i-1}\}^{k+1}}{(k+1)!}\|f^{(k+1)}\|_{\infty}$.

Newton-Côtes. Dados los $x_i$ tomamos el polinomio interpolador $p(x)=\sum_{i=0}^nf(x_i)\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}$ y se obtiene $A_i=\int_a^b\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}w(x)dx$. Llevar la regla de $[-1,1]$ a $[a,b]$ da $Q(f)=\frac{b-a}2\sum_{i=0}^nA_if\left(\frac{b-a}2x_i+\frac{a+b}2\right)$.

Trapecios. Es Newton-Côtes con $w=1$ y dos puntos: $\{-1,1\}$, es decir, $I(f)=f(-1)+f(1)$. (También está la versión abierta $\{-\frac12,\frac12\}$.) Error: si $f\in C^2[a,b]$, $R(f)=-\frac{(b-a)^3}{12}f''(\eta)$, con $\eta\in(a,b)$. Si partimos $[a,b]$ en $[x_i,x_i+h]$ con $h=\frac{b-a}n$ y $x_i=a+hi$ resulta $Q(f)=\left(\frac{f(a)}2+\sum_{i=1}^{n-1}f(x_i)+\frac{f(b)}2\right)h$ y $R(f)=-\frac{h^2}{12}(b-a)f''(\xi)$.

Prueba. $R(f)=\int_a^b(f-p)=\int_a^bf[a,b,x](x-a)(x-b)\,dx=\frac12f''(\eta)\int_a^b(x-a)(x-b)\,dx=$$-\frac{(b-a)^3}{12}f''(\eta)$ (usamos valor medio integral).∎

Simpson. Es Newton-Côtes con $w=1$ y tres puntos: $\{-1,0,1\}$, es decir, $I(f)=\frac{f(-1)+4f(0)+f(1)}3$. (También está la versión abierta $\{-\frac12,0,\frac12\}$ que da $\frac{4f(-\frac12)-2f(0)+4f(\frac12)}{3}$.) Error: si $f\in C^4[a,b]$, $R(f)=-\frac1{90}\left(\frac{b-a}2\right)^5f^{(4)}(\eta)$, con $\eta\in(a,b)$. Si partimos $[a,b]$ en $[x_i,x_i+2h]$ con $h=\frac{b-a}{2n}$ y $x_i=a+2hi$ resulta $Q(f)=\left(f(a)+2\sum_{i=1}^{n-1}f(x_i)+4\sum_{i=0}^{n-1}f(x_i+h)+f(b)\right)\frac{h}{3}$ y $R(f)=-\frac{h^4}{180}(b-a)f^{(4)}(\xi)$.

Prueba. La idea es ver el problema en $[-h,h]$ y poner $e(t)=\int_{-t}^t f-t\left(\frac{f(-t)+4f(0)+f(t)}{3}\right)$, calcular $e'''(t)=\frac{f'''(-t)-f'''(t)}3t=-\frac23t^2f^{(4)}(\xi)$, luego obtener $e''$, $e'$ y $e$ integrando y usando valor medio integral, lo que da $e(h)=-\frac{h^5}{90}f^{(4)}(\eta)$.∎

Cuadratura gaussiana. Sea $\la f, g\ra=\int_a^b fgw$; sean $p_n$ los polinomios ortonormales; $p_n$ tiene $n$ raíces distintas $x_1,\ldots,x_n$ en $(a,b)$. Newton-Côtes sobre $x_i$ da la cuadratura gaussiana $Q(f)=\sum_{i=1}^n A_if(x_i)$; $Q$ tiene grado de exactitud $2n-1$; el grado de exactitud es máximo para $n$ puntos y los $x_i$ son los únicos que lo alcazan; además $A_i>0$. Error: si $f\in C^{2n}[a,b]$, $R(f)=\frac{f^{(2n)}(\xi)}{(2n)!}\int_a^b q_n^2w$, donde $q_n=\prod_{i=1}^n(x-x_i)$.

Prueba. Primero, vemos que $p_n$ tiene raíz en $(a,b)$ (viendo $\la p_n,1\ra=0$); luego que no tiene raíces distintas (si no $p_n=(x-x_0)^2q$ y vemos $\la p_n,q\ra=0$); luego que tiene $n$ raíces (si no $p_n=qr$ con $\deg r>0$ sin raíces y vemos $\la p_n,q\ra=0$). Segundo, $Q$ tiene grado de exactitud $2n-1$: si $f\in\PP_{2n-1}$, $f=p_nq+r$ con $q,r\in\PP_{n-1}$, luego $r$ es el interpolador de $f$ en los $x_i$, luego $R(f)=\int_a^b(f-r)w=\la p_n,q\ra=0$; si $f=\prod_{i=1}^n(x-x_i)^2$, $r=0$ y $R(f)>0$. Tercero, si $x_0',\ldots,x_n'$ son otros puntos y $R(f)=0$ para todo $f\in\PP_{2n-1}$, si $p=\prod_{i=1}^n(x-x_i')$ y $q\in\PP_{n-1}$ tengo $0=R(pq)=\la p,q\ra$, luego $p$ es ortogonal a $\PP_{n-1}$ y tiene que ser $\lambda p_n$. Cuarto, si $f_i=\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}$, $0<\int_a^b f^2w=\sum_{i=1}^n A_if(x_i)^2=A_i$. Quinto, sea $f\in C^{2n}[a,b]$ y $p\in\PP_{2n-1}$ con $p(x_i)=f(x_i)$ y $p'(x_i)=f'(x_i)$; tengo $R(f)=R(f-p+p)=R(f-p)=\int_a^b(f-p)w=\int_a^b\frac{f^{(2n)}(\eta)}{(2n)!}q_n^2w=\frac{f^{(2n)}(\xi)}{(2n)!}\int_a^b q_n^2w$.∎

Fórmula de Rodrígues. Los polinomios ortogonales para $\la f, g\ra=\int_a^b fg$ se llaman polinomios de Legendre $\ell_n$. Rodrígues: $\ell_n=\frac{n!}{(2n)!}\frac{d^n}{dx^n}(x^2-1)^n$ y $\int_{-1}^1\ell_n^2=\frac{2^{2n+1}(n!)^4}{(2n+1)(2n)!^2}$. Error de la cuadratura gausiana en $[a,b]$: $R(f)=\frac{(n!)^4(b-a)^{2n+1}}{(2n+1)(2n)!^3}f^{(2n)}(\xi)$.

Prueba. Si $p\in \PP_{n-1}$ se ve que $\int_{-1}^1 \ell_np=(-1)^n\int_{-1}^1 \frac{n!}{(2n)!}(x^2-1)^np^{(n)}=0$ por integración por partes. Si $q=(x^2-1)^n$, $\int_{-1}^1 (\partial^n q)^2=(-1)^n\int_{-1}^1 q\partial^{2n}q=(2n)!\int_{-1}^1(1-x^2)^n$, que es, con $x=\sen y$, $\int_{-\frac\pi2}^{\frac\pi2}\cos^{2n+1}=\int_{-\frac\pi2}^{\frac\pi2}\cos^{2n}\sen'=\int_{-\frac\pi2}^{\frac\pi2}2n\cos^{2n-1}\sen^2=\int_{-\frac\pi2}^{\frac\pi2}2n\cos^{2n-1}(1-\cos^2)$. Entonces $\int_{-\frac\pi2}^{\frac\pi2}\cos^{2n+1}=\frac{2n}{2n+1}\int_{-\frac\pi2}^{\frac\pi2}\cos^{2n-1}$ y $\int_{-\frac\pi2}^{\frac\pi2}\cos^{2n+1}=\frac{(2^nn!)^2}{(2n+1)!}\int_{-1\frac\pi2}^{\frac\pi2}\cos=\frac{2^{2n+1}(n!)^2}{(2n+1)!}$ y se obtiene $\int_{-1}^1\ell_n^2=\frac{2^{2n+1}(n!)^4}{(2n+1)(2n)!^2}$.∎

Ecuaciones diferenciales ordinarias

Método de Taylor. Dada la ecuación $x’(t)=f(t,x(t))$, $x(t_0)=x_0$ queremos encontrar $x(T).$ Tomamos $n$, $h=\frac{T-t_0}n$, $t_i=t_0+hi$. Tenemos $x(t_{i+1})=x(t_i+h)=\sum_{k=0}^p \frac{h^k}{k!}x^{(k)}(t_i)+$ $\frac{h^{p+1}}{(p+1)!}x^{(p+1)}(\xi)$. Dada una aproximación $x_i\approx x(t_i)$ podemos aproximar $x'(t_i)=f(t_i,x(t_i))\approx f(t_i,x_i)$, $x''(t_i)=f_t(t_i,x(t_i))+f_x(t_i,x(t_i))x'(t_i)\approx$ $f_t(t_i,x_i)+f_x(t_i,x_i)f(t_i,x_i)$, etc, y obtenemos una aproximación de $x_{i+1}\approx x(t_{i+1})$, hasta $x_n\approx x(T)$.

Métodos de Runge-Kutta de orden 2. Dado $\alpha\in[0,1]$, $x_{i+1}=x_i+h\Phi(t_i,x_i,h)$, con $\Phi(t,x,h)=(1-\frac1{2\alpha}))f(t,x)+\frac1{2\alpha}f(t+\alpha h, x+\alpha hf(t,x))$. Vale por Taylor que $x(t+h)-x(t)-h\Phi(t,x,h)=O(h^3)$.

Error. Dado $n$, $h=\frac{T-t_0}n$, $x_{i+1}=x_i+h\Phi(t_i,x_i,h)$ con $K$ tal que $|\Phi(t,x,h)-\Phi(t,y,h)|\leqq K|x-y|$ para todos $x,y\in\R$, $t\in[t_0,T]$, y $\tau$ tal que $\left|\frac{x(t_{i+1})-x(t_i)}h-\Phi(t_i,x(t_i),h)\right|\leqq\tau$, se tiene $|x(T)-x_n|\leqq \frac{\tau}{K}(e^{K(T_0-t)}-1)$.

Prueba. Si $e_i=|x(t_i)-x_i|$ tenemos $e_1\leqq h\tau$; $x(t_{i+1})=x(t_i)+h\Phi(t_i,x(t_i),h)+h\tau_i$, con $|\tau_i|\leqq\tau$; y $x_{i+1}=x_i+h\Phi(t_i,x_i,h)$; luego $x(t_{i+1})-x_{i+1}=x(t_i)-x_i+$ $h(\Phi(t_i,x(t_i),h)-\Phi(t_i,x_i,h))+h\tau_i$ y $e_{i+1}\leqq e_i+hKe_i+h\tau=(1+Kh)e_i+\tau h$; luego $e_n\leqq \frac{\tau}K((1+Kh)^{n}-1)\leqq \frac{\tau}{K}(e^{K(T_0-t)}-1)$.∎

Métodos multipaso lineales. Consisten en aproximar $x_{n+k}$ dados $x_{n},\ldots,x_{n+k-1}$ por una fórmula $\sum_{i=0}^k \alpha_ix_{n+i}=h\sum_{i=0}^k \beta_if(t_{n+i},x_{n+i})$. Decimos que es de orden $p$ si $\tau_th=\sum_{i=0}^k \alpha_ix(t+hi)-h\sum_{i=0}^k \beta_ix'(t+hi)=O(h^{p+1})$; equivale a que $\sum_{i=0}^k\alpha_i=0$ y $\sum_{i=0}^k i^q\alpha_i=q\sum_{i=0}^k i^{q-1}\beta_i$ para $q=1,\ldots,p$. Sea $p(x)=\sum_{i=0}^k \alpha_ix^i$; si para toda raíz simple $r$ de $p$ vale $|r|\leqq 1$ y para toda raíz múltiple vale $|r|<1$, y además $\lim_{h\to0}\tau_t=0$, entonces el método converge.