La solución numérica elemental con Octave

LA SOLUCIÓN NUMÉRICA ELEMENTAL CON OCTAVE

Página 1 de 115
MATIAUDA MARIO EUGENIO
YACHECEN CARLOS GUSTAVO
 EDITORIAL UNIVERSITARIA DE MISIONES

San Luis 1870

Posadas - Misiones – Tel-Fax: (03752) 428601
Correos electrónicos:
edunam-admini@[Link]
edunam-direccion@[Link]
edunam-produccion@[Link]
edunam-ventas@[Link]

Colección: Cuadernos de Cátedra

Coordinación de la edición: Claudio Oscar Zalazar

ISBN 978-950-579-190-3
Impreso en Argentina
©Editorial Universitaria

Matiauda, Mario Eugenio

La solución numérica elemental con octaves. - 1a ed. - Posadas : EDUNAM - Editorial
Universitaria de la Universidad Nacional de Misiones, 2011.
Internet.
ISBN 978-950-579-190-3
1. Matemática . 2. Solución Numérica. I. Título.
CDD 510
Fecha de catalogación: 10/03/2011
 La solución numérica elemental con Octave

INDICE

UNIDAD 1 : RESOLVER f(x) =0

1.1 Método de la bisección......................................................................................... 6
1.2 Método del punto fijo ............................................................................................ 6
1.3 Método de Newton ............................................................................................... 7
1.4 Método de la secante ........................................................................................... 7
1.5 Método de la posición falsa .................................................................................. 8
1.6 Método de Steffensen .......................................................................................... 8
1.7 Polinomios. Raíces. Métodos ............................................................................... 8
1.7.1 Método de Horner para ubicar raíces de un polinomio P de grado n............. 9
1.7.2 Método de Muller ........................................................................................... 9
UNIDAD 2 : INTERPOLACIÓN
2.1 Interpolación. Empleo de polinomios.................................................................. 11
2.2 Polinomio de Lagrange ...................................................................................... 11
2.3 Polinomio de Hermite ......................................................................................... 12
2.4 Método de Neville .............................................................................................. 13
2.5 Diferencias divididas .......................................................................................... 14
2.6 Polinomios osculadores. Polinomios de Hermite ............................................... 16
2.7 Adaptador cúbico ............................................................................................... 16
UNIDAD 3 : DIFERENCIACIÓN E INTEGRACIÓN NUMÉRIA
3.1 Diferenciación numérica. Formulas de aproximación ......................................... 17
3.2 Diferenciación numérica. Técnica de Richardson .............................................. 18
3.3 Integración numérica .......................................................................................... 18
3.3.1 Fórmula de Newton...................................................................................... 18
3.3.2 Regla compuesta del trapezio ..................................................................... 19
3.3.3 Regla compuesta de Simpson ..................................................................... 19
3.3.4 Integración de Romberg .............................................................................. 20
3.3.5 Regla trapezoidal recursiva ......................................................................... 21
3.4 Fórmulas Gaussianas ........................................................................................ 21
3.4.1 Fórmula de Gauss-Legendre ....................................................................... 21
3.4.2 Fórmula de Gauss-Laguerré ........................................................................ 22
3.4.3 Fórmula de Gauss-Chebyshev .................................................................... 22
Página 2 de 115
 La solución numérica elemental con Octave

3.5 Cuadratura adaptativa utilizando la regla de Simpson’s..................................... 22

3.6 Regla de Simpson’s ........................................................................................... 23
3.7 Cuadratura de Gauss-Legendre ......................................................................... 23
UNIDAD 4 : ECUACIONES DIFERENCIALES ORDINARIAS
4.1 Métodos para resolver P.V.I. ............................................................................. 24
4.1.1 Método de Euler .......................................................................................... 24
4.1.2 Serie de Taylor ............................................................................................ 24
4.1.3 Método de Runge-Kutta ............................................................................... 25
4.1.4 Métodos multipasos ..................................................................................... 25
[Link] Método de Adams-Bashforth .......................................................... 26
[Link] Método de Adams-Moulton ............................................................. 27
4.1.5 Método de Milne .......................................................................................... 28
4.1.6 Predictor-Corrector ...................................................................................... 28
4.1.7 Uso de la extrapolación ............................................................................... 28
UNIDAD 5 : SISTEMAS LINEALES
5.1 Representación de un sistema lineal. Sustitución hacia atrás ............................ 29
5.2 Técnica de Gauss Jordan .................................................................................. 30
5.3 Pivoteo ............................................................................................................... 30
5.4 Factorización matricial........................................................................................ 31
5.4.1 Factorización LU .......................................................................................... 31
5.5 Técnicas de repetición para sistemas lineales .................................................. 34
5.5.1 Técnica de Jacobi ........................................................................................ 35
5.5.2 Técnica de Gauss-Seidel ............................................................................. 35
5.6 Técnica de relajación ......................................................................................... 36
5.7 Gradiente conjugado .......................................................................................... 38

UNIDAD 6 : APROXIMACIÓN
6.1 Polinomios de Chebyschev ................................................................................ 39
6.2 Aproximación mínimo-máximo (o de Chebyshev) .............................................. 39
6.3 Aproximación por valores propios ...................................................................... 39
6.4 Técnicas de aproximación.................................................................................. 39
6.4.1 Método de potencias.................................................................................... 39
6.4.2 Método de la potencia inversa ..................................................................... 40

Página 3 de 115
 La solución numérica elemental con Octave

6.4.3 Método simétrico de potencias .................................................................... 41

6.4.4 Método clásico de Jacobi ............................................................................. 41
6.4.5 Técnicas de deflación .................................................................................. 42
6.4.6 Técnicas de Householder ............................................................................ 42
6.4.7 Algoritmo QR ............................................................................................... 42
6.5 Ortogonalización de Gram-Schmidt ................................................................... 42

UNIDAD 7 : SISTEMAS DE ECUACIONES NO LINEALES

7.1 Método de Newton ............................................................................................. 43
7.2 Técnicas cuasi-Newton ...................................................................................... 43
7.3 Técnicas de mayor pendiente ............................................................................ 44
7.4 Problema de valor de frontera para EDO’s y en derivadas parciales ................. 44
7.5 Técnica de disparo para el problema lineal ........................................................ 44
7.6 Técnica de disparo para el problema no lineal ................................................... 45
7.7 Diferencias finitas para problemas lineales ........................................................ 45

UNIDAD 8 : ECUACIONES EN DERIVADAS PARCIALES

8.1 Ecuación general ............................................................................................... 47
8.2 Elíptica o de Poisson .......................................................................................... 47
8.3 Ecuación parabólica ........................................................................................... 48
8.4 Ecuación hiperbólica .......................................................................................... 49

Página 4 de 115
 La solución numérica elemental con Octave

Prólogo

En diversas ocasiones, situaciones del mundo físico se desean representar

matemáticamente pero, esos problemas no encuentran una resolución analítica
exacta o deben ser encarados desde la óptica Numérica.
Varias de esas situaciones se incluyeron en la Solucion Numérica Elemental del autor,
conteniendo la alternativa de búsqueda de resultados a través de la computadora.
En esta presentación se aborda la misma línea, simplemente iniciando al uso de
Octave, sofá libre y de permanente actualización y con amplia Potencialidad de
emulación de software bajo licencia (caso MatLab).
Por ello no se abunda en desarrollo de los métodos sino directamente su ejecución
bajo Octave.

Página 5 de 115
 La solución numérica elemental con Octave

RESOLVER f(x) =0

SOLUCIONES DE ECUACIONES UNIVARIABLES

Para una función f ( x)  0 , el problema más simple de la aproximación numérica será

hallar la raíz x que anula f .

1.1 METODO DE LA BISECCIÓN

Se parte de f definida en I real con extremos a y b pertenecientes a I, bajo el

supuesto de que f (a). f (b)0 .

ai  bi
pi  Con pi = punto medio
2
ab ab
Se toma el punto medio si f ( )0 ya hemos encontrado la raíz
2 2
ab ab ab
x . En caso contrario, si f ( ) f (b)  0 entonces hacemos a  y
2 2 2
volvemos a subdividir el nuevo intervalo [a, b]. Si, por el contrario,
ab ab
f (a) f ( )  0 entonces hacemos b  y volvemos a empezar. Las
2 2
sucesivas subdivisiones del intervalo [a, b]. Van aproximando la raíz.

ba
Criterio de detención | pn  pn1 |  n 1
2n
Algoritmo bisect

1.2 METODO DEL PUNTO FIJO

El problema de hallar raíces de h(q)  0 con q =Punto fijo

h( x)  x  f ( x) o h( x)  x   f ( x) ,   R

Si la función h(x) tiene un punto fijo en p , entonces la función definida por

f ( x)  x  h( x) tiene un cero en p .
La forma de punto fijo es más fácil de analizar que la búsqueda de raíces.

Página 6 de 115
 La solución numérica elemental con Octave

Si h  C  a, b y h( x)  a, b para toda x   a, b , h tiene un punto fijo en

 a, b  .
Si también g '( x) a, b con una constante m1 , de modo que
existe en

h '( x)  m1 x  a, b

Que garantiza que el punto fijo es único en  a, b  .

Algoritmo fixed_point

1.3 METODO DE NEWTON

Partiendo de un p0 de arranque (requiere la elección adecuada de p0 ), se genera
la sucesión  pn  definida a través de:
g( p )
pn  pn 1  ' n 1 para n  1
g ( pn 1 )
Es una técnica de iteración funcional, donde las aproximaciones se obtienen usando
tangentes sucesivas.

Criterio de detención pn  pN 1   0 o
pN  pN 1
 pN  0
pN
Algoritmo newton

1.4 METODO DE LA SECANTE

Se parte de las aproximaciones de arranque x0 y x1 , donde las posteriores

aproximaciones se obtienen usando secantes sucesivas; evitando así el problema de
conocer la derivada de la función como en la técnica de Newton, así sustituye la
derivada por la relación de las diferencias
.

g ( xn 1 )( xn 1  xn 2 )
xn  xn 1 
g ( xn 1 )  g ( xn 2 )

Algoritmo secant

Página 7 de 115
 La solución numérica elemental con Octave

1.5 METODO DE LA POSICIÓN FALSA

Genera aproximaciones del mismo modo que el método de la secante, pero incorpora
el acorralamiento de la raíz, asegurando que la misma quede entre dos iteraciones
sucesivas.

f (b)  f (a)
y  f (a )  ( x  a) Como x1 es el valor de x para y=0, se tiene
ba
f (a) af (b)  bf (a) x1  a x b
x1  a  ( x  a)  si  ó 1 
f (b)  f (a) f (b)  f (a) x1 x1

para una  pre-establecida, entonces x1 es la raíz buscada. De lo contrario, se

calcula f ( x1 ) y elegir entre a y b aquella cuya f tenga signo opuesto de f ( x1 ). Con
x1 *(en)ese punto calculado x2 mediante la fórmula de la secante. El proceso
iterativo debe continuar hasta obtener la raíz con la precisión pre-establecida.

Algoritmo regula

1.6 METODO DE STEFFENSEN

Genera la secuencia:

 
p0( 0) , p1( 0)  g ( p0( 0) ), p2( 0)  g ( p1( 0) ), p0(1)  2 ( p0( 0) ), p1(1)  g ( p0(1) ),.......
(pn ) 2 
donde  indica que se usa la ecuación p  pn  2 , para n  0.
2

 pn
la cual genera cada tercer término; los demás usan la iteración del punto fijo en el
término anterior.
Para encontrar rápidamente una solución de la ecuación de punto fijo x = g(x) dada
una aproximación inicial P0; donde se supone que tanto g(x) y g'(x) son continuas,
|g’(x)|<1, y que la iteración de punto fijo ordinario converge lentamente (lineal) para p

Algoritmo steff

1.7 POLINOMIOS. RAÍCES. MÉTODOS

P  x   an x n  an1 x n1  ...  a1x  a0 Polinomio de grado n

Con ai  coeficientes de P (son constantes)

Página 8 de 115
 La solución numérica elemental con Octave

1.7.1 Método de Horner

Sea bn  an
bnk  xbnk 1  ank , k  1,2,..., n . (Fórmula de recurrencia)

Basándonos en la fórmula de recurrencia dada, sustituyendo x por z, si z es una raíz

de P(x), podemos escribir que:
P( x)  ( x  z)Q( x)
P( x)
Donde Q( x)  bn x n1  bn1 x n2  ...  b1  en efecto
xz
n 1 n2
(bn x  bn 1 x  ...  b1 )( x  z )
 bn x n  (bn 1  zbn ) x n 1  ...
 (b1  zb2 ) x  (b0  zb1 )
 an x n  an 1 x n 1  ...  a1 x  a0  P( x)
Concluimos que cualquier raíz de Q(x) es también una raíz de P(x), esto nos permite
trabajar con un polinomio de grado n-1,o sea, con Q(x), para calcular las raíces
subsecuentes de P(x), este proceso recibe el nombre de deflación, con este se evita
que un mismo cero sea calculado varias veces.

Algoritmo horner

Algoritmo newtonhorner

1.7.2 Método de Muller

Es una extensión del método de la secante, utiliza tres aproximaciones iniciales: x0 y

x1 , x2 .

Considerando el polinomio cuadrático:

f 2 ( x)  a( x  x2 ) 2  b( x  x2 )  c

Con

Página 9 de 115
 La solución numérica elemental con Octave

c  f ( x2 )
( x0  x2 )[ f ( x1 )  f ( x2 )]  ( x1  x2 ) 2 [ f ( x0 )  f ( x2 )]
b
( x0  x2 )( x1  x2 )( x0  x1 )
( x1  x2 )[ f ( x0 )  f ( x2 )]  ( x0  x2 )[ f ( x1 )  f ( x2 )]
a
( x0  x2 )( x1  x2 )( x0  x1 )

Hallando la raíz, se implementa la solución convencional, pero debido al error de

redondeo potencial, se usará una formulación alternativa:
 2c  2c
x3  x 2  despejando x3  x 2 
b  b 2  4ac b  b 2  4ac

La gran ventaja de este método es que se pueden localizar tanto las raíces reales
como las imaginarias.

La elección de las fórmulas anteriores equivale a aproximar f (x) por la parábola

que pasa por los puntos ( xn3 , f ( xn3 ), ( xn2 , f ( xn2 ), ( xn1 , f ( xn1 ), y
calcular posteriormente las derivadas de dicha parábola.

Algoritmo muller

Algoritmo muller2

Página 10 de 115
 La solución numérica elemental con Octave

INTERPOLACION

2.1 INTERPOLACION. EMPLEO DE POLINOMIOS

El problema general de la interpolación de funciones consiste en, a partir del

conocimiento del valor de una función (y eventualmente de sus derivadas) en un
conjunto finito de puntos, aproximar el valor de la función fuera de ese conjunto
finito de puntos.

Según el teorema de aproximación

Con  > 0,
f ( x)  P( x)  x   a, b
El uso de los polinomios tiene consigo ventajas como su derivabilidad e integrabilidad,
dando también polinomio.
El planteo del problema será hallar un polinomio interpolante que brinde una
aproximación precisa en todo el intervalo, de allí la no adecuación de los polinomios
de Taylor que concentran la información alrededor de un punto, dígase x0 .

2.2 POLINOMIO DE LAGRANGE

f ( xk )  P( xk ) k  0,1, 2,....n
n
P( x)  f ( x0 ) Ln,0 ( x)  .....  f ( xn ) Ln,n ( x)   f ( xk )Ln,k ( x)
k 0
Con

Ln,k ( x) 
 x  x0  x  x1  ...  x  xk 1  x  xk 1  ...  x  xn 
 xk  x0  xk  x1  ...  xk  xk 1  xk  xk 1  ...  xk  xn 
Ln,k ( x)  
n
 x  xi  k  0,1, 2,....n
i  0  xk  xi 
ik

Uno de los inconvenientes que presenta la interpolación de Lagrange, es que el grado

del polinomio necesario, no se conoce, hasta la determinación de los cálculos
necesarios.

Algoritmo lagrange

Página 11 de 115
 La solución numérica elemental con Octave

2.3 POLINOMIO DE HERMITE

En ocasiones, resulta de interés interpolar no sólo el valor de la función en

ciertos puntos {xi }i=0,..,N , sino también el valor de sus derivadas. Un
ejemplo clásico de ello es el desarrollo de Taylor de una función en un punto a.
En este caso, aproximamos f (x) por un polinomio de grado N , PN (x) tal que f (x)
y PN (x) poseen las mismas derivadas en el punto a desde el orden 0 hasta el
orden N.

f ' (a) f ( N ) (a)

PN ( x)  f (a)  ( x  a)  ...  ( x  a) ( N )
1! N!
El error de interpolación viene dado por la fórmula

f ( N 1) ( )
f ( x)  PN ( x)  ( x  a) ( N 1)
( N  1)!
donde ξ es un valor intermedio entre x y a. En el caso general, donde
buscamos un polinomio P (x) tal que él y todas sus derivadas hasta un cierto
orden M coincidan con una función f (x) en los puntos { x i }i=0,..,N, se utilizan
los denominados polinomios base de Hermite Hi,j (x), que son polinomios de
grado menor o igual que
(N + 1)(M + 1) – 1 dados por las siguientes condiciones:

l H i, j 
1 si l  j y k  i
( xk )  
x l 
0 si l  j o k  i
A partir de los polinomios base de Hermite, el polinomio interpolador de Hermite se
define como:

N M
j f
P( x)    ( xi ) H i , j ( x)
i 0 j 0 x j

Algoritmo hermite

Página 12 de 115
 La solución numérica elemental con Octave

2.4 MÉTODO DE NEVILLE

Si N i , j con 0  i  j representa el polinomio interpolante de grado j en los

 j  1 números xi  j , , xi  j 1 ,..., xi 1, xi
Ni , j  Pi  j ,i  j 1,...,i 1,i

Se dispondrá de un diagrama para los polinomios interpolantes:

x0 P0  N 0,0


x1 P1  N1,0 P0,1  N1,1
x2 P2  N 2,0 P1,2  N 2,1 P0,1,2  N 2,2
i más grande
x3 P3  N30 P2,3  N3,1 P1,2,3  N3,2 P0,1,2,3  N3,3
x4 P4  N 4,0 P3,4  N 4,1 P2,3,4  N 4,2 P1,2,3,4  N 4,3 P0,1,2,3,4  N4,4
 Desplazarse hacia la derecha aumenta el grado del polinomio

( x  xi  j ) N i , j 1  ( x  xi ) N i 1, j 1
Ni, j 
xi  xi  j

Algoritmo neville1

Algoritmo neville2

Página 13 de 115
 La solución numérica elemental con Octave

2.5 DIFERENCIAS DIVIDIDAS

Para calcular las diferencias divididas de una función f ( x) en x0 , x1 , x2 ,..., xn

construimos una tabla de diferencias divididas

La primera columna contiene los puntos xk , k  0,1,..., n

La segunda contiene los valores de la función en xk , k  0,1,..., n
En las columnas 3, 4, 5,…., están las diferencias divididas de orden 1, 2, 3,…. Cada
una de estas es, una fracción cuyo numerador es siempre una diferencia entre dos
diferencias divididas consecutivas, y de orden inmediatamente inferior y, cuyo
denominador es la diferencia entre los dos extremos de los puntos involucrados.

f1  f 0
f ( x0 , x1 ) 
x1  x0
f ( x1 ,..., xn )  f ( x0 ,..., xn1 )
Para orden n : f ( x0 , x1 ,..., xn ) 
xn  x0

Fórmulas de diferencias divididas interpolantes de Newton

n
P( x)  f  x0    f  x0 ,..., xk   x  x0  ... x  xk 1  para cada k  0,1,..., n
k 1

Fórmula de Newton de diferencias divididas hacia adelante

Página 14 de 115
 La solución numérica elemental con Octave

Una simplificación de la fórmula de interpolación de diferencias divididas se obtiene

considerando igual espaciado entre los símbolo, x0 , x1 ,..., xn
Denotando h  xi 1  xi para i  0,1,..., n  1 con x  x0  sh .
Será x  xi   s  i  h , convirtiendo a
n
P( x)   s  s  1 ...  s  k  1 h k f  x0 , x1 ,..., xk 
k o
Fórmula de Newton de diferencias divididas hacia atrás

Si los nodos tienen espacios iguales con x  xn  sh y x  xi  (s  n  i)h

entonces

Pn ( x)  Pn ( xn  sh)  f [ xn ]  sh f [ xn , xn1 ]  .....

s( s  1)h 2 f [ xn , xn1 , xn2 ]  s( s  1) ( s  n  1)h n f [ xn ,..., x0 ].
x  xn
Con s 
h
Diferencias centradas de Stirling

Si se busca aproximar un valor de x que se ubica cerca del centro de la tabla, no son
adecuadas las expresiones de Newton.
Se emplean las fórmulas de diferencias centradas, según sea el caso, y dada la
amplitud de expresiones se mencionan la de Stirling.

P( x)  P2 m 1  x   f  x0  
sh
2
 f  x1 , x0   f  x0 , x1   s 2h 2 f  x1 , x1 , x0 

s  s 2  1 h3

2
 f x 1 , x0 , x1 , x2   f  x2 , x1 , x0 , x1   ...

 
 s 2  s 2  1 s 2  4  ... s 2   m  1 h 2 m f  x m ,....xm 
2

s  s 2  1 ...  s 2  m 2  h 2 m1

2
 f x m ,....xm 1   f  x m1 ,....xm 
para n=2m+1 impar, y si n=2m es par se elimina el último término de la misma
expresión.

Algoritmo divdiff1

Algoritmo newpoly

Página 15 de 115
 La solución numérica elemental con Octave

2.6 POLINOMIOS OSCULADORES. POLINOMIOS DE HERMITE

Polinomio con grado  2n  1 y oscilación de primer orden, expresado por

n n
P2 n1 ( x)   f ( x j ) Pn, j ( x)   f '( x j ) Pn, j ( x)
j 0 j 0

Siendo Pn, j ( x)  1  2  x  x j  L 'n, j ( x j )  L2n, j ( x j )

Pn, j ( x)   x  x j  L2n, j ( x j )
Con Ln , j el j-ésimo polinomio de coeficientes de Lagrange(grado n).

Dado lo laborioso de determinar los polinomios de Lagrange con sus derivadas,

alternativamente para aproximar por Hermite se hace uso de la fórmula de diferencia
para los polinomios de Lagrange.

2.7 ADAPTADOR CÚBICO

El empleo de polinomios cúbicos garantiza la derivabilidad en el intervalo y posee

derivadas segundas continuas, sin suponer que coincidan derivadas del polinomio
interpolantes con las de la función.

h j 1c j 1  2  h j 1  h j  c j  h j c j 1   a j 1  a j    a j  a j 1 
3 3
hj h j 1
Para j  0,1,..., n  1
h = al paso entre nodos para j  0,1,..., n  1 y an  f  xn 
Siendo:
c j 1  c j  3d j h j
h2j
 a j 1  a j    2c j  c j 1 
hj
 2c  c j 1  con b j 
1
a j 1  a j  b j h j  j
3 hj 3
b j 1  b j  h j  c j  c j 1 

el conjunto c 
n
j j 0
serán las incógnitas.

Al conocer los c j , se podrán hallar los coeficientes b j y d j , para poder generar los

A 
n 1
polinomios cúbicos j j 0
.
Algoritmo linear_spline
Algoritmo spline_eval

Página 16 de 115
 La solución numérica elemental con Octave

DIFERENCIACIÓN E INTEGRACIÓN NUMERICAS

3.1 DIFERENCIACIÓN NUMERICA. FÓRMULAS DE APROXIMACIÓN

Para un intervalo con (n+1) puntos diferentes  x0 , x1 ,..., xn  con f  C

n 1
sobre dicho
intervalo
f n 1 ( ( x)) n
  xk  x j 
n
f '( x)   f ( x j )L ' j ( x) 
j 0 (n  1)! j 0
j 0

que combina linealmente (n+1) valores de f ( x) , para j  0,1,..., n , conocida como

fórmula de (n+1) puntos.

Fórmulas de tres puntos

1 h2 3
f ( x)   3 f ( x0 )  4 f ( x0  h)  f ( x0  2 h)   f ( 0 )
2h 3
Con  0 entre x0 y x1  2h, y
1 h2 (3)
f ( x)   f ( x0  h)  f ( x0  h)  f ( 1 )
2h 6
Con 1 entre ( x0  h) y ( x0  h)

Fórmulas de cinco puntos

1 h4
f '( x0 )   f ( x0  2h)  8 f ( x0  h)  8 f ( x0  h)  f ( x0  2h)  f (5) ( )
12h 30

1 h4
f '( x0 )  25 f ( x0 )  48 f ( x0  h)  36 f ( x0  2h)  16 f ( x0  3h)  3 f ( x0  4h)  f (5) ( )
12h 5

Con 0 entre x0 y x0  4h

La diferenciación basada en N+1 puntos para aproximar f ' ( x0 ) numéricamente

mediante la construcción del polinomio de Newton de enésimo grado.

P( x)  a0  a1 ( x  x0 )  a2 ( x  x0 )( x  x1 )
 a3 ( x  x0 )( x  x1 )( x  x2 )  ...  aN ( x  x0 )...( x  xN 1 )

Y usando f ' ( x0 )  P' ( x0 ) como la aproximación esperada. El método debe ser

utilizado en x0. Los puntos se pueden cambiar xk , x0 ,..., xk 1 , xk 1 ,..., xN  para
calcular f ' ( xk )  P' ( xk ) .

Algoritmo diffnew
Página 17 de 115
 La solución numérica elemental con Octave

La diferenciación usando límite para aproximar f ' ( x0 ) , generando la secuencia

f ( x  10 k h)  f ( x  10 k h)
f ' ( x)  Dk  k  0,..., n
2(10k h)
Ya que | Dn1  Dn || Dn  Dn1 | o | Dn  Dn1 | tolerancia , que es un intento de
encontrar la mejor aproximación f ' ( x)  Dn .

Algoritmo difflim

3.2 DIFERENCIACIÓN NUMERICA. TÉCNICA DE RICHARDSON

La expresión general será:

m 1
R  R(h)   M j h j O(h m ) para c / j  2,3,..., m
j 1
Llamando R(h) la expresión que aproxima el valor no conocido R, con un error de
truncado O(h) , para ciertas constantes M1 , M 2 , M 3 ...:
Existe una aproximación para cada j del tipo:

Rh1  h / 2   R j 1  h 
R(h)  R j 1 
2 j 1  1
La técnica (extrapolante) es aplicable siempre y cuando el error de truncado para una
fórmula sea del tipo:
m 1

 M h
j 1
j
j
O(h m )

Para constantes M j y que se verifique 1 1 ...1  m

Algoritmo diffext

3.3 INTEGRACIÓN NUMERICA

b b

Si P(x) aproxima a f(x), se tendrá:  P( x)dx  f ( x)dx

a a

3.3.1 Fórmula de Newton

x1
h
 P( x)dx  2  f ( x )  f ( x ) 
x0
0 1

x2
h
 P( x)dx  3  f ( x )  4 f ( x )  f ( x ) 
x0
0 1 2

Página 18 de 115
 La solución numérica elemental con Octave

x3
 3h 
x P ( x ) dx     f ( x0 )  3 f ( x1 )  3 f ( x2 )  f ( x3 ) 
0
 8 

3.3.2 Regla Compuesta del Trapecio

Usa segmentos de recta como aproximación a f ( x) .

xn
h
 f ( x)dx  2  f ( x )  2 f ( x )  ...  2 f ( x
x0
0 1 n 1 )  f ( xn ) 

Incluye nodos del intervalo abierto  x0 , xn 

O también es posible expresarlo cómo:

M 1
h
f ( x)dx  ( f (a)  f (b))  h  f ( xk )
b
a 2 k 1

por muestreo f ( x) en los M  1 puntos equiespaci ados

xk  a  kh, para k  0,1,2...., M . Siendo x0  a y xM  b

Algoritmo traprl

3.3.3 Regla compuesta de Simpson

xn
h
 f ( x)dx  3  f ( x )  4 f ( x )  2 f ( x )  4 f ( x )  ...  2 f ( x
x0
0 1 2 3 n2 )  4 f ( xn1 )  f ( xn ) 

Incluye los puntos extremos del intervalo cerrado  x0 , xn  .

O también es posible expresarlo cómo:

h 2h M 1 4h M
 2k 3 
b
a
f ( x)dx 
3
( f (a)  f (b)) 
3 k 1
f ( x ) 
k 1
f ( x2 k 1 )

por muestreo f ( x) en los 2M  1 puntos equiespaci ados

xk  a  kh, para k  0,1,2....,2M . Siendo x0  a y x2 M  b

Algoritmo simprl

Página 19 de 115
 La solución numérica elemental con Octave

3.3.4 Integración compuesta de Romberg

Emplea la regla compuesta de los trapecios para una primaria aproximación y

posteriormente utiliza la técnica de Richardson para mejoramientos de las
aproximaciones.

1 k 2

Rk ,1   Rk 1,1  hk 1  f (a  (2i  1)hk )  k  2,3,..., n
2 i 1 

Para cada k  2,3, 4,..., n y j  2,..., k se puede reducir la notación expresando

Rk , j 1  Rk 1, j 1
Rk , j  Rk , j 1 
4 j 1  1

Conformando un diagrama característico

R1,1
R2,1 R2,2
R3,1 R3,2 R3,3
.
. O
Rn ,1 L L L Rn ,n

O también es posible expresarlo cómo:

b
a
f ( x)dx R( J , J )
mediante la generación de una tabla de aproximaciones R (J, K) para J≥K y con
R(J+1,J+1) como la respuesta final. Las aproximaciones R(J, K) se almacenan en una
matriz triangular inferior especial. Los elementos de R(J,0) de la columna 0 se
calculan usando la regla trapezoidal secuencial basado en 2 J subintervalos de [a, b],
entonces R (J, K) se calcula utilizando la regla de Romberg.
Los elementos de la fila J son:

R( J , K  1)  R( J  1, K  1)
R( J , K )  R( J , K  1) 
4K 1
para 1  K  J .

Algoritmo de romber

Página 20 de 115
 La solución numérica elemental con Octave

3.3.5 Regla trapezoidal recursiva

Para aproximar:
J
h 2
f ( x)dx   f ( xk 1 )  f ( xk )
b

a 2 k 1

Mediante el uso de la regla del trapecio y sucesivamente cada vez mayor el número
de subintervalos de [a, b]. En la iteración J-ésima para f (x) en los 2J+1 puntos
igualmente espaciados.

Algoritmo trapezoidal

3.4 FÓRMULAS GAUSSIANAS

Presentación general

b n

  ( x) f ( x)dx   ai f ( xi )
a i 1

Con  ( x) una función de peso, que en el caso de  ( x)  1 se presenta la forma más

simple.
b
ai    ( x) Li ( x)dx Li ( x) = función multiplicadora de Lagrange.
a

Los x1 ,..., xn representan los ceros del polinomio Pn ( x) , de grado n, de una familia que
verifica ortogonalidad:

  ( x)P ( x)P ( x)dx  0

a
n m para mn

con los coeficientes dependientes de  ( x), entonces  ( x) ejerce influencia sobre

ai , xi aunque no esté taxativamente presente en la fórmula de Gauss.

3.4.1 Fórmula de Gauss-Legendre

 ( x)  1
 a, b    1,1
Los polinomios ortogonales son polinomios de Legendre

1 dn 2
n 
x  1
n
Pn ( x)  n P0 ( x)  1
2 n ! dx
Página 21 de 115
 La solución numérica elemental con Octave

2 1  xi2 
xi raíces del polinomio y ai  2
n 2 Pn 1 ( xi )
Tanto los xi como los ai están tabulados para reemplazar en

b n

 f ( x)dx   ai f ( xi )
a i 1

Algoritmo gausslegendre

3.4.2 Formula de Gauss-Laguerre

 n

e
x
f ( x)dx  ai f '( xi )
0 i 1

d n x n
Ln ( x)  e
dx n
e x 
x

Y
(n !) 2
ai 
xi  L '( xi )
2

3.4.3 Formula de Gauss-Chebyshev

  n
1
f ( x)
 dx     f ( xi )
1 1 x 2
 n  i 1

El polinomio n-simo de Chebyshev: Tn ( x)  cos(n arccos( x))

xi los ceros del polinomio Tn ( x) .

3.5 Cuadratura adaptativa utilizando la regla de Simpson’s

Para aproximar la integral:

M
f ( x)dx   ( f ( x4 k 4 )  4 f ( x4 k 3 )  2 f ( x4 k 2 )
b
a
k 1

 4 f ( x4 k 1 )  f ( x4 k )

Página 22 de 115
 La solución numérica elemental con Octave

La regla compuesta de Simpson’s es ampliada para los 4M subintervalos [x 4k-4, x4k],

donde [a,b]=[x0,x4M] y x4k-4+j=x4k-4+jhk, para cada k=1,…M y J=1,…4.

Algoritmo adapt

3.6 Regla de Simpson’s

El algoritmo “srule”, es una modificación de la regla de Simpson’s. El resultado es un

vector Z que contiene los resultados de la regla de Simpson en el intervalo [a 0, b0]. El
algoritmo “adapt” utiliza srule para llevar a cabo la regla de Simpson’s en cada uno de
los subintervalos generado por el proceso de cuadratura adaptativa.

Algoritmo srule

3.7 Cuadratura de Gauss-Legendre

Para aproximar la integral:

ba N
WN ,k f (t N ,k )
b

a
f ( x)dx 
2 k 1
Por muestreo f(x) en los N puntos desigualmente espaciados t  N
N , k k 1 .
Los cambios de variable:
ab ba ba
t  x y dt  dx
2 2 2

Son usados. Las abscisas x  N

N , k k 1 y el correspondiente peso W  N
N , k k 1 se deben
obtener de una tabla de valores conocidos.

Algoritmo gauss_q

Página 23 de 115
 La solución numérica elemental con Octave

ECUACIONES DIFERENCIALES ORDINARIAS

4.1 METODOS PARA RESOLVER P.V.I.

P.V.I. y '  f (t , y) at b y(a)  

En esta ecuación f es una función real dada de en t , donde y es una función

incógnita de variable independiente t . Tanto y como pueden ser vectores, en cuyo
caso tendremos un sistema de ecuaciones diferenciales de primer orden.

Una ecuación diferencial conjuntamente con una condición inicial constituye un

problema de valor inicial; siendo:
 y '  f (t , y )

 y ( x0 )  y0

4.1.1 Método de Euler

Se basa en la ecuación diferencial, en el tamaño de paso y en el paso particular de la

aproximación, con el error local de truncamiento = (h) .

wi 1  wi  hf (ti , wi ) i  0,1,2,..., N  1 , (ecuación de diferencias)

h=tamaño de paso, con distribución uniforme en [a,b].

ti  a  ih, para cada i  0,1,2,..., N .
wi  y(ti ) denota el valor exacto de la solución en t i
se suele emplear para el método una expresión del tipo:
0     0
i 1  i  hf (ti , wi )  i 1 i  0,1, 2,..., N 1

Siendo i error de redondeo asociado con i .

El mínimo valor de E (h) se da para h  2
M
O sea h superiores inducirán a incrementar el error total en la aproximación.

4.1.2 Serie de Taylor

El método de Euler es el método de Taylor de orden uno.

Para órden dos:
w0  0.5,
 h 
wi 1  wi  h (1  )( wi  ti2  1)  hti 
 2 
Página 24 de 115
 La solución numérica elemental con Octave

Para órden cuatro:

w0  0.5,
 h h 2 h3 h h2 h h 2 h3 
wi 1  wi  h (1    )( wi  ti )  (1   )hti  1    
2

 2 6 24 3 12 2 6 24 
i  0,1,..., N  1.
para  i en (ti , ti 1 )
yi 1  yi ( n ) hn
 i 1 (h)   T (ti , yi )  f ( n )  i , y(si )  i  0,1, 2,..., N 1
h  n  1!
con el error local de truncamiento = (h ) n  derivada admitida .
n

4.1.3 Métodos de Runge-Kutta

Fórmulas más difundidas

1  hf (t , y )
 h 1 
 2  hf  t  , y   O(h 2 )
 2 2
 h 2 
 3  hf  t  , y   O (h3 )
 2 2 
 4  hf  t  h, y   3  O (h 4 )

1
y (t  h)  y (t )  1  2 2  2 3   4  cuarto orden
6

Algoritmo rk4

4.1.4 MÉTODOS MULTIPASO

Dado el P.V.I. y '  f (t , y) at b y (a )  b

wi 1 en ti 1  m  Z , y1 será:

La ecuación de diferencias
wi 1  cm1wi  cm2 wi 1  ...  c0 wi 1m  h bm f (ti 1, wi 1 )  bm1 f (ti , wi )  ...  b0 f (ti 1m , wi 1m ) 
(b  a)
m  paso>1; ti = paso de red; h  ; a0 , a1 ,..., am1 y b0 , b1 ,..., bm son
N
constantes.

i  m  1, m,..., N 1

Página 25 de 115
 La solución numérica elemental con Octave

Valores iniciales dados w0   , w1  1 , w2   2 ,..., wm1   m1

[Link] Métodos de Adams-Bashforth: EXPLICITOS

A. De dos pasos
w0   w1  1
h
wi 1  wi  3 f (ti , wi )  f (ti 1 , wi 1 ) i  1,2,..., N  1
2
5 3
El error local de truncamiento es  i 1 (h)  y ( i ) h 2 con algún i entre
12
( ti 1 , ti 1 ).

B. de tres pasos
w0   w1  1 w2   2
h
wi 1  wi  [23 f (ti , wi )  16 f (ti 1 , wi 1 )  5 f (ti 2 , wi 2 )]
12
i  2,3,..., N  1
3
El error local de truncamiento es  i 1 (h)  y ( 4) ( i )h3 con algún i entre
8
(ti 2 , ti 1 ) .

C. De cinco pasos
w0   , w1  1 , w2  2 , w3  3 , w4  4
h
wi 1  wi  1901 f (ti , wi )  2774 f (ti1, wi1 )  2616 f (ti2 , wi2 ) 1274 f (ti3 , wi3 )  251 f (ti4 , wi4 )
720

i  4,5,..., N  1
El error local de truncamiento es  i 1 (h) 
95 ( 6)
y ( i ) h 5 con algún i entre
288
(ti 4 , ti 1 ) .

Para los implícitos se debe usar  ti1 , f (ti1 , y(ti1 ))  como un nodo de interpolación
ti 1

extra para aproximar la integral 

ti
f (t , y (t ))dt .

Página 26 de 115
 La solución numérica elemental con Octave

[Link] Métodos de Adams-Moulton: IMPLICITOS

A. De dos pasos

Las diferencias ordinarias de orden 1 y 2, están dadas por:

f n  f n 1  f n ,
2 f n  f n  2  2 f n 1  f n

para sustituir a las diferencias ordinarias en la expresión anterior y agrupando

términos semejantes, obtenemos:
Dados w0   y w1  1
h
wi 1  wi  5 f (ti1 , wi1 )  8 f (ti , wi )  f (ti1, wi1 )
12
i  1, 2,..., N  1

El error local de truncamiento es:

1 (4)
 i 1 (h)   y ( i )h3 con ti 1  i ti 1
24
B. De tres pasos
w0   w1  1 w2   2
h
wi 1  wi  [9 f (ti 1 , wi 1 )  19 f (ti , wi )  5 f (ti 1 , wi 1 )  f (ti 2 , wi 2 )]
24
i  1,2,..., N  1
19 (5)
El error local de truncamiento es:  i 1 (h)   y ( i )h4 con ti 2  i ti 1
720

C. De cuatro pasos
w0   w1  1 w2   2 w3  3

h
wi 1  wi   251 f (ti 1 , wi 1 )  646 f (ti , wi )  264 f (ti 1, wi 1 )  106 f (ti 2 , wi 2 ) 19 f (ti 3 , wi 3 ) 
720
i  3, 4,..., N  1

El error local de truncamiento es:

3 (6)
 i 1 (h)   y ( i )h5 con ti 3  i ti 1
160

Algoritmo abm

Página 27 de 115
 La solución numérica elemental con Octave

4.1.5 Método de Milne

Técnicas denominadas de predicción-corrección, empleando el explicito para predecir

la aproximación y la implícita para corregirla.
Otros métodos multipasos se generan por integración o interpolación de polinomios
en intervalos t j , ti 1  con j  i  1 parar aproximar a y(ti 1 ) .
Así está el método de Milne, entre ti 3 , ti 1  :
4h
wi 1  wi 3   2 f (ti , wi )  f (ti1 , wi1 )  2 f (ti2 , wi2 )
3

Algoritmo milne

4.1.6 Predictor-Corrector

Si usamos el primer orden (explícito) el método de Euler hacia adelante para

inicializar el método de Crank-Nicolson, obtenemos el método de Heun (también llama
método de Euler mejorado), que es un método de segundo orden explícito de Runge-
Kutta

un*1  un  hf n
h
un 1  un  [ f n  f (t n 1 , un*1 )]
2
El paso explícito constituye un factor de predicción, mientras que la implícita se llama
corrector. Este tipo de métodos disfrutan de la orden de precisión del método
corrector. Sin embargo, el ser explícito, se somete a una restricción de estabilidad que
suele ser el mismo que el del método predictor. Por lo tanto, no son adecuados para
integrar un problema de Cauchy en intervalos no acotados.

Algoritmo heun

4.1.7 USO DE LA EXTRAPOLACIÓN

condición inicial w0    y(a) y para w1 se utiliza la técnica de Euler.

wi 1  wi 1  2hf (ti , wi ) i 1
Al llegar al valor t, se corrige el punto extremo que comprende las dos aproximaciones
del punto medio, generando una aproximación w(t , h) a y (t ) dada por.

y (t )  w(t , h),  k h 2 k
k 1

Con k ctes., ligadas a las derivadas de y (t ) , pero es independiente de h.

Página 28 de 115
 La solución numérica elemental con Octave

SISTEMAS LINEALES

Ax  b
Planteo generalizado:
R1 a11 x1  a12 x2  ...  a1n xn  b1
R2 a21 x1  a22 x2  ...  a2 n xn  b2

Rn an1 x1  an 2 x2  ...  ann xn  bn

5.1 REPRESENTACIÓN DE UN SISTEMA LINEAL. SUSTITUCIÓN HACIA ATRÁS

a11 x1  a12 x2  ...  a1n xn  b1

a21 x1  a22 x2  ...  a2 n xn  b2
Un sistema puede llevarse a la forma:

an1 x1  an 2 x2  ...  ann xn  bn

 a11 a12 a1n a1,n 1 

0 a22 a2 n a2,n 1 
A ai ,n1  bi c  1, 2,..., n
  con para
i
 
0 0 ann an,n 1 

Permitiendo la sustitución hacia atrás.

Sustitución hacia atrás, solamente si todos los elementos de la diagonal son no ceros.
Primero calcula xN  bN / aNN y entonces realiza:

bk   j k 1 akj x j
N

xk  para k  N  1, N  2,...,1.
akk
Algoritmo backsub

Triangularización superior, seguido por la sustitución hacia atrás. Para construir la

solución AX = B, primero reducción de la matriz aumentada
[A | B] a la forma triangular superior y luego realizar la sustitución hacia atrás

Algoritmo uptrbk (utilizado por backsub)

Página 29 de 115
 La solución numérica elemental con Octave

5.2 TÉCNICA DE GAUSS-JORDAN

Partiendo de una matriz aumentada, por hipótesis tenemos

a11  0, pues det( A)  0
a i(,in)1
La solución se obtiene usando xi  para cada i  1,2,3,..., n.
a ii( i )
Prescindiendo de la sustitución hacia atrás en la eliminación gaussiana.
 n3 n  n3 n 
Requiere   n   productos/divisiones y    sumas/ restas.
2

 2 2  2 2

Algoritmo gjelim

5.3 PIVOTEO

Una pista para seleccionar es tomar el menor p  k tal que:

a (pk,k)  max ai(,kk) y efectuar  Rk    Rp 

k i  h
(Técnica parcial de pivoteo)

Otra técnica es la de reescalonamiento de columnas. El primer paso del

procedimiento consiste en definir, para cada renglón, un factor escalar si
Por medio de:

si  máx | aij | con 1  j  n cuando si =0 (solución múltiple); de no ser así, el

intercambio adecuado de renglones para poner ceros en la primer columna se
determina seleccionando el menor entero k con:

ak1 ai , j
 max
sk sj
j 1,2,..., n

Y se realiza  R1    Rk 
Algoritmo gauss con pivoteo

Página 30 de 115
 La solución numérica elemental con Octave

5.4 FACTORIZACION MATRICIAL

Si A es no singular, admite el producto SI (triangular superior por triangular inferior), la

factorización es de gran provecho, para la resolución de sistemas lineales.
Si la eliminación gaussiana para Ax  b se efectúa sin cambios de filas, A podrá
factorizarse como el producto IS.

A=IS

Con

 1 0 0
m 1 
 21 
I=  0  Matriz triangular inferior cij  0 para i  j
 
 mn1 mn ,n 1 1
 a11(1) a12(1) a1(1)n 
 (2) 
 0 a22 
S=  ( n 1)  Matriz triangular superior c  0 para i  j
an 1,n ij
 
 0 0 an( n,n) 

La solución para ISx=B puede ser obtenida definiendo y=Sx, entonces se resuelven
los dos sistemas:
1°) Iy=B para y
2°) Sx=y para x

5.4.1 Factorización LU

Tenemos que demostrar que la descomposición se puede hacer a una matriz de

orden n = k. Es entonces la matriz de orden k. Empezamos esta serie en sub-matrices
de la forma
 Ak 1 r 
A 
 st akk 

Donde r y s son vectores columna, ambos con k - 1 componentes.
Tenga en cuenta que la matriz Ak-1 es de orden k-1. Por lo tanto, por hipótesis de
inducción esta puede ser descompuesta en forma A k-1 = Lk-1Uk-1. Utilizando las
matrices Lk-1 y Uk-1 formamos las siguientes matrices:

 Lk 1 0  p 
L  ; U  U k 1 
 mt 1 0 ukk 
 
Página 31 de 115
 La solución numérica elemental con Octave

Donde m y p son vectores columna, ambos con k-1 componentes. Tenga en cuenta
que m, p y ukk son desconocidos.
Así, se establece que la matriz A se puede descomponer en LU trataremos de
determinarlos.
Realización del producto LU, se deduce que:

 Lk 1U k 1 Lk 1 p 
LU   
 mt S m p  ukk 
t
 k 1

Estudiemos ahora la ecuación LU=A, esto es:

 Lk 1U k 1 Lk 1 p   Ak 1 r 
  
 mtU m p  ukk   s t
t
akk 
 k 1

De esta igualdad concluimos

Lk 1U k 1  Ak 1 ,
Lk 1 p  r ,
m tU k 1  s t ,
m t p  ukk  akk .

Tenga en cuenta que la primera ecuación tiene por hipótesis la inducción, y por lo
tanto Lk 1 y U k 1 son determinadas de manera única. Además, ni Lk 1 ni U k 1
son singulares (o Ak 1 también seria singular, contradiciendo la hipótesis). Por lo
tanto:

Lk 1 p  r  p  Lk11r ,
mtU k 1  s t  mt  s tU k11 ,
mt p  ukk  akk  ukk  akk  mt p

Por lo tanto p, m y ukk son determinadas unívocamente en este orden, L y U son

determinados únicamente. Al final:

det( A)  det(L) * det(U )  1* det(U k 1 ) * ukk

 u11u22......uk 1,k 1ukk ,

Página 32 de 115
 La solución numérica elemental con Octave

Cabe señalar que la descomposición LU constituye uno de los algoritmos más

eficientes para el cálculo del determinante de una matriz.

La resolución del sistema Ax = b implica hacer Ly = b, Ux = y

Factorización con pivoteo

Algoritmo lufact

Según determinadas variaciones para la técnica se conocen los métodos:

 Doolittle: requiere que i11  i22  ...  inn  1

 Crout: requiere que los elementos de la diagonal principal de S sean todos 1.
 Choleski: requiere iii  sii para cada i.

Este tipo de factorización presentado es ventajoso cuando no se requiere el

intercambio de filas para ajustar el error de redondeo.

Matrices de permutación

Se utiliza cuando se requieren intercambios de renglones

 Si P es matriz de permutación tiene inversa y P1  Pt .
 Para cualquier A, no singular, habrá una P tal que PAx  Pb se resuelve sin
intercambiar filas, con lo que PA se puede factorizar como PA=IS
P1  Pt  A   P I  S
t
Y al ser

( P t I ) , no será triangular inferior (salvo que P=Identidad).

Ejemplo:
 1 0 0 3 3 1 
   
I   1 / 3 1 0 , S   0  2 14 / 3 
 2 / 3 0 1  0 0  5 / 3
   
P será
1 0 0
 
P  0 0 1
0 1 0
 

Página 33 de 115
 La solución numérica elemental con Octave

Matrices características

a) Choleski : para una definida positiva requiere menos operaciones que IDI t
pero implica hallar n raíces cuadradas; con I con aii  1 y D diagonal con
aii  0 . A puede factorizarse como IDI t .
b) Doolittle o de Crout : para matrices de banda que disponen sus elementos no
nulos en derredor de la diagonal. Las más comunes, p=q=2 y los de p=q=4,
pentadiagonales.
Las primeras, asociadas a aproximaciones de splines cúbicos o
aproximaciones lineales por parte en problemas de valores de contorno, los
pentadiagonales de utilidad para problemas de valores de frontera.
El trabajar con matrices de banda facilita los algoritmos de factorización (por el
número de operaciones).

5.5 TÉCNICAS DE REPETICIÓN PARA SISTEMAS LINEALES

 
Transformar el sistema original Ax  b en uno equivalente x  Wx  c para cierta
matriz fija W y un vector c.
La sucesión se va calculando como x ( k )  Wx ( k 1)  c k  1, 2,3...,

Métodos iterativos de resolución de sistemas lineales

Estas técnicas consisten en transformar un sistema de la forma
Au = b
en una ecuación de punto fijo de la forma: u = Mu + c
de tal manera que, al hacer iteraciones de la forma:
n n−1
u = Mu +c
n
se obtenga que u converge hacia u, la solución del sistema original.
Consideremos el sistema de ecuaciones
2x – y = 1
-x + 2y – z = 0
−y + 2z = 1

Buscar la solución de este sistema es equivalente a buscar un vector u = (x, y,

z) que verifique que
1 y xz 1 y
x ,y ,z 
2 2 2

Hacer iteraciones de esta ecuación de punto fijo consiste en partir de una

aproximación inicial (x1 , y1 , z1 ) y hacer iteraciones de la forma:
1  yn1 x z 1  yn1
xn  , yn  n 1 n1 , zn 
2 2 2
Página 34 de 115
 La solución numérica elemental con Octave

En este caso, la solución exacta del sistema es u = (1, 1, 1).

Si hacemos iteraciones del esquema anterior a partir de la aproximación inicial
1
u = (0, 0, 0), obtenemos que
1 0 1 00 1 0 1
x2   , y2   0 , z2  
2 2 2 2 2
las sucesivas iteraciones se van aproximando a la solución u = (1, 1, 1).
n n−1
Si el esquema iterativo u = M u + c converge hacia un vector u, entonces u
verifica que u = Mu + c

Existen diferentes métodos para convertir un sistema de la forma Au = b en

una ecuación de punto fijo u = Mu + c. Todas se basan en descomponer A
de la forma A = L + D + U, donde D es la matriz diagonal que corresponde a la
parte diagonal de A, L es la matriz triangular inferior que corresponde a la parte
de A situada por de- bajo de la diagonal, y U es la matriz triangular superior
que corresponde a la parte de A situada por encima de la diagonal.

5.5.1 Técnica de Jacobi

Este método consiste en tomar

−1
Mj = D (−L − U )
−1
cj = D b
Es el que se ha utilizado en el ejemplo anterior. El paso de una iteración a
otra del método de Jacobi puede expresarse de la siguiente forma:

 a1 2u 2n 1  ...  a1 N u N
n 1
 b1
u n
1 
a1 1
 a2 1u1n 1  a2 3u3n 1 ...  a2 N u N
n 1
 b2
u n
2 
a2 2
 a N 1u1n 1  a N 2u 2n 1...  a NN 1u N
n 1
1  bN
u n
N 
a NN
Algoritmo jacobi

5.5.2 Técnica de Gauss-Seidel

Este método consiste en tomar:

M GS  ( D  L) 1 (U )
cGS  ( D  L) 1 b

Página 35 de 115
 La solución numérica elemental con Octave

A efectos prácticos, la aplicación de este método no requiere el cálculo directo

de la matriz inversa (D + L)−1 ,puesto que el paso de una iteración a otra
puede hacerse de la siguiente forma:

 a12u2n 1  ...  a1N u Nn 1  b1

u 
n
1
a11
 a21u1n  a23u3n 1...  a2 N u Nn 1  b2
u 
n
2
a22
 a N 1u1n  a N 2u2n ...  a NN 1u Nn 1  bN
u 
n
N
a NN

Si hacemos un barrido para el cálculo de la solución de arriba hacia abajo, y

vamos actualizando las componentes del vector aproximación según las
vamos calculando, obtenemos el método de Gauss-Seidel. Por tanto,
básicamente, podemos decir que la diferencia entre el método de Gauss-
Seidel y el método de Jacobi es que en el método de Gauss-Seidel se actualiza
el vector aproximación después del cálculo de cada componente, y en el caso
de Jacobi se actualiza sólo al final, después de haber calculado todas las
componentes por separado.

Algoritmo gseid

5.6 Técnica de relajación

El objetivo de este método es intentar mejorar el método de Gauss-Seidel

introduciendo un parámetro de relajación w. Se toman, en este caso:
M = (D + wL)−1 ((1 − w)D − wU )
w

cw = w (D + wL)−1 b

Estas nuevas matrices permiten realizar un promediado entre el resultado

obtenido por Gauss-Seidel y el estado de la solución en la etapa anterior,
de la forma siguiente:

Página 36 de 115
 La solución numérica elemental con Octave

 a12u 2n1  ...  a1N u Nn1  b1

u w
n
1  (1  w)u1n1
a11
 a21u1n ...  a2 N u Nn1  b2
u w
n
2  (1  w)u 2n1
a22
 a N 1u1n ...  a NN 1u Nn 1  bN
u Nn  w  (1  w)u 2n1
a NN

La elección del parámetro w es, en general, un problema difícil. Sin embargo, en

el caso de matrices tridiagonales, es decir, matrices con todos los elementos
nulos salvo la diagonal principal y sus codiagonales, el siguiente resultado
muestra la forma de calcular el valor óptimo de w.
Si A es una matriz tridiagonal y ρ(Mj) < 1,entonces el valor de w que optimiza la
velocidad de convergencia del método es:
2
wopt 
1  1   (M j ) 2
Como puede observarse de la expresión anterior, el valor de wopt se
encuentra siempre entre 1 y 2.

CONCLUSIÓN: el teorema se concluye teniendo en cuenta que cualquier norma

de una matriz es siempre mayor o igual que su radio espectral.
Este resultado se puede generalizar un poco al caso de matrices irreducibles
de la siguiente forma:
A) Una matriz A es irreducible si un sistema de la forma Au = b no puede
descomponerse en dos sub-sistemas independientes de dimensión
menor.
B) Una matriz es irreducible si el cambio de cualquier valor del vector b del
sistema Au = b afecta a todos los elementos del vector u.
Si A es una matriz irreducible y se verifica que

aii   j i aij i

a jj  i  j aij j
ó

Con la desigualdad estricta en al menos una fila o columna, entonces los métodos
iterativos convergen.

OTRA FORMA DE PLANTEAR LA TÉCNICA DE RELAJACIÓN

( k 1) rii( k )
xi
(k )
x i q
aii
Página 37 de 115
 La solución numérica elemental con Octave

vector residuo ri (k1) asociado al vector, definido por r  b  Ax entonces habrá un

ri ( k )   r1i( k ) , r2(ik ) ,..., rni( k ) 

t
que se busca converja a cero.

 Si 0 q1 , técnicas de subrrelajación, para sistemas no convergentes por G-S.

 q1 , para acelerar convergencia de sistemas convergentes por G-S (métodos
SOR o sobrerelajación).
Para el cálculo se usará:
 i 1
( k 1) 
   ij j 

n
q
 i  ij j
( k 1)
xi
(k )
 (1  q) x
i  b  a x (k )
 a x
aii  j 1 j i 1 
Matricialmente
x ( k )  ( D  qI )1 (1  q) D  qS  x ( k 1)  q( D  qI ) 1 b

5.7 GRADIENTE CONJUGADO

Se busca una secuencia de n direcciones conjugadas, y luego se calculan los

coeficientes.

La solución x* de Ax  b es también el único minimizador de:
1 T
f ( x)  x Ax  bT x , x R n
2
rk=b-Axk

rk es el negativo del gradiente de f en x = xk, con lo que se puede mostrar, ya que las
direcciones pk son conjugadas entre sí, que una expresión para tales direcciones es:

pkT Ark
pk 1  rk  T pk
pk Apk

Algoritmo conjgrad

Página 38 de 115
 La solución numérica elemental con Octave

APROXIMACIÓN

6.1 POLINOMIOS DE CHEBYSCHEV

p( x)  a0 H 0 ( x)  ...  am H m ( x) polinomio de aproximación

6.2 APROXIMACIÓN MÍNIMO –MÁXIMO (O DE CHEBYSHEV)

P( x)  mx  c
y (b)  y(a) y (a)  y ( x2 ) (a  x2 )  y (b)  y(b)
m c 
ba 2 2(b  a)

6.3 APROXIMACIÓN POR VALORES PROPIOS

El teorema de descenso Gerschgorin es útil para determinar los límites globales a los
valores propios de una matriz Anxn. El término "global", los límites que encierran todos
los valores propios.

 
 n

Dada Anxn Ri   z  C / z  aii   aij  Ri=es el círculo en el plano complejo con
 j 1 
 j i  
aii como centro.
n

Los auto-valores de A están dentro de R= Ri (círculo de Gerschgorin)

i 1
El función Gerschgorin devuelve la parte inferior y superior de los límites globales a
los valores propios de una matriz tridiagonal simétrica A.

Algoritmo Gerschgorin

6.4 Técnicas de Aproximación

6.4.1 Métodos de Potencias

Sea una matriz A que posee una base de autovectores tal que en módulo su
autovalor máximo λmax es único. Sea un vector u1 no ortogonal al subespacio
engendrado por los autovectores del autovalor λ max, entonces, si definimos la
secuencia

Limn sign((u n , u n1 )) u n  max

Página 39 de 115
 La solución numérica elemental con Octave

un
Limn ( sign((u , u n n 1
))) n
es un autovector de max
un
u n 1
un  A
kun 1k

Además, dicho autovector tiene norma 1.

sign((u n , u n1 )) es el signo del producto escalar de u n u n1 , es decir
sign((u n , u n1 ))  1 si (u n , u n1 )  0 y sign((u n , u n1 ))  1 si (u n , u n1 )  0
Para calcular el valor propio dominante 1 y su vector propio asociado V1 para la
matriz Anxn . se supone que los valores propios n tienen la propiedad de posición

dominante 1  2  3  ...  n  0 .

Algoritmo power1

6.4.2 Método de la potencia inversa

El método anterior también se puede utilizar para el cálculo del autovalor de

módulo menor λmin , teniendo en cuenta que:
1
min 
max{  i autovalores de A -1 }
'

Por tanto, si aplicamos el método anterior a A−1 , obtenemos que la

u n 1 1
secuencia uA 1
verifica que Limn sign((u n , u n1 )) u n 
u n 1 min
un
Limn sign((u n , u n1 )) es un autovector de min
un
En los casos prácticos, se evita calcular directamente A−1 , y se obtiene un
resolviendo el sistema
u n 1
Au  n 1
n

Para calcular el valor propio dominante 1 y su vector propio asociado Vj para la

matriz Anxn. Se supone que los n valores propios tienen la propiedad 1  2  ...  n y
que α es un número real tal que:

Página 40 de 115
 La solución numérica elemental con Octave

 j    i   para cada i  1,2,..., j  1, j  1,..., n

Algoritmo invpow

6.4.3 Método simétrico de Potencias

x (0)t Ax (0)
:q 
(0)
x
x (0)t x (0)
Para elegir q se considera que si x es vector propio de A respecto al valor propio ,
x t Ax
se dará que Ax   x , entonces x Ax   x x y   t
t t
xx
6.4.4 METODO CLASICO DE JACOBI

El clásico método de Jacobi, o simplemente el método de Jacobi, es un método

numérico utilizado para determinar la auto-valores y la auto-vectores de matrices
simétricas. Teniendo en cuenta la matriz A, se realiza una secuencia de rotación
A1  A ; A2  U1t A1U1  A3  U 2t A2U 2 
 ...  Ak 1  U kt AkU k  D
donde U i , i  1,2,..., k son matrices de rotación, y D es una matriz diagonal.
La secuencia de las matrices de Ak es calculada a través de una repetición de:
Ak 1  U kt AkU k ( K  1,2,...)
Como A1  A , obtenemos :
Ak 1  U kt U kt 1...U 2tU1t AU 1U 2 ...U k 1U k  V t AV

Donde V  U1U 2 ...U k 1U k .

Con la hipótesis de que Ak ≈ D, obtenemos que D = VtAV, donde V es una matriz

ortogonal, porque la matriz V es producto de matrices ortogonales. Así, D contiene la
auto-valores de A y V contiene su correspondientes
auto-vectores (en columnas), es decir, la j-ésima columna de V es un auto-vector
correspondiente a auto-valor  j .

Para calcular tanto los auto-valores como los auto-vectores  ,V 

j
n
j j 1 de la matriz

simétrica Anxn :

Algoritmo jacobi1

Página 41 de 115
 La solución numérica elemental con Octave

6.4.5 Técnicas de Deflación

Se emplean para la determinación de aproximaciones a los otros valores propios de

una matriz después de haber encontrado la aproximación al valor propio dominante.
Consiste en generar una nueva matriz B, con valores propios iguales a los de A, salvo
que el valor propio dominante de A se sustituya por el valor propio 0 en B.
La técnica es muy sensible a errores de redondeo.

6.4.6 Técnica de Householder

Para reducir la matriz simétrica Anxn a la forma tridiagonal mediante el uso de n - 2

transformaciones de householder.
La transformación puede ser expresada como:
QA' Q  A'wu T  uwT
que nos da el siguiente procedimiento de cálculo que ha de llevarse a cabo con:
i  1,2,..., n  2 .

Algoritmo house

6.4.7 Algoritmo QR

Supongamos que A es una matriz simétrica real, En la sección anterior vimos cómo el
método de Householder se utiliza para construir una matriz tridiagonal similar. El
método QR se utiliza para encontrar todos los valores propios de una matriz
tridiagonal. rotaciones similares a los que utiliza el método de Jacobi se utilizan para
construir una matriz ortogonal Q1 = Q y una triangular superior R = R1 de manera que
A1 = A tiene la factorización:

A(i 1) = R(i )Q(i ) = Q(i )t A(i )  Q(i )  Q(i )t A(i )Q(i )
Algoritmo qr1

Algoritmo qr2

6.5 Ortogonalización de Gram-Schdmidt

Este método permite obtener de una matriz Anxn una matriz Q con columnas
ortonormales y R, una matriz triangular superior de tal forma que A = Q*R.

Algoritmo grams

Página 42 de 115
 La solución numérica elemental con Octave

SISTEMAS DE ECUACIONES NO LINEALES

Forma general:
f1 ( x1 , x2 ,..., xn )  o
f 2 ( x1 , x2 ,..., xn )  o

f n ( x1 , x2 ,..., xn )  o

7.1 MÉTODO DE NEWTON

x ( k )  F ( x ( k 1) )  x ( k 1)  J ( x ( k 1) )( 1) G( x ( k 1) ) k 1

Método de Newton para sistemas no lineales, esperando convergencia cuadrática, si

se conoce x
(0)
preciso y exista J (q )1 .

 g1 ( x ) g1 ( x ) g1 ( x ) 

 x x2 xn 
 1

J ( x)   
 
 g n ( x ) g n ( x ) g n ( x ) 
 x1 x2 xn 

En la búsqueda de un punto fijo, vimos la transformación x  g (x) . Para el caso de

G de D  R en R tiene un punto fijo en
n n
funciones de R n en R n => Una función
p  D si G( p)  p

Algoritmo newtons

7.2 TÉCNICAS CUASI-NEWTON

Método de Broyden

G ( x (1) )  G ( x (1) )  J ( x (1) )  x (1)  x (0)    x (1)  x (0) 

t

A1  J ( x )  
(0)
z

x (1)  x (0)
e

 (1) ( 0 )
x  x es el vector ortogonal

Algoritmo broyden
Página 43 de 115
 La solución numérica elemental con Octave

7.3 TÉCNICAS DE MAYOR PENDIENTE

x (1)  x (0)  ag ( x (0) ), a0

La dirección de la máxima disminución del valor de g en x es la dirección dada por

 g (x)

7.4 PROBLEMA DE VALOR DE FRONTERA PARA ECUACIONES IFERENCIALES

ORDINARIAS Y EN DERIVADAS PARCIALES

y ''  f ( x, y, y ') en  a, b
Sujeta a: y(a)    y(b)  

Cuando y ''  f ( x, y, y ') en a, b tiene la forma:

y ''  p( x) y ' q( x) y  r ( x) en a, b el problema es lineal,
Si además Satisface:
1. p( x), q( x), y r ( x) continuas en  a, b 
2. q( x) 0 en  a, b 
El problema tiene solución única.

7.5 TECNICA DE DISPARO PARA EL PROBLEMA LINEAL

y1 ( x)  y ''  p( x) y '  q( x) y  r ( x), a  x  b, y(a)   , y ' (a)  0

y2 ( x )  y ''  p( x) y '  q( x) y, a  x  b, y(a)  0, y ' (a)  1

Ambos problemas tienen solución única y se cumple:

  y1 (b)
y ( x)  y1 ( x)  y2 ( x) (Solución única de la ecuación diferencial
y2 (b)
con valor en la frontera).

Algoritmo disparo Lineal

Página 44 de 115
 La solución numérica elemental con Octave

7.6 TECNICA DE DISPARO PARA EL PROBLEMA NO LINEAL

y’’=f(x,y,y’), a≤x≤b, y(a)=ɑ, y’(a)=t

donde se eligen los parámetros

lim
t  tk de manera que y(b, tk )  y(b)  
k 
la solución en x como en t, requerirá la forma de los PVI como:
 y’’(x,t)=f(x,y(x,t),y’(x,t)), a≤x≤b, y(a,t)=ɑ, y’(a,t)=t y

 z’’(x,t)=(∂f/∂y)(x,y,y’)z(x,t)+(∂f/∂y’)(x,y,y’)z’(x,t) , a≤x≤b, z(a,t)=0,

z’(a,t)=1
simbolizando z(x,t) a (∂y/∂t)(x,t)
entonces en cada iteración se deberá resolver ambos PVI, con:
y (b, tk 1 )  
tk  tk 1 
z (b, tk 1 )

Algoritmo disparo no Lineal

7.7 DIFERENCIAS FINITAS PARA PROBLEMAS LINEALES

para la ecuación general de orden dos y ( xi )  p( xi ) y ( xi )  q( xi ) y( xi )  r ( xi )

'' '

en [a,b] con alpha=y(a), beta=y(b) y N el número de etapas.

Al desarrollar y en el tercer polinomio de Taylor evaluando xi entre xi 1 y xi 1 ,

suponiendo que y  C 4 [ xi 1 , xi 1 ] tendremos:

h 2 '' h3 ''' h 4 ( 4) 
y( xi 1 )  y( xi  h)  y( xi )  hy ( xi ) 
'
y ( xi )  y ( xi )  y (i )
2 6 24

h 2 '' h3 ''' h 4 ( 4) 
y( xi 1 )  y( xi  h)  y( xi )  hy ( xi ) '
y ( xi )  y ( xi )  y (i )
2 6 24
Sumando estas dos ecuaciones y aplicando el teorema del valor intermedio, podemos
''
llegar a las fórmulas de diferencias centradas para y ( xi ) e y ' ( xi ) (aquí no
se detalla).
La utilización de las fórmulas de diferencias centradas de:
y '' ( xi )  p( xi ) y ' ( xi )  q( xi ) y( xi )  r ( xi ) genera la ecuación:

Página 45 de 115
 La solución numérica elemental con Octave

y ( xi 1 )  2 y ( xi )  y ( xi 1 )  y ( xi 1 )  y ( xi 1 ) 
 p( xi )    ...

2
h 2h

q( xi ) y ( xi )  r ( xi ) 
h2
12

2 p( xi ) y ''' (i )  y ( 4) ( i ) 
Este método de deferencias finitas con (h 2 ) para el error de

truncamiento, junto con las condiciones de frontera alpha=y(a), beta=y(b) se emplea

para definir: w0   , wN 1   y
 wi 1  2wi  wi 1  wi 1  wi 1 
 p ( xi )    ...
h2 2h
q( xi ) wi  r ( xi ) para i  1,2,..., N

El sistemas de ecuaciones resultantes se expresa en forma de una matriz tridiagonal

de N*N como Aw  b .

Página 46 de 115
 La solución numérica elemental con Octave

ECUACIONES EN DERIVADAS PARCIALES

8.1 ECUACIÓN GENERAL

Af xx 
 2 Bf xy  Cf yy  Df x  Ef y  Ff  G

con A, B,..., G funciones reales de x e y , en base a su analogía con el

discriminante de secciones cónicas:
Ax 2  Bxy  Cy 2  Dx  Ey  F  0 ,
se las describirá como hiperbólicas, elípticas y parabólicas.

El objetivo es buscar un método numérico de aproximación de la solución exacta.

Se formularán las EDP clásicas en base a su aproximación por diferencias finitas

8.2 ELÍPTICA O DE POISSON

2 f 2 f
( x, y )  2 ( x, y )  g ( x, y ) Si g ( x, y)  0 se obtiene la ecuación de Laplace.
x 2 y

Si la función representa un campo de temperaturas, las restricciones del problema a

través de la distribución de temperatura en el borde de la región, son las condiciones
de frontera de Dirichlet, dadas por:

f ( x, y)  r ( x, y) ( x, y) en S, frontera de la región U.

U  ( x, y) / a xb, c y d 
El método de resolución consiste en una adaptación del método de diferencias finitas
para problemas con valor de frontera tomando:

( a  b) (d  c)
h y k obteniendo una cuadrícula con los puntos de red, en
n m
cada punto de red utilizamos la serie de Taylor en la variable x alrededor de xi para
generar la fórmula de la diferencias centradas.

2 f f ( xi 1 , y j )  2 f ( xi , y j )  f ( xi 1 , y j ) h2  4 f
( xi , yi )   ( i , y j )
x 2 h2 12 x 4
 i  ( xi 1 , xi 1 )
Lo mismo para la variable y alrededor de yi
2 f f ( xi , y j 1 )  2 f ( xi , y j )  f ( xi 1 , y j 1 ) h2  4 f
( xi , yi )   ( xi , j )
y 2 k2 12 y 4
 j  ( y j 1 , y j 1 )
Página 47 de 115
 La solución numérica elemental con Octave

El uso de estas fórmulas permite expresar la ecuación de Poisson en los puntos

( xi , yi ) como:

f ( xi 1 , y j )  2 f ( xi , y j )  f ( xi 1 , y j ) f ( xi , y j 1 )  2 f ( xi , y j )  f ( xi 1, y j 1 )
 
h2 k2
h2  4 f h2  4 f
 g ( xi , y j )  ( i , y j )  ( xi , j )
12 x 4 12 y 4
i  1, 2,...(n  1) j  1, 2,...( m  1)
Para condiciones de frontera:
f ( x0 , y j )  r ( x0 , y j ) j  0,1,..., m
f ( xn , y j )  r ( xn , y j ) j  0,1,..., m
f ( xi , y0 )  r ( xi , y0 ) i  1, 2,..., n  1
f ( xi , ym )  r ( xi , ym ) i  1, 2,..., n  1
Para un error de truncado del orden O(h2  k 2 ) , se tendrá:

con i , j aproximando a f ( xi , y j )
 h  2  h
2

2    1 i , j  ( i 1, j  i 1, j )    ( i , j 1  i , j 1 )  h2 g ( xi , y j ) (A)

 k   k
(diferencias centradas)
i  1, 2,..., n  1 j  0,1,..., m
0, j  r ( x0 , y j ) j  0,1,..., m
n, j  r ( xn , y j ) j  0,1,..., m
0, j  r ( xi , y0 ) i  1, 2,..., n  1
0, j  r ( xi , ym ) i  1, 2,..., n  1
A través de (A) se aproxima af ( x, y) en los puntos
( xi 1 , y j ),( xi , y j ),( xi 1 , y j ),( xi , y j 1 ) y ( xi , y j 1 ) .

Algoritmo dirich (Método de Dirichlet para la ecuación de Laplace’s)

Para aproximar la solución de u xx ( x, y )  u yy ( x, y )  0 más
R  {( x, y ) : 0  x  a, 0  y  b) con u ( x,0)  f1 ( x), u ( x, b)  f 2 ( x),
para 0  x  a, y u (0, y )  f 3 ( y ), u (a, y )  f 4 ( y ), para 0  y  b.
Asumiendo que x  y  h y que los enteros n y m existen
tal que a  nh y b  mh.

Página 48 de 115
 La solución numérica elemental con Octave

8.3 ECUACION PARABÓLICA

Para la ecuación de Calor (parabólica)- unidimensional

u  2u
2 2 su fórmula en diferencias regresivas (incondicionalmente estable)
t x
será:

ui , j  ui , j 1 ui 1, j  2ui , j  ui 1, j

 2 0
k h2
para i  1,2,..., m  1; y j  1,2,...

este método de diferencias tiene la representación matricial

(1  2 )   0 0   u1, j   u1, j 1 

   u   u 
   2, j   2, j 1 
 0 0  .    . 
    
    .   . 
 0 0   (1  2 ) um 1  um 1, j 1 

O sea Aw(j) = w(j-1) para toda j = 1, 2, ….

Indicando la necesidad de resolver un sistema lineal para hallar w(j) a partir de w(j−1). ,
siendo estable para cualquier λ .
Dado que λ > 0, la matriz A es definida positiva y con dominancia diagonal estricta,
tridiagonal, pudiendo resolverse por algoritmos ya vistos (caso Crout).
Para evitar la falta de precisión, generada por el error de truncado, se requiere que los
intervalos de tiempo sean mucho más pequeños que los de espacio (h > k), haciendo
necesario un método que permita tomar valores similares para h y k, y estable para
todo λ.

Algoritmo crnich (Crank-Nicholson método para la ecuación del calor)

Para aproximar la solución de ut ( x, t )  c 2u xx ( x, t ) más
R  {( x, t ) : 0  x  a, 0  t  b) con u ( x,0)  f ( x), para 0  x  a.
y u (0, t )  c1 , u (a, t )  c2 , para 0  t  b.

8.4 ECUACIÓN HIPERBÓLICA

Se presenta la ecuación hiperbólica a través de una situación problema

 2u 2  u
2

Sea la EDP c
t 2 x 2
con las condiciones de frontera u(0, t) = 0 = u(1, t), 0 < t < 1
y con las condiciones iniciales u(x, 0) = sen(πx), 0 ≤ x ≤ 1 y ut(x, 0) = 0,
0≤x≤1

Página 49 de 115
 La solución numérica elemental con Octave

Sn exacta: u(x, t) = sen(πx)cos(2πt)

La división rectangular del dominio (x,t): : 0≤ x ≤ a, 0 ≤ t ≤ b, denominando

u(x, t) = ui,j
u(x + h, t) = ui+1,j
u(x − h, t)) = ui−1,j
u(x, t + k)) = ui,j+1
u(x, t − k)) = ui,j−1
la ecuación en diferencias será
ui , j 1  2ui , j  ui , j 1 ui 1, j 1  2ui , j 1  ui 1, j 1
 c2
k2 h2
Haciendo ɑ=ck/h
ui , j 1  2ui , j  ui , j 1 =ɑ2( ui 1, j 1  2ui , j 1  ui 1, j 1 )
Colocando en términos de ui,j+1 y reordenando, se tiene
ui , j 1  (2  2 2 )ui , j   2 (ui 1, j  ui 1, j )  ui , j 1
La ecuación (arriba) es aplicable para i = 2, 3, ..., n − 1 y j = 2, 3, ...,m − 1, con el paso
(j + 1 de tiempo requiriendo de los j-ésimo y los (j − 1) ésimo pasos; los valores del (j −
1)ésimo paso vienen dados por las condiciones iniciales ui,1 = f(xi) pero, los valores del
j ésimo paso no
se suelen proporcionar, por lo que se usa g(x) para conseguir los valores de este
paso.
Para asegurar la estabilidad de este método es necesario que:
ɑ = ck/h ≤ 1.

Algoritmo finedif (diferencia finita para la ecuación de onda)

Para aproximar la solución de utt ( x, t )  c 2u xx ( x, t ) más
R  {( x, t ) : 0  x  a, 0  t  b) con u (0, t )  0, u (a, t )  0,
para 0  t  b, y u ( x,0)  f ( x), ut ( x,0)  g ( x), para 0  x  a.

Página 50 de 115
La solución numérica elemental con Octave

ANEXO
GNU OCTAVE

Página 51 de 115
 La solución numérica elemental con Octave

Octave es un lenguaje de ordenador de alto nivel, orientado al cálculo numérico. Y

su vez es un programa que permite interpretar este lenguaje de forma interactiva,
mediante órdenes o comandos. Estas órdenes pueden ser introducidas a través de un
entorno en forma de terminal o consola de texto. Por supuesto, Octave incluye,
además, la posibilidad de ser utilizado de forma no interactiva, interpretando las
órdenes contenidas en un fichero.
Estas características se pueden resumir diciendo que Octave es un lenguaje
interpretado orientado al cálculo numérico matricial. Por supuesto, existe una gran
variedad de lenguajes interpretados orientados a las matemáticas, algunos de ellos
con licencia libre, como Máxima, Axiom, Scilab, EULER, Sage o R y otros con licencia
privativa, como Matlab, Maple, Mathematica o S.
Los lenguajes interpretados, entre los cuales destaca la menor velocidad de
proceso que lenguajes compilados como C, C++, Fortran, etc. Sin embargo, Octave y
el resto de los lenguajes interpretados presentan una serie de ventajas muy
significativa: la razón entre el esfuerzo empleado en el desarrollo del programa y el
resultado final es muy pequeña, es decir, permiten obtener resultados razonables con
un esfuerzo no demasiado importante. Este hecho, que está relacionado con la
eliminación de la segunda etapa dentro del ciclo clásico de desarrollo de programas
(escritura-compilado-depurado), es vital en:
a) Desde sus orígenes, Octave es software libre, gracias a lo cual, su comunidad de
desarrolladores y usuarios ha ido creciendo hasta llegar a ser muy significativa.
Es un proyecto en colaboración cuyos miembros nos prestarán ayuda cuando la
solicitemos, a través de sus listas de correo. Éstas son públicas y una gran
fuente de información.
b) Su lenguaje de programación es altamente compatible con el de Matlab, el
conocido entornode cálculo numérico (con licencia privativa) desarrollado por
The MathWorks, Inc. De hecho, la mayor parte de la materia desarrollada en las
siguientes páginas es, de hecho, un análisis del amplio lenguaje que es común a
Matlab y Octave y al que, indistintamente, denominaremos “lenguaje Octave” o
“lenguaje Matlab”.
c) Está disponible en numerosos plataformas, entre ellas sistemas Unix (para los
que fue originalmente desarrollado), Windows y MacOsX. Su código fuente
contiene una potente librería de cálculo matricial para C++. De hecho, ha sido
diseñado de forma que sea extensible dinámicamente, a través de nuevos
procedimientos escritos en C++. Estas extensiones son programas “especiales”
que tienen la extensión .oct y, como programas compilados, cuentan con un
rendimiento mucho mayor que las funciones escritas en el lenguaje interpretado
de Octave.
d) interfaces gráficas de usuario que, en la actualidad, ofrece prestaciones más
interesantes para un usuario medio: qtOctave. Este programa recibe su nombre
de las bibliotecas QT. desarrolladas con la finalidad de facilitar la creación de
interfaces gráficas.

Página 52 de 115
 La solución numérica elemental con Octave

Interfaz qtOctave

La interfaz gráfica permite, entre otras cosas, trabajar con el interprete de órdenes de
Octave, ver el listado de variables, editar un archivo .m, y ver el histórico de órdenes.

qtOctave está disponible:

A través de su página web, [Link] Esta página es el

centro de comunicación con el autor y la comunidad de usuarios de qtOctave.

Página 53 de 115
 La solución numérica elemental con Octave

function[c,yc,err,P]=bisect(a,b,delta) METODO DE LA BISECCIÓN

# Ejemplo de carga
# [c,yc,err,P]=bisect(0,2,0.001)
# función 'x*sin(x)-1'
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');
err=0;
P=[a b err];
ya=feval(f,a);
yb=feval(f,b);
if ya*yb>0,break,end
max1=1+round((log(b-a)-log(delta))/log(2));
for k=1:max1,
c=(a+b)/2;
yc=feval(f,c);
if yc==0,a=c;
b=c;
elseif yb*yc>0,
b=c;
yb=yc;
else
a=c;
ya=yc;
end
err=abs(b-a)/2;
P=[P;a b err];
if b-a<delta,break,end
end
c=(a+b)/2;
yc=feval(f,c);
err=abs(b-a)/2;
disp(c);
disp(yc);
disp(err);
disp(' X0 X1 ERROR');
disp(P);
end

Dada la función f(x)= xsenx-1, hallar una raíz en el intervalo [0,2] por el
método de la bisección y una tolerancia de 10-3. Igual para x3+4x2-10 en
-4
[1,2] y tolerancia de 10

octave:4> [c,yc,err,P]=bisect(0,2,0.001)
ingresef(x) entre comillas en términos de x 'x*sin(x)-1'
'x*sin(x)-1'
1.1138
-5.3832e-004
4.8828e-004
X0 X1 ERROR
0.00000 2.00000 0.00000
1.00000 2.00000 0.50000
1.00000 1.50000 0.25000
1.00000 1.25000 0.12500
1.00000 1.12500 0.06250
1.06250 1.12500 0.03125
1.09375 1.12500 0.01563
1.10938 1.12500 0.00781
1.10938 1.11719 0.00391

Página 54 de 115
 La solución numérica elemental con Octave
1.11328 1.11719 0.00195
1.11328 1.11523 0.00098
1.11328 1.11426 0.00049
c = 1.1138
yc = -5.3832e-004
err = 4.8828e-004
P =
0.00000 2.00000 0.00000
1.00000 2.00000 0.50000
1.00000 1.50000 0.25000
1.00000 1.25000 0.12500
1.00000 1.12500 0.06250
1.06250 1.12500 0.03125
1.09375 1.12500 0.01563
1.10938 1.12500 0.00781
1.10938 1.11719 0.00391
1.11328 1.11719 0.00195
1.11328 1.11523 0.00098
1.11328 1.11426 0.00049

function fixed_point(p0,N) METODO DEL PUNTO FIJO

%Fixed_Point(p0, N) aproxima la raíz de la ecuación f(x) = 0
%reescrita en la forma x = g(x), partiendo con una aproximación inicial p0.
%La técnica iterativa se implementa N veces.
%El usuario tiene que entrar g(x)abajo

g=input('ingrese f(x) entre comillas en términos de x');

g=inline(g,'x');
i = 1;
p(1) = p0;
tol = 1e-05;
while i <= N
p(i+1) = g(p(i));
if abs(p(i+1)-p(i)) < tol %stopping criterion
disp('el procedimiento fue exitoso luego de k iteraciones')
k = i
disp('la raíz de la ecuación es')
p(i+1)
return
end
i = i+1;
end

if abs(p(i)-p(i-1)) > tol | i > N

disp('el procedimiento falló')

Página 55 de 115
 La solución numérica elemental con Octave
disp('Condición |p(i+1)-p(i)| < tol no se satisfizo')
tol
disp('Por favor, examine la secuencia de iteraciones')
p = p'
disp('Si observa convergencia, incremente el número máximo de
iteraciones')
disp('Si hay divergencia, intente otro p0 o reescriba g(x)')
disp('si |g''(x)|< 1 cerca de la raíz')
end

Dada la función f(x)=x-x3+4x2-10/3x2+8x, encontrar la raíz por el método de

punto fijo
Para la función y = x - (x.^3 + 4*x.^2 - 10)/(3*x.^2 + 8*x);

octave:7> fixed_point(1.5, 10)

ingresef(x) entre comillas en términos de x 'x- x^3-4*x^2+10'
'x- x^3-4*x^2+10'
el procedimiento falló
Condición |p(i+1)-p(i)| < tol no se satisfizo
tol = 1.0000e-005
Por favor, examine la secuencia de iteraciones
p =
1.5000e+000
-8.7500e-001
6.7324e+000
-4.6972e+002
1.0275e+008
-1.0849e+024
1.2771e+072
-2.0827e+216
NaN
NaN
NaN
Si observa convergencia, incremente el número máximo de iteraciones
Si hay divergencia, intente otro p0 o reescriba g(x)
si |g'(x)|< 1 cerca de la raíz

Tomando su equivalente
y = sqrt(10./x - 4*x);

octave:8> fixed_point(1.5, 10)

ingresef(x) entre comillas en términos de x 'sqrt(10./x - 4*x)'
'sqrt(10./x - 4*x)'
el procedimiento falló
Condición |p(i+1)-p(i)| < tol no se satisfizo
tol = 1.0000e-005
Por favor, examine lasecuencia de iteraciones
p =
1.50000 - 0.00000i
0.81650 - 0.00000i
2.99691 - 0.00000i
0.00000 - 2.94124i
2.75362 + 2.75362i
1.81499 - 3.53453i
2.38427 + 3.43439i
2.18277 - 3.59688i
2.29700 + 3.57410i
2.25651 - 3.60656i
2.27918 + 3.60194i

Página 56 de 115
 La solución numérica elemental con Octave

function [p0,y0,err,P] = newton(p0,delta,max1) METODO DE NEWTON

# Ejemplo carga [p0,y0,err,P] = newton(0.25*pi,0.01,12)
# función 'cos(x)-x'
# derivada '-sin(x)-1'
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');
epsilon = 1.0842e-19;
df = input('ingrese la derivada de f(x) entre comillas, en términos de x');
df = inline(df,'x');
P(1) = p0;
P(2) = 0;
P(3) = 0;
P(4) = 0;
P(5) = 0;
y0 = feval(f,p0);
for k=1:max1,
df0 = feval(df,p0);
if df0 == 0,
dp = 0;
else
dp = y0/df0;
end
p1 = p0 - dp;
y1 = feval(f,p1);
err = abs(dp);
relerr = err/(abs(p1)+eps);
p0 = p1;
y0 = y1;
P = [P;p1 y0 df0 err relerr];
if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
end
disp(' Xi F(Xi) F(Xi) ErrorAbs ErrorRel ');
disp(P);
disp('Error Absoluto: '),disp(err);
disp(' Error Relativo: '), disp(relerr);
end

Dada la función f(x)=x-x3+4x2-10, hallar una raíz por el método de Newton

con aprox. inicial 1.5, tolerancia para la raíz de 10-4 y 10-5 para los
valores de f, con 10 iteraciones.
Sea la función x^3+4*x^2-10
Su derivada es 3*x^2+8*x

octave:9> newton(1.5,0.0001,10)
ingresef(x) entre comillas en términos de x 'x^3+4*x^2-10'
'x^3+4*x^2-10'
ingrese la derivada de f(x) entre comillas, en términos de x '3*x^2+8*x'
'3*x^2+8*x'
Xi F(Xi) F(Xi) ErrorAbs ErrorRel
1.50000 0.00000 0.00000 0.00000 0.00000
1.37333 0.13435 18.75000 0.12667 0.09223
1.36526 0.00053 16.64480 0.00807 0.00591
1.36523 0.00000 16.51392 0.00003 0.00002
ErrorAbsoluto:
3.2001e-005
Error Relativo:
2.3440e-005
ans = 1.3652

Página 57 de 115
 La solución numérica elemental con Octave

function [p1,y1,err,P] = secant(p0,p1,delta,max1) METODO DE LA SECANTE

# Ejemplo de carga [p1,y1,err,P] = secant(0.5,pi/4,0.01,40)
# función 'cos(x)-x'
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');
epsilon = 1.0842e-19;
P(1) = p0;
P(2) = p1;
P(3) = 0;
P(4) = 0;
P(5) = 0;
P(6) = 0;
y0 = feval(f,p0);
y1 = feval(f,p1);
for k=1:max1,
df = (y1-y0)/(p1-p0);
if df == 0,
dp = 0;
else
dp = y1/df;
end
p2 = p1 - dp;
y2 = feval(f,p2);
err = abs(dp);
relerr = err/(abs(p2)+eps);
p0 = p1;
y0 = y1;
p1 = p2;
y1 = y2;
P = [P;p0 p2 p1 y1 err relerr];
if (err<delta)|(relerr<delta)|(abs(y2)<epsilon), break, end
end
disp('P0 P2 Xi F(Xi) ErrAbs ErrorRelativo');
disp(P);
disp('Error Absoluto:'),disp(err),
disp('Error Relativo:'),disp(relerr);
end

Para la función cosx-x , hallar una raíz con aprox. iniciales 0.5 y π/4,
tolerancia para p1=0.01, para f de 0 01 y 10 iteraciones, por el método de la
secante
octave:10> [p1,y1,err,P]=secant(0.5,pi/4,0.01,10)
ingresef(x) entre comillas en términos de x 'cos(x)-x'
'cos(x)-x'
P0 P2 Xi F(Xi) ErrAbs ErrorRelativo
0.50000 0.78540 0.00000 0.00000 0.00000 0.00000
0.78540 0.73638 0.73638 0.00452 0.04901 0.06656
0.73638 0.73906 0.73906 0.00005 0.00267 0.00362
ErrorAbsoluto:
0.0026740
ErrorRelativo:
0.0036181
p1 = 0.73906
y1 = 4.5177e-005
err = 0.0026740
P =
0.50000 0.78540 0.00000 0.00000 0.00000 0.00000
0.78540 0.73638 0.73638 0.00452 0.04901 0.06656
0.73638 0.73906 0.73906 0.00005 0.00267 0.00362

Página 58 de 115
 La solución numérica elemental con Octave

function [c,yc,err,P] = regula(a,b,delta,max1) METODO DE LA POSICIÓN FALSA

# Ejemplo carga [c,yc,err,P] = regula(1.38,1.0,0.001,10)

# función 'x^3+4*x^2-10'
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');
epsilon = 1.0842e-19;
P = [a b 0 0 0];
ya = feval(f,a);
yb = feval(f,b);
if ya*yb > 0, break, end
for k=1:max1,
dx = yb*(b - a)/(yb - ya);
c = b - dx;
ac = c - a;
yc = feval(f,c);
end
if yc == 0,
break;
elseif yb*yc > 0,
b = c;
yb = yc;
else
a = c;
ya = yc;
end
dx = min(abs(dx),ac);
err = abs(dx);
if abs(dx) < delta, break, end
if abs(yc) < epsilon, break, end
end
err = abs(dx);
P = [a b c yc err];
disp(' X0 X1 Raiz FunctValor Error')
disp(P)
end

Utilizar la falsa posición para a) x3+4x2-10 en [ 1.38 1], tolerancias de 10-3

y 10 iteraciones

octave:12> [c,yc,err,P] = regula(1.38,1.0,0.001,10)

ingresef(x) entre comillas en términos de x 'x^3+4*x^2-10'
'x^3+4*x^2-10'
X0X1 Raiz FunctValor Error
1.380000 1.362203 1.362203 -0.049906 0.017797
c = 1.3622
yc = -0.049906
err = 0.017797
P =
1.380000 1.362203 1.362203 -0.049906 0.017797

Página 59 de 115
 La solución numérica elemental con Octave

function [p,Q]=steff(p0,delta,epsilon,max1) METODO DE STEFFENSEN

%Ingreso - f es la function ingresada como string 'f'

% - df es la derivada f ingresada como string 'df'
% - p0 es la aproximación inicial a un cero de f
% - delta es la tolerancia para p0
% - epsilon es la tolerancia para los valores de la función y
% - max1 es el número máximo de iteraciones
%Salida - p es la aproximación de Steffensen al cero
% - Q matriz conteniendo la secuencia de Steffensen
%Inicializar la matriz R
f = input('ingrese f(x) entre comillas en términos de x');
f = inline(f,'x');
df = input('ingrese la derivada de f(x) entre comillas, en términos de x');
df = inline(df,'x');
R = zeros(max1,3);
R(1,1)=p0;
for k=1:max1
for j=2:3

%Cálculo del Denominador en el método de Newton-Raphson

nrdenom=feval(df,R(k,j-1));

%Condicional calcula las aproximaciones de Newton-Raphson

if nrdenom==0
'división cero en el método de Newton-Raphson '
break
else
R(k,j)=R(k,j-1)-feval(f,R(k,j-1))/nrdenom;
end

%Cálculo de Denominador en el proceso de aceleración de Aitkens

aadenom=R(k,3)-2*R(k,2)+R(k,1);

% condicional en las aproximaciones de Aitkens

if aadenom==0
'división por cero en aceleración de Aitkens '
break
else
R(k+1,1)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;
end
end

%Corte y salida del programa si ocurre división por cero

if (nrdenom==0)|(aadenom==0)
break
end

%criterio de parada; se determinan p y matriz Q

err=abs(R(k,1)-R(k+1,1));
relerr=err/(abs(R(k+1,1))+delta);
y=feval(f,R(k+1,1));
if (err<delta)|(relerr<delta)|(y<epsilon)
p=R(k+1,1);
Q=R(1:k+1,:);
break
end
end

Página 60 de 115
 La solución numérica elemental con Octave
Aplicar el método de Steffensen. Recordando el problema de f(x)=cos x –x con
newton.

octave:14> [p,Q]=steff(0.25*pi,0.01,0.01,4)
ingrese f(x) entre comillas en términos de x 'cos(x)-x'
'cos(x)-x'
ingrese la derivada de f(x) entre comillas, en términos de x '(x)-sin(x)-1'
'(x)-sin(x)-1'
p = 0.73819
Q =
0.78540 0.70046 0.76834
0.73819 0.00000 0.00000

function [y,b] = horner(a,z) Método de Horner

%HORNER Horner algorithm

% Y=HORNER(A,Z) calcula
% Y = A(1)*Z^N + A(2)*Z^(N-1) + ... + A(N)*Z + A(N+1)
% usando el algoritmo de division sintética de Horner.
n = length(a)-1;
b = zeros(n+1 ,1);
b(1) = a(1);
for j=2:n+1
b(j) = a(j)+b(j -1)*z;
end
y = b(n+1);
b = b(1:end -1);
return

function [roots ,iter]= newtonhorner(a,x0 ,tol ,nmax) Método de Horner

%NEWTONHORNER Newton -Horner

% [roots ,ITER ]= NEWTONHORNER(A,X0) calcula las raíces del
% polinomio
% P(X) = A(1)*X^N + A(2)*X^(N-1) + ... + A(N)*X +
% A(N+1)
% usando el método de Newton -Horner partiendo de los datos iniciales X0. El
método para dada raíz
% luego de 100 iteraciones o luego que el valor absoluto de la diferencia
entre dos iteraciones consecutivas es menor que
% 1.e -04.
% [roots ,ITER ]= NEWTONHORNER(A,X0 ,TOL ,NMAX) permite definir la
tolerancia en el criterio de parada y el número máximo de iteraciones.
if nargin == 2
tol = 0.0001; nmax = 100;
elseif nargin == 3
nmax = 100;
end
n=length(a)-1; raiz = zeros(n ,1); iter = zeros(n ,1);

Página 61 de 115
 La solución numérica elemental con Octave
for k = 1:n
% iteraciones de
niter = 0; x = x0; diff = tol + 1;
while niter <= nmax & diff >= tol
[pz ,b] = horner(a,x); [dpz ,b] = horner(b,x);
xnew = x - pz/dpz; diff = abs(xnew -x);
niter = niter + 1; x = xnew;
end
if niter >= nmax
fprintf(' falla de convregencia con máximo ' ,...
'número de iteraciones\n');
end
% Deflación
[pz ,a] = horner(a,x); raiz(k) = x; iter(k) = niter;
end
return

octave:36>[roots ,iter]= newtonhorner([-3 1 0],2,0.001,50)

roots =
3.3334e-001
-2.8529e-006
iter =
6
2

function [p2,y2,err,P]=muller(p0,p1,p2,delta,max1) Método de Muller

# Ejemplo carga [p2,y2,err,P]=muller(0.5,pi/3.5,pi/4,0.01,12)

# función 'cos(x)-x'
f=input('ingrese f(x):');
f=inline(f,'x');
epsilon = 1.0842e-19;
P(1) = p0;
P(2) = p1;
P(3) = p2;
P(4) = 0;
P(5) = 0;
y0 = feval(f,p0);
y1 = feval(f,p1);
y2 = feval(f,p2);
for k=1:max1;
h0 = p0 - p2;
h1 = p1 - p2;
c = y2;
e0 = y0 - c;
e1 = y1 - c;
det1 = h0*h1*(h0-h1);
a = (e0*h1 - h0*e1)/det1;
b = (h0^2*e1 - h1^2*e0)/det1;
if b^2 > 4*a*c,
disc = sqrt(b^2 - 4*a*c);
else
Página 62 de 115
 La solución numérica elemental con Octave
disc = 0;
end
if b < 0
disc = - disc;
end
z = - 2*c/(b + disc);
p3 = p2 + z;
if abs(p3-p1) < abs(p3-p0)
u = p1;
p1 = p0;
p0 = u;
v = y1;
end
end

Emplear el método de Muller para cosx–x con po=0.5 p1=π/3.5, p2= π/3.5
octave:17> [p2,y2,err,P]=muller(0.5,pi/3.5,pi/4,0.01,0.001,12)
ingresef(x): 'cos(x)-x '
'cos(x)-x '
p2 = 0.78540
y2 = -0.078291
err = [](0x0)
P =
0.50000 0.89760 0.78540 0.00000 0.00000

function [p,y,err]=muller2(p0,p1,p2,delta,epsilon,max1) Método de Muller

# Ejemplo carga [p2,y2,err,P]=muller2(0.5,pi/3.5,pi/4,0.01,0.001,12)

# funcion 'cos(x)-x'
%Ingresos - f es la function como string 'f'
% - p0, p1, p2 son las aproximaciones iniciales
% - delta tolerancia para p0, p1, p2
% - epsilon la tolerancia para los valores de la función y
% - max1 número máximo de iteraciones
%salida- p aproximación de Muller al cero de f
% - y valor de la function y = f(p)
% - err error en la aproximación de p.
%Inicializar las matrices P e Y
f=input('ingrese f(x):');
f=inline(f,'x');
P=[p0 p1 p2];
Y=feval(f,P);
%Calcula a y b de formula (15)
for k=1:max1
h0=P(1)-P(3);h1=P(2)-P(3);e0=Y(1)-Y(3);e1=Y(2)-Y(3);c=Y(3);
denom=h1*h0^2-h0*h1^2;
a=(e0*h1-e1*h0)/denom;
b=(e1*h0^2-e0*h1^2)/denom;

%Suprime las raíces complejas

if b^2-4*a*c > 0
disc=sqrt(b^2-4*a*c);
else
disc=0;
Página 63 de 115
 La solución numérica elemental con Octave
end
%halla la menor raíz de (17)
if b < 0
disc=-disc;
end

z=-2*c/(b+disc);
p=P(3)+z;
%clasifica las entradas de p para hallar las dos más cercanas a p
if abs(p-P(2))<abs(p-P(1))
Q=[P(2) P(1) P(3)];
P=Q;
Y=feval(f,P);
end
if abs(p-P(3))<abs(p-P(2))
R=[P(1) P(3) P(2)];
P=R;
Y=feval(f,P);
end
%Reemplaza las entradas de P más alejadas de p con p
P(3)=p;
Y(3) = feval(f,P(3));
y=Y(3);

%Determina criterio de parada

err=abs(z);
relerr=err/(abs(p)+delta);
if (err<delta)|(relerr<delta)|(abs(y)<epsilon)
break
end
end

[p2,y2,err,P]=muller2(0.5,pi/3.5,pi/4,0.01,0.001,12)
ingresef(x): 'cos(x)-x'
'cos(x)-x'
p2 = 0.73897
y2 = 1.9532e-004
err = 0.046430

function [C,L]=lagrange(X,Y) POLINOMIO DE LAGRANGE

# Ejemplo carga [C,L]=lagrange([1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860

0.4554022 0.2818186 0.1103623])
%Ingreso - X es el vector de abcisas
% - Y es el vector de ordenadas
%salida - C es una matriz que contiene los coeficientes del
% polinomio interpolatorio de Lagrange
% - L es una matriz que contiene los coeficientes polinomiales de
% Lagrange

Página 64 de 115
 La solución numérica elemental con Octave

w=length(X);
n=w-1;
L=zeros(w,w);

%Forma los coeficientes polinomiales de Lagrange

for k=1:n+1
V=1;
for j=1:n+1
if k~=j
V=conv(V,poly(X(j)))/(X(k)-X(j));
end
end
L(k,:)=V;
end
%Determina los coeficientes del polinomio interpolatorio de Lagrange
C=Y*L;

Dados los siguientes valores x-y:

x=[1.0 1.3 1.6 1.9 2.2];
y=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
encontrar el polinomio interpolante de Lagrange

octave:19> [C,L]=lagrange([1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860

0.4554022
0.2818186 0.1103623])
C =
0.0018251 0.0552928 -0.3430466 0.0733913 0.9777351
L =
5.1440 -36.0082 93.3642 -106.2243 44.7243
-20.5761 137.8601 -338.2716 358.6008 -137.6132
30.8642 -197.5309 460.1852 -461.2346 167.7160
-20.5761 125.5144 -279.0123 268.2305 -94.1564
5.1440 -29.8354 63.7346 -59.3724 20.3292

function hp=hermite(x,y,dy) POLINOMIO DE HERMITE

% busca los polinomios interpolantes de Hermite para múltiples subintervalos

%Ingreso : [x,y],dy – puntos y sus derivadas
%salida: H = coeficientes de polinomio interpolante cúbico de Hermite
%
%Polinomio interpolante de Hermite. Considere el problema de hallar
%el polinomio de interpolaciónn para los N + 1 = 4 datos
%{(0, 0), (1, 1), (2, 4), (3, 5)}
%>>x = [0 1 2 3]; y = [0 1 4 5]; dy = [2 0 0 2]; xi = [0:0.01:3];
%>>H = hermite(x,y,dy);
%yi = ppval(mkpp(x,H), xi);

Página 65 de 115
 La solución numérica elemental con Octave

for n = 1:length(x)-1
H(n,:) = hermit(0,y(n),dy(n),x(n + 1)-x(n),y(n + 1),dy(n + 1));
end

function H = hermit(x0,y0,dy0,x1,y1,dy1)
% para obtener los coef. de Hermite para un conjunto de múltiples
subintervalos.
A = [x0^3 x0^2 x0 1; x1^3 x1^2 x1 1;
3*x0^2 2*x0 1 0; 3*x1^2 2*x1 1 0];
b = [y0 y1 dy0 dy1]'; %
H = (A\b)';

octave:20> H = hermite([0 1 2 3],[0 1 4 5],[2 0 0 2])

H =
0.00000 -1.00000 2.00000 0.00000
-6.00000 9.00000 0.00000 1.00000
0.00000 1.00000 0.00000 4.00000
octave:21> yi = ppval(mkpp([0 1 2 3],H), [0:0.1:3])
yi =
Columns 1 through 7:
0.00000 0.19000 0.36000 0.51000 0.64000 0.75000 0.84000
Columns 8 through 14:
0.91000 0.96000 0.99000 1.00000 1.08400 1.31200 1.64800
Columns 15 through 21:
2.05600 2.50000 2.94400 3.35200 3.68800 3.91600 4.00000
Columns 22 through 28:
4.01000 4.04000 4.09000 4.16000 4.25000 4.36000 4.49000
Columns 29 through 31:
4.64000 4.81000 5.00000

function p = neville1(x, f, xdach) MÉTODO DE NEVILLE

% Algoritmo Aitken- Neville-

%
% usar: p = neville(x, f, xdach)
%
% ingreso: vector x - ,x-coordenadas
% vector f - ,y-coordenadas
% escalar xdach – punto a evaluar
%
% salida: escalar p – valor interpolado en xdach

n = length(x);

Página 66 de 115
 La solución numérica elemental con Octave
for k = n-1:-1:1
f(1:k) = f(2:k+1) + ...
( xdach - x(n-k+1:n) ) ./ ...
( x(n-k+1:n) - x(1:k) ) .* ...
( f(2:k+1) - f(1:k) );
end
p = f(1);

Dados los siguientes valores x-y, generar el valor de x=1.4 por Neville
x=[1.0 1.3 1.6 1.9 2.2];
y=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];

octave:23> p=neville1([1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860 0.4554022

0.2818186 0.1103623],1.4)
p = 0.56685

function y = neville2(xx,n,x,Q) MÉTODO DE NEVILLE

% Algoritmo de Neville como una funcióon (guardar como "nev.m")

%
% ingresos:
% n = orden de interpolación (n+1 = # de puntos)
% x(1),...,x(n+1) x coordenadas
% Q(1),...,Q(n+1) y coordenadas
% xx=punto de evaluación para polinomio interpolante p
%
% salida: p(xx)
% x=[Link];
% Q=[132165 151326 179323 203302 226542 249633];
% neville2(2000,5,x,Q)

%rpta = 251654
for i = n:-1:1
for j = 1:i
Q(j) = (xx-x(j))*Q(j+1) - (xx-x(j+n+1-i))*Q(j);
Q(j) = Q(j)/(x(j+n+1-i)-x(j));
end
end

y = Q(1);

Recordando el ejercicio de neville1.

Con el algoritmo de neville2, para el mismo punto, grado 4

octave:23> p=neville1([1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860 0.4554022

0.28
18186 0.1103623],1.4)
p = 0.56685
octave:24> y=neville2(1.4,4,[1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860
0.4554022 0.2818186 0.1103623])
y = 0.56685

Página 67 de 115
 La solución numérica elemental con Octave

function nf = divdiff1(xi,fi) DIFERENCIAS DIVIDIDAS

%DIVDIFF calcula los the coeficientes para la diferencia hacia
% adelante de Newton,del polinomio interpolante asociado
% con un dado conjunto de puntos interpolantes
%
% secuencias de llamado:
% nf = divdiff1 ( xi, fi )
% divdiff1 ( xi, fi )
%
% ingresos:
% xi vector conteniendo los puntos interpolantes
% fi vector conteniendo los valores de la función
% la entrada I en este vector es el valor de la función
% asociado con la entrada i del vector'xi'
%
% salida:
% nf vector conteniendo coeficientes de la diferencia hacia
% adelante de Newton,del polinomio interpolante asociado
% con un dado conjunto de puntos interpolantes
% NOTA:
% para evaluar la forma de Newton , applicar rutina NF_EVAL
%ejemplo:> fi=[0.8;1.6;2.25;2.8;3.35];
% xi=[0.5;0.8;1.0;1.4;1.75];
%divdiff1 ( xi, fi );ans = 0.8000 2.6667 1.1667 -4.7685
%6.6669;

n = length ( xi );
m = length ( fi );

if ( n ~= m )
disp ( 'error divdiff : número de ordendas y número de valores de la
función deben ser iguales' )
return
end

nf = fi;
for j = 2:n
for i = n:-1:j
nf(i) = ( nf(i) - nf(i-1) ) / ( xi(i) - xi(i-j+1) );
end
end

Dados los siguientes valores x-y.

x=[1.0 1.3 1.6 1.9 2.2];
y=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
generar las diferencias divididas y el polinomio de Newton

octave:29> xi=[1.0 1.3 1.6 1.9 2.2];

octave:30> fi=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
octave:31> divdiff1 ( xi, fi )
ans =
0.7651977 -0.4837057 -0.1087339 0.0658784 0.0018251

Página 68 de 115
 La solución numérica elemental con Octave

function [C,D] = newpoly(X,Y) DIFERENCIAS DIVIDIDAS

%Ingreso - X vector de abcisas

% - Y vector de ordenadas
%salida - C es un vector que contiene los coeficientes
% del polinomio interpolante de Newton
% - D es la tabla de diferencias divididas

n=length(X);
D=zeros(n,n);
D(:,1)=Y';

for j=2:n
for k=j:n
D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));
end
end

%Determina los coeficientes del polinomio interpolante de Newton

C=D(n,n);

for k=(n-1):-1:1
C=conv(C,poly(X(k)));
m=length(C);
C(m)=C(m)+D(k,k);
end

Dados los siguientes valores x-y.

x=[1.0 1.3 1.6 1.9 2.2];
y=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
generar los coeficientes del polinomio interpolante de Newton
y la tabla de diferencias divididas

octave:29> xi=[1.0 1.3 1.6 1.9 2.2];

octave:30> fi=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
octave:32> [C,D] = newpoly(xi,fi)
C =
0.0018251 0.0552928 -0.3430466 0.0733913 0.9777351
D =
0.76520 0.00000 0.00000 0.00000 0.00000
0.62009 -0.48371 0.00000 0.00000 0.00000
0.45540 -0.54895 -0.10873 0.00000 0.00000
0.28182 -0.57861 -0.04944 0.06588 0.00000
0.11036 -0.57152 0.01182 0.06807 0.00183

function ls = linear_spline ( xi, fi ) SPLINE

%LINEAR_SPLINE calcula el interpolante lineal a trozos asociados a

Página 69 de 115
 La solución numérica elemental con Octave
% un conjunto dado de puntos de interpolación y los valores
% de la función
%
% forma de la llamada:
% ls = linear_spline ( xi, fi )
% linear_spline ( xi, fi )
%
% Entrada:
% xi vector que contiene los puntos de interpolación
% (orden ascendente)
% fi vector que contiene los valores de la función
%

%
% Salida:
% ls matriz de tres columnas
% - primer columna: puntos de interpolación
% - segunda columna: valores de la función
% - tercer columna: piezas lineales
%
% Nota:
% para evaluar el interpolante lineal a trozos aplicar el
% algoritmo SPLINE_EVAL para la salida de esta rutina
%

n = length ( xi );
m = length ( fi );

if ( n ~= m )
disp ('Número de ordenadas y valores de la función, debe ser igual')
return
end

b = zeros (n,1);
b(1:n-1) = ( fi(2:n) - fi(1:n-1) ) ./ ( xi(2:n) - xi(1:n-1) );

if ( nargout == 0 )
disp ( [ xi' fi' b ] )
else
ls = [ xi' fi' b ];
end

Dados los valores de x-f(x), encontrar el trazador lineal:

xi=[0.1 0.2 0.3 0.4]; fi=[-0.62049958 -0.2839668 0.00660095 0.24842440];
octave:7> xi=[0.1 0.2 0.3 0.4];
octave:8> fi=[-0.62049958 -0.2839668 0.00660095 0.24842440];
octave:9> linear_spline ( xi, fi )
0.10000 -0.62050 3.36533
0.20000 -0.28397 2.90568
0.30000 0.00660 2.41823
0.40000 0.24842 0.00000

La tercera columna da las pendientes de los trozos lineales

Los valores del pol se obtienen con “spline_eval”. Ver a continuación
function y = spline_eval ( sp, x ) ADAPTADOR CÚBICO

%SPLINE_EVAL evalúa una dada función spline function en un dado

Página 70 de 115
 La solución numérica elemental con Octave
% conjunto de valores de la variable independiente
%
% secuencia de llamado:
% y = spline_eval ( sp, x )
% spline_eval ( sp, x )
%
% ingresos:
% sp matriz conteniendo la información que define la
% función spline
% - típicamente será la la salida de
% LINEAR_SPLINE, CUBIC_NAK o CUBIC_CLAMPED
% x valor(es) de la variable independiente donde se
% evalúa la función spline
% - escalar o vector
%
% salida:
% y valor(es) de la función spline
%
% NOTA:
% rutina que acompaña a las CUBIC_NAK,
% CUBIC_CLAMPED y LINEAR SPLINE
%

[p c] = size ( sp );
npts = length ( x );

for i = 1 : npts
piece = max ( find ( sp(1:p-1) < x(i) ) );
if ( length ( piece ) == 0 )
piece = 1;
end

temp(i) = sp(piece,2) + sp(piece,3) * ( x(i) - sp(piece,1) );

if ( c == 5 )
temp(i) = temp(i) + sp(piece,4) * ( x(i) - sp(piece,1) )^2;
temp(i) = temp(i) + sp(piece,5) * ( x(i) - sp(piece,1) )^3;
end
end
if ( nargout == 0 )
disp ( temp )
else
y = temp;
end

Los valores del polinomio obtenidos con spline_eval

octave:10> sp=[3.3653 2.9057 2.4182 0];

octave:11> x=[0.1 0.2 0.3 0.4];
octave:12> spline_eval( sp, x )
-4.9904 -4.7486 -4.5068 -4.2650
function [A,df]=diffnew(X,Y) DIFERENCIACIÓN NUMERICA

%Ingreso - X vector abcisa 1xn

% - Y vector ordenada 1xn
%salida - A vector 1xn que contiene los coeficientes del polinomio de grado
% n de Newton
% - df derivada aproximada
A=Y;
N=length(X);
for j=2:N
Página 71 de 115
 La solución numérica elemental con Octave
for k=N:-1:j
A(k)=(A(k)-A(k-1))/(X(k)-X(k-j+1));
end
end

x0=X(1);
df=A(2);
prod=1;
n1=length(A)-1;

for k=2:n1
prod=prod*(x0-X(k));
df=df+prod*A(k+1);
end

function [L,n]=difflim(x,toler) DIFERENCIACIÓN NUMERICA

%Ingreso - x punto de diferenciación

% - toler tolerancia deseada
%Salida - L=[H' D' E']: H vector de tamaños de etapas
% D vector de derivadas aproximadas

% - E vector de cotas de error

% - n coordenda de la mejor aproximación

f=input('ingrese f(x) entre comillas en términos de x');

f=inline(f,'x');
max1=15;
h=1;
H(1)=h;
D(1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(1)=0;
R(1)=0;

for n=1:2
h=h/10;
H(n+1)=h;
D(n+1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(n+1)=abs(D(n+1)-D(n));
R(n+1)=2*E(n+1)*(abs(D(n+1))+abs(D(n))+eps);
end

n=2;

while((E(n)>E(n+1))&(R(n)>toler))&n<max1
h=h/10;
H(n+2)=h;
D(n+2)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(n+2)=abs(D(n+2)-D(n+1));
Página 72 de 115
 La solución numérica elemental con Octave
R(n+2)=2*E(n+2)*(abs(D(n+2))+abs(D(n+1))+eps);
n=n+1;
end

n=length(D)-1;
L=[H' D' E'];

Dada la función x2-3x, hallar su derivada en x=3 con tolerancia 0.01

difflim(3,0.01)
ingrese f(x) entre comillas en términos de x 'x^2-3*x'
'x^2-3*x'
ans =
1.00000 3.00000 0.00000
0.10000 3.00000 0.00000
0.01000 3.00000 0.00000

La salida es L=[H' D' E']: H vector de tamaños de paso;D el de [Link],E error

function [D,err,relerr,n]=diffext(f,x,delta,toler) TÉCNICA DE RICHARDSON

%Ingreso - f es la función
% - x es el punto de diferenciación
% - delta es la tolerancia para el error
% - toler es la tolerancia para el error relativo
%Salida - D es la matriz de derivadas aproximadas
% - err es la cota de error
% - relerr es la cota de error relativo
% - n es la coordenada de la "mejor aproximación"
% Si f es archivo M usar el llamado
%[D,err,relerr,n]=diffext(@f,x,delta,toler).
% Si f es una función anónima llamar como
%[D,err,relerr,n]=diffext(f,x,delta,toler).
%
err=1;
relerr=1;
h=1;
j=1;
D(1,1)=(f(x+h)-f(x-h))/(2*h);

while relerr > toler & err > delta & j <12
h=h/2;
D(j+1,1)=(f(x+h)-f(x-h))/(2*h);
for k=1:j
D(j+1,k+1)=D(j+1,k)+(D(j+1,k)-D(j,k))/((4^k)-1);
end
err=abs(D(j+1,j+1)-D(j,j));
relerr=2*err/(abs(D(j+1,j+1))+abs(D(j,j))+eps);
j=j+1;
end

Página 73 de 115
 La solución numérica elemental con Octave
[n,n]=size(D);

Dada la función x2-3x, hallar su derivada en x=3 con tolerancia 0.0001

octave:23> f=@(x)x^2-3*x;
octave:24> [D,err,relerr,n]=diffext(f,3,0.01,0.0001)
D =
3 0
3 3
err = 0
relerr = 0
n = 2

function s=traprl(a,b,M) Regla Compuesta del Trapecio

%Ingreso - f es el integrando ingresado como string 'f'

% - a y b son los límites de integración
% - M es el número de subintervalos
%Salida - s es la suma de la regla trapezoidal

f=input('ingrese f(x) entre comillas en términos de x');

f=inline(f,'x');
h=(b-a)/M;
s=0;

for k=1:(M-1)
x=a+h*k;
s=s+feval(f,x);
end
s=h*(feval(f,a)+feval(f,b))/2+h*s;

Hallar la integral de la función f(x)=sen x entre 0 y π/3

Trapezoidal compuesta
octave:29> traprl(0,pi/3,6)
ingrese f(x) entre comillas en términos de x 'sin(x)'
'sin(x)'
ans = 0.49873

function s=simprl(a,b,M) Regla compuesta de Simpson

# Ejemplo carga s=simprl(0,pi/3,6)

# funcion 'sin(x)'
# - a Limite inferior de integración
# - b Limite superior de integración
# - M Numero de subintervalos
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');

Página 74 de 115
 La solución numérica elemental con Octave
h=(b-a)/(2*M);
s1=0;
s2=0;
for k=1:M
x=a+h*(2*k-1);
s1=s1+feval(f,x);
end
for k=1:(M-1)
x=a+h*2*k;
s2=s2+feval(f,x);
end
s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;

Hallar la integral de la función f(x)=sen x entre 0 y π/3

Simpon’s compuesta
octave:30> simprl(0,pi/3,6)
ingrese f(x) entre comillas en términos de x 'sin(x)'
'sin(x)'
ans = 0.50000

function [R,quad,err,h]=romber(a,b,n,tol) Integración compuesta de Romberg

%Ingreso - f es el integrando entrado como string 'f'

% - a y b son los límites de integración
% - n es el número máximo de filas en la tabla
% - tol is the tolerance
%salida - R es la tabla de Romberg
% - quad es el valor de cuadratura
% - err es la estimación de error
% - h es el menor tamaño de paso usado
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;

while((err>tol)&(J<n))|(J<4)
J=J+1;
h=h/2;
s=0;
for p=1:M
x=a+h*(2*p-1);
s=s+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*s;
M=2*M;
for K=1:J
R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
Página 75 de 115
 La solución numérica elemental con Octave
end
err=abs(R(J,J)-R(J+1,K+1));
end

Hallar la integral de la función f(x)=sen x entre 0 y π/3

Romberg compuesta
octave:31> romber(0,pi/3,6,0.001)
ingrese f(x) entre comillas en términos de x 'sin(x)'
'sin(x)'
ans =
0.45345 0.00000 0.00000 0.00000 0.00000
0.48852 0.50022 0.00000 0.00000 0.00000
0.49714 0.50001 0.50000 0.00000 0.00000
0.49929 0.50000 0.50000 0.50000 0.00000
0.49982 0.50000 0.50000 0.50000 0.50000

function Ih = trapezoidal(func,a,b,I2h,k) Regla trapezoidal recursiva

% regal trapezoidal recursiva.

% USAR: Ih = trapezoid(func,a,b,I2h,k)
% func = handle de function a integrar.
% a,b = límites de integración.
% I2h = integral con 2ˆ(k-1) paneles.
% Ih = integral con 2ˆk paneles.
if k == 1
fa = feval(func,a); fb = feval(func,b);
Ih = (fa + fb)*(b - a)/2.0;
else
n = 2^(k -2 ); % Número de nuevos puntos
h = (b - a)/n ; % espaciado de nuevos puntos
x = a + h/2.0; % Coord. Del primer punto nuevo
suma = 0.0;
for i = 1:n
fx = feval(func,x);
suma = suma + fx;
x = x + h;
end
Ih = (I2h + h*suma)/2.0;
end

Hallar la integral de la función f(x)=sen x entre 0 y π/3

trapezoidal recursive
octave:32> f=@(x)sin(x);
octave:33> Ih = trapezoidal(f,pi/3,6,1,0)
Ih = 0.50000

Página 76 de 115
 La solución numérica elemental con Octave

function I = gausslegendre(func,a,b,n) Fórmula de Gauss-Legendre

% Cuadratura de Gauss-Legendre.
% USAR I = gaussQuad(func,a,b,n)
% Ingreso:
% func = handle de la function a integrar.
% a,b = límites de integración.
% n = orden de integración.
% Salida:
% I = integral
c1 = (b + a)/2; c2 = (b - a)/2; % constants de mapeo
[x,A] = gaussNodos(n); % pesos y abcisas nodales
sum = 0;
for i = 1:length(x)
y = feval(func,c1 + c2*x(i)); % Función en nodo i
sum = sum + A(i)*y;
end
I = c2*sum;

function [x,A] = gaussNodos(n,tol)

% Calcula abcisas nodales x y pesos A de
% cuadratura del punto n de Gauss-Legendre.
% USAR: [x,A] = gaussNodes(n,epsilon,maxIter)
% tol = tolerancia de error (default = 1.0e4*eps).
if nargin < 2; tol = 1.0e4*eps; end
A = zeros(n,1); x = zeros(n,1);
nRoots = fix(n + 1)/2; % Número de raíces no-neg.
for i = 1:n
t = cos(pi*(i - 0.25)/(n + 0.5)); % raíces aproximadas
for j = i:30
[p,dp] = legendre(t,n); % método de Newton
dt = -p/dp; t = t + dt; % búsqueda de raíz
if abs(dt) < tol %
x(i) = t; x(n-i+1) = -t;
A(i) = 2/(1-t^2)/dp^2; %
A(n-i+1) = A(i);
break
end
end
end

function [p,dp] = legendre(t,n)

% Evalua polinomio de Legendre p de grado n
% y su derivada dp en x = t.
p0 = 1.0; p1 = t;
for k = 1:n-1
p = ((2*k + 1)*t*p1 - k*p0)/(k + 1); % Eq. (6.19)
p0 = p1;p1 = p;
end
dp = n *(p0 - t*p1)/(1 - t^2);

Hallar la integral de la función f(x)=sen x entre 0 y π/3

Gauss-Legendre
octave:35> f=@(x)sin(x);
octave:36> I = gausslegendre(f,0,pi/3,25)
I = 0.50000

Página 77 de 115
 La solución numérica elemental con Octave

function [SRmat,quad,err]=adapt(a,b,tol) Cuadratura adaptativa

%Ingresos - f es el integrando como string 'f'
% - a y b límites inferior y superior de integración
% - tol es la tolerancia
%salida - SRmat tabla de valores
% - quad valor de cuadratura
% - err error estimado
%Initialize values
f=input('ingrese f(x) entre comillas en términos de x');
f=inline(f,'x');
SRmat = zeros(30,6);
iterating=0;
done=1;
SRvec=zeros(1,6);
SRvec=srule(f,a,b,tol);
SRmat(1,1:6)=SRvec;
m=1;
state=iterating;

while(state==iterating)
n=m;
for j=n:-1:1
p=j;
SR0vec=SRmat(p,:);
err=SR0vec(5);
tol=SR0vec(6);
if (tol<=err)

%intervalo Bisecccionado,aplica regla de Simpson

%recursivamente, y determina error
state=done;
SR1vec=SR0vec;
SR2vec=SR0vec;
a=SR0vec(1);
b=SR0vec(2);
c=(a+b)/2;
err=SR0vec(5);
tol=SR0vec(6);
tol2=tol/2;
SR1vec=srule(f,a,c,tol2);
SR2vec=srule(f,c,b,tol2);
err=abs(SR0vec(3)-SR1vec(3)-SR2vec(3))/10;

%prueba de exactitud
if (err<tol)
SRmat(p,:)=SR0vec;
SRmat(p,4)=SR1vec(3)+SR2vec(3);
SRmat(p,5)=err;
else
SRmat(p+1:m+1,:)=SRmat(p:m,:);
m=m+1;
SRmat(p,:)=SR1vec;
SRmat(p+1,:)=SR2vec;
state=iterating;
end
end
end
end
quad=sum(SRmat(:,4));
err=sum(abs(SRmat(:,5)));
SRmat=SRmat(1:m,1:6);
Página 78 de 115
 La solución numérica elemental con Octave

function Z=srule(f,a0,b0,tol0) Regla de Simpson’s

%Ingreso - f es el integrando como string 'f'

% - a0 y b0 límites inferior y superior de integración
% - tol0 es la tolerancia
% salida - Z es un vector 1 x 6 [a0 b0 S S2 err tol1]
h=(b0-a0)/2;
C=zeros(1,3);
C=feval(f,[a0 (a0+b0)/2 b0]);
S=h*(C(1)+4*C(2)+C(3))/3;
S2=S;
tol1=tol0;
err=tol0;
Z=[a0 b0 S S2 err tol1];

Hallar la integral de la función f(x)=sen x entre 0 y π/3

Cuadratura adaptativa de Simpson’s
octave:37> adapt(0,pi/3,0.0001)
ingrese f(x) entre comillas en términos de x 'sin(x)'
'sin(x)'
ans = 0.00000 1.04720 0.50022 0.50001 0.00002 0.00010

function I= gauss_q(f,a,b,n) Cuadratura de Gauss-Legendre

p=Legen_pw(n);
x=roots(p)';
x=sort(x);
for j=1:n
y=zeros(1,n);
y(j)=1;
p=polyfit(x,y,n-1);
P=poly_itg(p);
w(j)=polyval(P,1)- polyval(P,-1);
end
x=0.5*((b-a)*x+a+b);
y= feval(f,x);
I=sum(w.*y)*(b-a)/2;
fprintf('\n x y w\n')
for j=1:n
fprintf('%e %e %e\n', x(j) , y(j) , w(j))
end

function pn=legen_pw(n)
pbb=[1]; if n==0, pn=pbb; break; end
pb=[1 0]; if n==1, pn=pb; break; end
for i=2:n;
pn=((2*i-1)*[pb,0]-(i-1)*[0, 0, pbb])/i;
pbb=pb; pb=pn;
end

Página 79 de 115
 La solución numérica elemental con Octave

function py = poly_itg(p)
% poly_itg(p) integra un polinomio p que es una
% serie de potencias. El resultado también es una serie de potencias.
n=length(p);
py = [p.*[n:-1:1].^(-1),0];

Hallar la integral de la función f(x)=sen x entre 0 y π/3

Cuadratura de Gauss-Legendre
octave:37> I= gauss_q(@(x)sin(x),0,pi/3,20)
x y w
3.597857e-003 3.597849e-003 1.761401e-002
1.886425e-002 1.886314e-002 4.060143e-002
4.595395e-002 4.593777e-002 6.267205e-002
8.423816e-002 8.413857e-002 8.327674e-002
1.328203e-001 1.324301e-001 1.019301e-001
1.905618e-001 1.894106e-001 1.181945e-001
2.561094e-001 2.533188e-001 1.316886e-001
3.279267e-001 3.220809e-001 1.420961e-001
4.043304e-001 3.934032e-001 1.491730e-001
4.835296e-001 4.649070e-001 1.527534e-001
5.636680e-001 5.342903e-001 1.527534e-001
6.428672e-001 5.994927e-001 1.491730e-001
7.192708e-001 6.588363e-001 1.420961e-001
7.910881e-001 7.111187e-001 1.316886e-001
8.566357e-001 7.556433e-001 1.181945e-001
9.143772e-001 7.921827e-001 1.019301e-001
9.629594e-001 8.208853e-001 8.327674e-002
1.001244e+000 8.421423e-001 6.267205e-002
1.028333e+000 8.564397e-001 4.060143e-002
1.043600e+000 8.642209e-001 1.761401e-002

I =

0.5000

function R=rk4(a,b,ya,M) Métodos de Runge-Kutta

%Ingresos - f es la función como string 'f'

% - a y b extremos derecho e izquierdo
% - ya condición inicial y(a)
% - M es el número de etapas
%salida - R = [T' Y'] con T vector de abcisas
% e Y vector de ordenadas

f=input('ingrese f entre comillas en términos de t,y');

f=inline(f,'t','y');
h=(b-a)/M;
T=zeros(1,M+1);
Y=zeros(1,M+1);
T=a:h:b;

Página 80 de 115
 La solución numérica elemental con Octave
Y(1)=ya;
for j=1:M
k1=h*feval(f,T(j),Y(j));
k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
k4=h*feval(f,T(j)+h,Y(j)+k3);
Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end

R=[T' Y'];

Sea la ecuación y’=0.5*(t-y) con y(0)=1 en [0,3], se plantean los métodos de

Adams Basforth-Moulton( 4 pasos);Milne- Simpson ( 3 pasos), Hamming(3
pasos); todos requieren el cálculo previo de coordenadas T e Y ( por ej con
Runge Kutta de cuarto orden)
Seleccionando un paso de 0.25
octave:10> R=rk4(0,3,1,12)
ingrese f entre comillas en términos de t,y '0.5*(t-y)'
'0.5*(t-y)'
R =
0.00000 1.00000
0.25000 0.89749
0.50000 0.83640
0.75000 0.81187
1.00000 0.81959
1.25000 0.85579
1.50000 0.91710
1.75000 1.00059
2.00000 1.10364
2.25000 1.22396
2.50000 1.35952
2.75000 1.50852
3.00000 1.66939

function A=abm(T,Y) Métodos de Adams-Moulton

%Ingreso - f función como string 'f'

% - T vector de abcisas
% - Y vector de ordenadas
% Las primeras cuatro coordendaas de T e Y deben
% tener valores iniciales obtenidos con RK4
%salida - A=[T' Y'] con T vector de abcisas e
% Y vector de ordenadas

f=input('ingrese f entre comillas en términos de t,y');

f=inline(f,'t','y');
n=length(T);

Página 81 de 115
 La solución numérica elemental con Octave
if n<5, return, end;

F=zeros(1,4);
F=feval(f,T(1:4),Y(1:4));
h=T(2)-T(1);

for k=4:n-1
%Predictor
p=Y(k)+(h/24)*(F*[-9 37 -59 55]');
T(k+1)=T(1)+h*k;
F=[F(2) F(3) F(4) feval(f,T(k+1),p)];
%Corrector
Y(k+1)=Y(k)+(h/24)*(F*[1 -5 19 9]');
F(4)=feval(f,T(k+1),Y(k+1));
end

A=[T' Y'];

Empleando los datos optenidos con el algoritmo rk4.

Desde la ventana de comandos hacemos:
octave:11> T=zeros(1,13);
octave:12> Y=zeros(1,13);
octave:13> T=0:1/4:3;
octave:14> Y(1:4)=[1 0.8975 0.8364 0.8119];
octave:15> abm(T,Y)
ingrese f entre comillas en términos de t,y '0.5*(t-y) '
'0.5*(t-y) '
ans =
0.00000 1.00000
0.25000 0.89750
0.50000 0.83640
0.75000 0.81190
1.00000 0.81962
1.25000 0.85580
1.50000 0.91711
1.75000 1.00060
2.00000 1.10365
2.25000 1.22396
2.50000 1.35952
2.75000 1.50852
3.00000 1.66939

function M=milne(T,Y) Método de Milne

%Ingreso - f función como 'f'

% - T vector de abcisas
% - Y vector de ordenadas
%. Las primeras cuatro coordendaas de T e Y deben
% tener valores iniciales obtenidos con RK4

%salida - M=[T' Y'] con T vector de abcisas e

Página 82 de 115
 La solución numérica elemental con Octave
% Y vector de ordenadas

f=input('ingrese f entre comillas en términos de t,y');

f=inline(f,'t','y');
n=length(T);
if n<5, break, end;

F=zeros(1,4);
F=feval(f,T(1:4),Y(1:4));
h=T(2)-T(1);
pold=0;
yold=0;

for k=4:n-1
%Predictor
pnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');
%Modifier
pmod=pnew+28*(yold-pold)/29;
T(k+1)=T(1)+h*k;
F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];
%Corrector
Y(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]');
pold=pnew;
yold=Y(k+1);
F(4)=feval(f,T(k+1),Y(k+1));
end

Empleando los datos optenidos con el algoritmo rk4.

Desde la ventana de comandos hacemos:
octave:11> T=zeros(1,13);
octave:12> Y=zeros(1,13);
octave:13> T=0:1/4:3;
octave:14> Y(1:4)=[1 0.8975 0.8364 0.8119];

octave:15> milne(T,Y)
ingrese f entre comillas en términos de t,y '0.5*(t-y) '
'0.5*(t-y) '
ans =
0.00000 1.00000
0.25000 0.89750
0.50000 0.83640
0.75000 0.81190
1.00000 0.81958
1.25000 0.85582
1.50000 0.91709
1.75000 1.00062
2.00000 1.10362
2.25000 1.22399
2.50000 1.35950
2.75000 1.50855
3.00000 1.66937

Página 83 de 115
 La solución numérica elemental con Octave

function H=heun(a,b,ya,M) Método Predictor-Corrector

# - a y b son los puntos finales de izquierda y derecha

# - ya condicion inicial y(a)
# - M Tamaño de paso
#Salidas - H=[T' Y'] T: vector de abcisas Y: vector de ordenadas
# Ejemplo carga H=heun(0,1,1,5)
# funcion '3*y+3*t'
f=input('ingrese f entre comillas en términos de t,y');
f=inline(f,'t','y');
h=(b-a)/M;
T=zeros(1,M+1);
Y=zeros(1,M+1);
T=a:h:b;
Y(1)=ya;
for j=1:M
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j+1),Y(j)+h*k1);
Y(j+1)=Y(j)+(h/2)*(k1+k2);
end
H=[T' Y'];

Dada la ecuación diferencial en PVI

y' =3t+3y
y(0) = 1 , en [0,1]
Solución Exacta: -1/3-t+4/3*exp(3*t)
H1 = heun(0,1,1,5)
ingrese f entre comillas en términos de t,y '3*y+3*t'
'3*y+3*t'
H1 =
0.00000 1.00000
0.20000 1.84000
0.40000 3.49120
0.60000 6.58634
0.80000 12.25168
1.00000 22.49199

function X = backsub(A,B) SUSTITUCIÓN HACIA ATRÁS

# --------
# Llamada
# X = backsub(A,B)
# Parametros
# A Matriz de coeficientes triangular
# superior obtenida de gauss(A)
# B Vector lado derecho de la ecuacion
# Devuelve
# X Vector de Solucion
#Ejemplo carga X = backsub([4 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3],[20;-
7;4;6])
n = length(B);
det1 = A(n,n);
X = zeros(n,1);
X(n) = B(n)/A(n,n);

Página 84 de 115
 La solución numérica elemental con Octave
for r = n-1:-1:1,
det1 = det1*A(r,r);
if det1 == 0, break, end
X(r) = (B(r) - A(r,r+1:n)*X(r+1:n))/A(r,r);
end
end

Resolver un sistema triangular (cuadrada), dada la matriz de coefcientes A y

el vector independiente B
octave:17> A=[4 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3];
octave:18> B=[20;-7;4;6];
octave:19> backsub(A,B)
ans =
3
-4
-1
2

function X = uptrbk(A,B) SUSTITUCIÓN HACIA ATRÁS

%Ingreso - A matriz N x N no singular

% - B N x 1 matriz
%salida - X es matriz N x 1 con la solución de AX=B.

%Iniciaalizar X y la matriz temporario de almacenaje C

[N N]=size(A);
X=zeros(N,1);
C=zeros(1,N+1);

%Forma la matriz aumentada: Aug=[A|B]

Aug=[A B];
for p=1:N-1
%Pivoteo parcial para columna p
[Y,j]=max(abs(Aug(p:N,p)));
%Intercambio fila p y j
C=Aug(p,:);
Aug(p,:)=Aug(j+p-1,:);
Aug(j+p-1,:)=C;

if Aug(p,p)==0
'A es singular. No hay solución única'
break
end
%etapa de eliminación para columna p
for k=p+1:N
m=Aug(k,p)/Aug(p,p);
Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);
end
end
%sustituación hacia atrás en [U|Y]
X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));

Página 85 de 115
 La solución numérica elemental con Octave
Resolver un sistema triangular (cuadrada), dada la matriz de coefcientes A y
el vector independiente B , empleando triangularización superior más
sustitución hacia atrás
octave:22> A=[1 2 2 4;2 0 4 3;4 2 2 1;-3 1 3 2]; % no singular
octave:23> B=[13 28 20 6]';
octave:24> uptrbk(A,B)
ans =
3.00000
-1.33333
5.00000
0.66667

function Gjelim(P) TÉCNICA DE GAUSS-JORDAN

%Gjelim
%eliminación de Gauss-Jordan.
%Opciones: formato en numero racional
% conteo de operaciones
% todas las etapas
%forma de llamado: Gjelim(A)

%puede usarse tolerancia de 1e-20

%o cambiar a una propia
%MATLAB calcula hasta 16 dígitos decimales
%
%valores iniciales
adds=0;totadds=0;mults=0;totmults=0;swaps=0;totswaps=0;
ops=[0];

%hold off
%
tol=1e-20;
[n,m]=size(P);
format compact

disp(' ')
h=input('Números racionales? y/n: ','s');
q=input('Conteo de operaciones? y/n: ','s');
g=input('todas las etapas? y/n: ','s');

disp(' ')
disp('matriz inicial')
if h=='y'
disp(rats(P))
else
disp(P)
end

if g=='y';
disp('[presione Enter en cada paso para seguir]')
disp(' ')
pause
end
Página 86 de 115
 La solución numérica elemental con Octave

%busca un pivot
j=1;
for i=1:n,
if j <= m
found=0;
if abs(P(i, j)) <= tol %fin está en línea 101
while (found == 0) %
%busca de unos e intercambio de filas si hace falta
for s=i:n,
if (abs(P(s, j)) > tol)
if (found == 0)
found=1;
if s~=i
for r=j:m,
temp=P(i, r);
P(i, r)=P(s, r);
P(s, r) = temp;
end
swaps = m-j+1;
totswaps = totswaps + swaps;
if g=='y'; %todas las etapas
disp('intercambio filas')
if h=='y'
disp(rats(P))
else
disp(P)
end
if q=='y'
disp('intercambio de elmentos:')
disp(swaps)
end
disp('--------------------')
pause
end %todas las etapas
end
end
end
end

if (found==0)
if (j <= m)
j = j + 1;
end
end

if j>m
found=1;
end
end %
if j > m
found = 0;
end
else
found = 1;
end %parte línea 51

%normaliza elemento lider en la fika ajustándo el resto

if found == 1
k=i;
if (P(k, j) ~= 1)
Página 87 de 115
 La solución numérica elemental con Octave
if (abs(P(k, j)) > tol)
y = P(i, j);
for l=j:m,
P(k, l) = P(k, l)/y ;
end
mults = m-j;
totmults = totmults + mults;
if g=='y'; %todas las etapas
disp('normaliza')
if h=='y'
disp(rats(P))
else
disp(P)
end
if q=='y'
disp('multiplicaciones:')
disp(mults)
end
disp('--------------------')
pause
end %todas las etapas
end
end
for r=1:n,
if (abs(P(r, j)) >tol)
if (r ~= i)
z=P(r, j);
for c=j:m,
P(r, c)=P(r, c) - z * P(i, c);
end
adds = m-j;
mults = m-j;
totadds = totadds + adds;
totmults = totmults + mults;
if g=='y'; %todas las etapas
disp('crea cero')
if h=='y'
disp(rats(P))
else
disp(P)
end
if q=='y'
disp('sumas, multiplicaciones:')
ops=[adds mults];
disp(ops)
end
disp('--------------------')
pause
end %todas las etapas
end
end
end
end

j = j + 1;

end
end %fin loop i

%muestra matriz final

disp('-forma escalonada reducida-')

Página 88 de 115
 La solución numérica elemental con Octave

if h=='y'
disp(rats(P))
else
disp(P)
end

if q=='y'
disp('Total de sumas, multiplicaciones, intercambio elementos:')
ops=[totadds totmults totswaps];
disp(ops)
disp('--------------------')
end

disp(' ')
format loose

function x = gaussPiv(A,b) PIVOTEO

% resuelve A*x = b por eliminación de Gauss con pivoteo de filas
% USAR: x = gaussPiv(A,b)
if size(b,2) > 1; b = b'; end
n = length(b); s = zeros(n,1);
%----------establce arreglo factor escala----------
for i = 1:n; s(i) = max(abs(A(i,1:n))); end
%---------intercambio filas si es necesario----------
for k = 1:n-1
[Amax,p] = max(abs(A(k:n,k))./s(k:n));
p = p + k - 1;
if Amax < eps; error('Matriz es singular'); end
if p ~= k b = swapRows(b,k,p); s = swapRows(s,k,p); A = swapRows(A,k,p);
end
%--------------Eliminación---------------
for i = k+1:n
if A(i,k) ~= 0
lambda = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - lambda*A(k,k+1:n);
b(i) = b(i) - lambda*b(k);
end
end
end
%------------fase de sustitución hacia atrás----------
for k = n:-1:1
b(k) = (b(k) - A(k,k+1:n)*b(k+1:n))/A(k,k);
end
x = b;

function v = swapRows(v,i,j)
% Swap rows i and j of vector or matrix v.
% USAGE: v = swapRows(v,i,j)
temp = v(i,:);
v(i,:) = v(j,:);
v(j,:) = temp;

Página 89 de 115
 La solución numérica elemental con Octave

resuelve A*x = b por eliminación de Gauss con pivoteo de filas, dada la

matriz de coefcientes A y el vector independiente B , empleando
triangularización superior más sustitución hacia atrás

octave:49> A=[4 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3];

octave:50> B=[20;-7;4;6];
octave:51> x = gaussPiv(A,B)
x =
3
-4
-1
2

function X = lufact(A,B) Factorización LU

%Ingreso - A es matriz N x N
% - B matriz N x 1
%salida - X matriz N x 1 conteniendo solución de AX = B.

%Inicializar X, Y,matriz temporario de almacenaje C, y la matriz de

información de permutación fila R

[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
C=zeros(1,N);
R=1:N;

for p=1:N-1

%Busca fila pivot para columna p

[max1,j]=max(abs(A(p:N,p)));

%Intercambia fila p y j
C=A(p,:);
A(p,:)=A(j+p-1,:);
A(j+p-1,:)=C;
d=R(p);
R(p)=R(j+p-1);
R(j+p-1)=d;

if A(p,p)==0
'A es singular. Sin solución única'
break
end

%Calcula multiplicadores y ubica en la subdiagonal de A

for k=p+1:N
mult=A(k,p)/A(p,p);
Página 90 de 115
 La solución numérica elemental con Octave
A(k,p) = mult;
A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);
end
end

%resuelve para Y
Y(1) = B(R(1));
for k=2:N
Y(k)= B(R(k))-A(k,1:k-1)*Y(1:k-1);
end

%resuelve para X
X(N)=Y(N)/A(N,N);
for k=N-1:-1:1
X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end

Resolver un sistema triangular (cuadrada), dada la matriz de coefcientes A y

el vector independiente B, con factorización con pivoteo
octave:49> A=[4 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3];
octave:50> B=[20;-7;4;6];
octave:55> X = lufact(A,B)
X =
3.00000
-1.33333
5.00000
0.66667

function [P,dP,Z] = jacobi(A,B,P,delta,max1) Método de Jacobi

%JACOBI iteración de Jacobi para resolver un sistema lineal.

% forma de llamado
% [X,dX] = jacobi(A,B,P,delta,max1)
% [X,dX,Z] = jacobi(A,B,P,delta,max1)
% Ingresos
% A matriz de coeficientes
% B vector lado derecho
% P vector de arranque
% delta telerancia de convergencia
% max1 maximo número de iteraciones
% devuelve
% X vector solución
% dX vector de estimación error
% Z matriz de historia de iteraciones
%
Z = P';
n = length(B);
Pnew = P;
for k=1:max1,
for r = 1:n,
Sum1 = B(r) - A(r,[1:r-1,r+1:n])*P([1:r-1,r+1:n]);
Pnew(r) = Sum1/A(r,r);
Página 91 de 115
 La solución numérica elemental con Octave
end
dP = abs(Pnew-P);
err = norm(dP);
relerr = err/(norm(Pnew)+eps);
P = Pnew;
Z = [Z;P'];
if (err<delta)|(relerr<delta), break, end
end

Resolución por técnicas iterativas , dada la matriz de coeficientes y el

vector independiente
octave:35> A=[1 -5 -1;4 1 -1;2 -1 -6];
octave:36> B=[-8 13 -2]';
octave:37> P=[0 0 0]';
octave:38> jacobi(A,B,P,0.001, 10) con 10 iteraciones
ans =
9.1711e+006
5.8774e+006
-8.0714e+005
octave:39> jacobi(A,B,P,0.001, 20) con 20 iteraciones
ans =
-3.0646e+013
-1.4564e+013
2.2794e+012

function [P,dP,Z] = gseid(A,B,P,delta,max1) Método de Gauss-Seidel

# GSEID Iteracion para resolver sistemas lineales Gauss-Seidel

# Llamado de la funcion
# [X,dX] = gseid(A,B,P,delta,max1)
# [X,dX,Z] = gseid(A,B,P,delta,max1)
# Entradas
# A Matriz de coeficientes
# B vector de soluciones
# P vector de inicio
# delta tolerancia
# max1 maximo numero de iteraciones
# Devuelve
# P vector solución
# dP vector de estimación error
# Z matriz de historia de las iteraciones
# Ejemplo carga [P,dP,Z] = gseid([1 2 3;2 4 3;3 6 4],[-
1;3;0],[0;0;0],0.0001,20)
Z = P';
n = length(B);
Pold = P;
for k=1:max1,
for r = 1:n,
Sum1 = B(r) - A(r,[1:r-1,r+1:n])*P([1:r-1,r+1:n]);
P(r) = Sum1/A(r,r);

Página 92 de 115
 La solución numérica elemental con Octave
end
dP = abs(Pold-P);
err = norm(dP);
relerr = err/(norm(P)+eps);
Pold = P;
Z = [Z;P'];
if (err<delta)|(relerr<delta), break, end
end
end

Resolución por técnicas iterativas , dada la matriz de coeficientes y el

vector independiente
octave:40> A=[1 -5 -1;4 1 -1;2 -1 -6];
octave:41> B=[-8 13 -2]';
octave:42> P=[0 0 0]';
octave:43> gseid(A,B,P,0.001,10)
ans =
3.8556e+012
-1.5624e+013
3.8892e+012
function x = conjgrad(A,b,tol) GRADIENTE CONJUGADO

% CONJGRAD Metodo del gradiente conjugado.

% X = CONJGRAD(A,B) busca resolver el sistema lineal A*X=B
% para X. La matriz de coeficientes ,NxN,A debe ser simétrica y el vector
% columna B de longitud N.
%
% X = CONJGRAD(A,B,TOL) especifica la tolerancia del método, por
% default es 1e-10.
%
% Ejemplo:
%{
n = 6000;
m = 8000;
A = randn(n,m);
A = A * A';
b = randn(n,1);
tic, x = conjgrad(A,b); toc
norm(A*x-b)
%}

if nargin<3
tol=1e-10;
end
r = b + A*b;
y = -r;
z = A*y;
s = y'*z;
t = (r'*y)/s;
x = -b + t*y;

for k = 1:numel(b);
r = r - t*z;
if( norm(r) < tol )
return;
end
B = (r'*z)/s;
y = -r + B*y;
z = A*y;
s = y'*z;

Página 93 de 115
 La solución numérica elemental con Octave
t = (r'*y)/s;
x = x + t*y;
end
end

Utilizar el método del Gradiente conjugado( directo no es tan bueno como

elim. Gauss con pivoteo, sí para iterativos)
octave:44> A=[4 3 0;3 4 -1;0 -1 4]; % debe ser simétrica
octave:45> b=[24 30 -24]';
octave:46> conjgrad(A,b,0.001)
ans =
3.0000
4.0000
-5.0000

%El método del gradiente conjugado ayuda a resolve el sistema Ax=b, con A
simétrica, sin calcula inversa de A. Sólo requiere poca memoria, siendo
adecuado para sistemas de gran escala.

%es más rápido que otros como eliminación Gaussiana, si A está bien
condicionada. Por ejemplo,

%n=1000;
%[U,S,V]=svd(randn(n));
%s=diag(S);
%A=U*diag(s+max(s))*U'; % A es simétrica, bien condicionada
%b=randn(1000,1);
%tic,x=conjgrad(A,b);toc
%tic,x1=A\b;toc
%norm(x-x1)
%norm(x-A*b)

%El GC es de dos a tres veces más rápido que A\b, que emplea eliminación
Gaussiana

function [C]=gerschgorin(A) VALORES PROPIOS

% Gershgorin's circulos C of the matrix A.

d = diag(A);
cx = real(d);
cy = imag(d);
B = A - diag(d);
[m, n] = size(A);
r = sum(abs(B'));
C = [cx cy r(:)];
t = 0:pi/100:2*pi;
c = cos(t);
s = sin(t);
[v,d] = eig(A);
d = diag(d);
u1 = real(d);
v1 = imag(d);
hold on
Página 94 de 115
 La solución numérica elemental con Octave
grid on
axis equal
xlabel('Re')
ylabel('Im')
h1_line = plot(u1,v1,'or');
set(h1_line,'LineWidth',1.5)
for i=1:n
x = zeros(1,length(t));
y = zeros(1,length(t));
x = cx(i) + r(i)*c;
y = cy(i) + r(i)*s;
h2_line = plot(x,y);
set(h2_line,'LineWidth',1.2)
end
hold off
title('Gershgorin circles and the eigenvalues of a')

Para A, acotar los autovalores empleando los círculos de Gerschgorin

octave:59> [C] = gerschgorin([4 1 1;0 2 1;-2 0 9])
C =
4 0 2
2 0 1
9 0 2

function [lambda,V]= power1(A,X,epsilon,max1) Métodos de Potencias

%Ingreso - A matriz nxn

% - X vector de arranque nx1
% - epsilon es la tolerancia
% - max1 número máximo de iteraciones
%salida - lambda eigenvalue dominante
% - V eigenvector dominante

%Inicializar parámetros

lambda=0;
cnt=0;
err=1;
state=1;

while ((cnt<=max1)&(state==1))
Y=A*X;

%Normaliza Y
[m j]=max(abs(Y));
c1=m;
dc=abs(lambda-c1);
Y=(1/c1)*Y;

%actualiza X y lambda y chequea para convergencia

dv=norm(X-Y);
Página 95 de 115
 La solución numérica elemental con Octave
err=max(dc,dv);
X=Y;
lambda=c1;
state=0;
if (err>epsilon)
state=1;
end
cnt=cnt+1;
end

V=X;

Dada A=[0 11 -5;-2 17 -7;-4 26 -10]; hallar el autovalor y autovector

dominante
octave:58> A=[0 11 -5;-2 17 -7;-4 26 -10];
octave:59> X=[1 1 1]';
octave:60> power1(A,X,0.001,10)
ans = 4.0016

function[lambda,V]=invpow(A,X,alpha,epsilon,max1) Método de la potencia inversa

%Ingreso - A matriz nxn

% - X vector de arranque nx1
% - alpha dado
% - epsilon tolerancia
% - max1 número máximo de iteraciones
%salida - lambda es el eigenvalue dominante
% - V eigenvector dominante

%Inicializar matriz A-alphaI y parámetros

[n n]=size(A);
A=A-alpha*eye(n);
lambda=0;
cnt=0;
err=1;
state=1;

while ((cnt<=max1)&(state==1))

%resuelve sistema AY=X

Y=A\X;

%Normaliza Y
[m j]=max(abs(Y));
c1=m;
dc=abs(lambda-c1);
Y=(1/c1)*Y;
Página 96 de 115
 La solución numérica elemental con Octave

%actualiza X y lambda y chequea la convergencia

dv=norm(X-Y);
err=max(dc,dv);
X=Y;
lambda=c1;
state=0;
if (err>epsilon)
state=1;
end
cnt=cnt+1;
end

lambda=alpha+1/c1;
V=X;

Dada A=[0 11 -5;-2 17 -7;-4 26 -10]; hallar el autovalor y autovector

dominante
Por el método de inverso de potencias
octave:61> A=[0 11 -5;-2 17 -7;-4 26 -10];
octave:62> X=[1 1 1]';
octave:63> alpha=4.1;
octave:64> invpow(A,X,alpha,0.001,10)
ans = 4.2000

function [V,D]=jacobi1(A,epsilon) METODO CLASICO DE JACOBI

%Ingreso - A matriz nxn

% - epsilon tolerancia
%salida - V matriz nxn de eigenvectors
% - D matriz diagonal nxn de eigenvalues
%Inicializar V, D, y parámetros
D=A;
[n,n]=size(A);
V=eye(n);
state=1;

%Calcula fila p y columna q del elementodiagonal más grande deelement

% A

[m1 p]=max(abs(D-diag(diag(D))));
[m2 q]=max(m1);
p=p(q);

while (state==1)
%elim. cero Dpq y Dqp
t=D(p,q)/(D(q,q)-D(p,p));
c=1/sqrt(t^2+1);
s=c*t;
R=[c s;-s c];
D([p q],:)=R'*D([p q],:);
D(:,[p q])=D(:,[p q])*R;
V(:,[p q])=V(:,[p q])*R;
[m1 p]=max(abs(D-diag(diag(D))));
[m2 q]=max(m1);
p=p(q);

Página 97 de 115
 La solución numérica elemental con Octave
if (abs(D(p,q))<epsilon*sqrt(sum(diag(D).^2)/n))
state=0;
end
end
D=diag(diag(D));

Para A simétrica, hallar sus autovalores y autovectores

octave:65> A=[8 -1 3 -1; -1 6 2 0; 3 2 9 1; -1 0 1 7];
octave:66> jacobi1(A,0.001)
ans =
0.528779 -0.572958 0.582194 0.230569
0.591968 0.472072 0.176060 -0.629067
-0.536039 0.282024 0.792519 -0.070981
0.287453 0.607725 0.044347 0.738968
La salida es una matriz 4x4 de eigenvectores y diagonal de eigenvalores

function T=house(A) Técnica de Householder

%Ingreso - A matriz simétrica nxn

%salida - T matriz tridiagonal

[n,n]=size(A);

for k=1:n-2
s=norm(A(k+1:n,k));
if (A(k+1,k)<0)
s=-s;
end
r=sqrt(2*s*(A(k+1,k)+s));
W(1:k)=zeros(1,k);
W(k+1)=(A(k+1,k)+s)/r;
W(k+2:n)=A(k+2:n,k)'/r;
V(1:k)=zeros(1,k);
V(k+1:n)=A(k+1:n,k+1:n)*W(k+1:n)';
c=W(k+1:n)*V(k+1:n)';
Q(1:k)=zeros(1,k);
Q(k+1:n)=V(k+1:n)-c*W(k+1:n);
A(k+2:n,k)=zeros(n-k-1,1);
A(k,k+2:n)=zeros(1,n-k-1);
A(k+1,k)=-s;
A(k,k+1)=-s;
A(k+1:n,k+1:n)=A(k+1:n,k+1:n) ...
-2*W(k+1:n)'*Q(k+1:n)-2*Q(k+1:n)'*W(k+1:n);
end
T=A;

Página 98 de 115
 La solución numérica elemental con Octave

Reducir una matriz simétrica a la forma tridiagonal(por transformación de

Householder)
octave:67> A=[8 -1 3 -1; -1 6 2 0; 3 2 9 1; -1 0 1 7];
octave:68> house (A)
ans =
8.00000 3.31662 0.00000 0.00000
3.31662 6.90909 -2.46630 0.00000
0.00000 -2.46630 8.24308 -0.79311
0.00000 0.00000 -0.79311 6.84783

function D=qr1(A,epsilon) Algoritmo QR

%Ingreso - A matriz tridiagonal simétrica nxn

% - epsilon tolerancia
%salida - D vector de eigenvalues

%Inicializar parámetros

[n,n]=size(A);
m=n;

while (m>1)
S=A(m-1:m,m-1:m);
if abs(S(2,1))<epsilon
A(m,m-1)=0;
A(m-1,m)=0;

else
shift=eig(S);
[j,k]=min([abs(A(m,m)-shift(1)) abs(A(m,m)-shift(2))]);
end
[Q,U]=qr(A-shift(k)*eye(n));
A=U*Q+shift(k)*eye(n);
m=m-1;
end
D=diag(A);

Aproximar los autovalores con el método QR

octave:69> A=[8 -1 3 -1; -1 6 2 0; 3 2 9 1; -1 0 1 7];
octave:70> qr1(A,0.001)
ans =
3.2976
8.4547
11.6541
6.5936

Página 99 de 115
 La solución numérica elemental con Octave

function D=qr2(A,epsilon) Algoritmo QR

%Ingreso- A matriz simétrica tridiagonal nxn

% - epsilon tolerancia
%salida - D nx1 vector de eigenvalues

%Inicializar parametros

[n,n]=size(A);
m=n;
D=zeros(n,1);
B=A;

while (m>1)
while (abs(B(m,m-1))>=epsilon)

%Calcula shift
S=eig(B(m-1:m,m-1:m));
[j,k]=min([abs(B(m,m)*[1 1]'-S)]);

%factorización QR de B
[Q,U]=qr(B-S(k)*eye(m));

%Calcula siguiente B
B=U*Q+S(k)*eye(m);
end

%ubica el m-simo eigenvalue en A(m,m)

A(1:m,1:m)=B;

%Repite proceso en m-1 x m-1 submatriz de A

m=m-1;
B=A(1:m,1:m);
end
D=diag(A);

Aproximar los autovalores con el método QR

octave:71> A=[8 -1 3 -1; -1 6 2 0; 3 2 9 1; -1 0 1 7];
octave:72> qr2(A,0.001)
ans =
11.7043
3.2957
8.4077
6.5923
Con A tridiagonal simétrica cuadrada

Página 100 de 115

 La solución numérica elemental con Octave

function [Q, R] = grams(A) Ortogonalización de Gram-Schdmidt

% gschmidt.m: Ortogonalización de las columnas de A

% vía el proceso de Gram-Schmidt.
% Se asumen que las columnas de A son LI
%
% Forma de uso:
% Q = grams(A) regresa una martiz Q (de orden m x n) cuyas
% columnas forman una base ortonormal para el espacio de A.
%
% [Q, R] = grams(A) regresa una matriz Q con columnas ortonormales
% y R como una matriz trinagular superior de tal forma que A = Q*R.
%
[m, n] = size(A);
Asave = A;
for j = 1:n
for k = 1:j-1
mult = (A(:, j)'*A(:, k)) / (A(:, k)'*A(:, k));
A(:, j) = A(:, j) - mult*A(:, k);
end
end
for j = 1:n
if norm(A(:, j)) < sqrt(eps)
error('Las columnas de A son linalmente independientes.')
end
Q(:, j) = A(:, j) / norm(A(:, j));
end
R = Q'*Asave;

Dados la matriz A aplicar la ortogonalización de Gram Schdmit

Grams Schmidt
A=[8 -1 3 -1; -1 6 2 0; 3 2 9 1; -1 0 1 7];
octave:75> [Q, R] = grams(A)
Q =
0.923760 -0.023148 -0.322415 0.205377
-0.115470 0.930114 -0.334071 0.099754
0.346410 0.366154 0.841005 -0.196575
-0.115470 -0.016835 0.277769 0.953534
R =
8.6603e+000 -9.2376e-001 5.5426e+000 -1.3856e+000
2.0817e-017 6.3361e+000 5.0693e+000 2.7146e-001
5.5511e-017 6.6613e-016 6.2114e+000 3.1078e+000
-1.1102e-016 -2.7756e-016 -2.2204e-016 6.2728e+000

Página 101 de 115

 La solución numérica elemental con Octave
OBSERVACION: La versión de octave a emplear para el algoritmo de broyden
deberá contar con el paquete “simbols”, dado que sin el mismo no puede
operarse con variables simbólicas.

function [x,k] = broyden(fcn1,fcn2,x0,maxits,tol) CUASI-NEWTON

% sintaxis: [X] = BROYDEN(FCN1,FCN2,X0,MAXITS)
% Una función para calcular un sistema 2x2 de ec. nolineales
% los ingresos mínimos requeridos son fcn1 y fcn2
% tener cuidado al ingresar fcn1,fcn2 please usar sólo x1 y x2
%
%
% ingresos - FCN1: primera ecuación como string
% - FCN2: segunda ecuación como string
% - X0: aproximation inicial de la solución, default es [1 1]
% - MAXITS: máximo número de iteraciones, default es 50
% - TOL: tolerancia, default es 1e-8
% salida - X: the solución del sistema
% - K: número de iteraciones para llegar a la solución
%
%

if nargin < 5
tol = 1e-8;
end
if nargin < 4
maxits = 50;
end
if nargin < 3
x0 = [1 1]';
end
if nargin < 2
help broyden
error('dos entradas como mínimo')
end
syms x1 x2
B = [diff(fcn1,x1) diff(fcn1,x2);diff(fcn2,x1) diff(fcn2,x2)];
x1 = x0(1);
x2 = x0(2);
h = inline(fcn1,'x1','x2');
g = inline(fcn2,'x1','x2');
f(1:2,1) = [h(x1,x2);g(x1,x2)];
B = eval(B);
x = [x1 x2]';
for k=1:maxits
s = B\(-f);
x = x + s;
fnew = [h(x(1),x(2));g(x(1),x(2))];
y = fnew-f;
if abs(fnew-f) < tol
break
end
f = fnew;
B = B + ((y-B*s)*s')/(s'*s);
end

Página 102 de 115

 La solución numérica elemental con Octave

function [x,fx,xx]=newtons(f,x0,TolX,MaxIter,varargin) MÉTODO DE NEWTON

%newtons.m resuelve conjunto de ec. nolineales f1(x)=0, f2(x)=0,..

%ingreso: f = vector1^st-orden ftn equivalente a un conjunto de ecuaciones
% x0 = suposición inicial de solución
% TolX = límite superior de |x(k) - x(k - 1)|
% MaxIter = máximo número de iteraciones
% salida: x = punto donde el algoritmo llegó
% fx = f(x(última))
% xx = historia de x
h = 1e-4; TolFun = eps; EPS = 1e-6;
fx = feval(f,x0,varargin{:});
Nf = length(fx); Nx = length(x0);
if Nf ~= Nx; error('dimensiones incompatible de f y x0!'); end
if nargin < 4, MaxIter = 100; end
if nargin < 3, TolX = EPS; end
xx(1,:) = x0(:).'; %Inicializar la solución como vector fila inicial
%fx0 = norm(fx); %(1)
for k = 1: MaxIter dx = -jacob(f,xx(k,:),h,varargin{:})\fx(:);%-[dfdx]ˆ-1*fx
%para l = 1: 3 %damping para eviatr divergencia %(2)
%dx = dx/2; %(3)
xx(k + 1,:) = xx(k,:) + dx.';
fx = feval(f,xx(k + 1,:),varargin{:}); fxn = norm(fx);
% if fxn < fx0, break; end %(4)
%end %(5)
if fxn < TolFun | norm(dx) < TolX, break; end
%fx0 = fxn; %(6)
end
x = xx(k + 1,:);
if k == MaxIter;
fprintf('la mejor en %d iteraciones\n',MaxIter);
end

function g = jacob(f,x,h,varargin) %Jacobiano de f(x)

if nargin < 3; h = 1e-4; end
h2 = 2*h; N = length(x); x = x(:).'; I = eye(N);
for n = 1:N
g(:,n) = (feval(f,x + I(n,:)*h,varargin{:}) ...
-feval(f,x - I(n,:)*h,varargin{:}))'/h2;
end

%Método del disparo Lineal TECNICA DE DISPARO PARA EL PROBLEMA LINEAL

% - - - - - - - - - - - - - - - - - - - - - - - - -
%
% este programa implementa el meétodo de disparo lineal
%
% para resolver el PVF
%
% x'' = p(t)x'(t) + q(t)x(t) + r(t)
%
% sobre [a,b] con x(a) = alpha y x(b) = beta
%
Página 103 de 115
 La solución numérica elemental con Octave
% Usa el cambio de variable y = x' y las
%
% funciones p(t)x'(t) + q(t)x(t) + r(t)
%
% p(t)x'(t) + q(t)x(t)
%
% para formar fn1.m y fm2.m
% - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
% Aplica el método de Runge-Kutta para un orden má alto para resolver:
%
% u'' = p(t)u'(t) + q(t)u(t) + r(t)
% con u(a) = alpha y u'(a) = 0
%
% v'' = p(t)v'(t) + q(t)v(t)
% con v(a) = 0 y v'(a) = 1
%
% Define y almacena la función fn1(t,Z)y fn2(t,Z)
% en los archivos m fn1.m y fn2.m
% function W = fn1(t,Z)
% x = Z(1); y = Z(2);
% W = [y, (2*t*y/(1 + t^2) - 2*x/(1 + t^2) + 1)];
% function W = fn2(t,Z)
% x = Z(1); y = Z(2);
% W = [y, (2*t*y/(1 + t^2) - 2*x/(1 + t^2))];
%.......................................................................
% inicia una sección que entra la/s funcion(es) necesarias para el ejemplo
% en archivos m(s) ejecutando el comando diary command en este script.
% el método de programación preferido no usa estas etapas.
% uno debe entrar la función8es) en los M con un editor.
delete output;
delete fn1.m;
diary fn1.m; disp('funcion W = fn1(t,Z)');...
disp('x = Z(1); y = Z(2);');...
disp('W = [y, (2*t*y/(1 + t^2) - 2*x/(1 + t^2) + 1)];');...
diary off;
delete fn2.m;
diary fn2.m; disp('funcion W = fn2(t,Z)');...
disp('x = Z(1); y = Z(2);');...
disp('W = [y, (2*t*y/(1 + t^2) - 2*x/(1 + t^2))];');...
diary off;
fn1(0,[0,0]);
fn2(0,[0,0]);
% . fn1.m fn2.m rks4.m se usan

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
% Usa el método disparo lineal para resolver
% el PVF x'' = 2t/(1+t^2) x' - 2/(1+t^2) + 1
% con x(0) = 1.25 y x(4) = -0.95.
%
% Entrar los extremos de [a,b].
% Entrar valor inicial x(a) en alpha.
% Entrar valor terminal x(b) en beta.
% Entrar número de subintervalos en m.

a = 0;
b = 4;
alpha = 1.25;
beta = -0.95;

Página 104 de 115

 La solución numérica elemental con Octave
m = 20;
% - - - - - - - - - - - - - - - - - - - - - - - -
%
% resolver u'' = p(t)u'(t) + q(t)u(t) + r(t)
% con u(a) = alpha
% y u'(a) = 0

Za = [alpha 0];
[T,Z] = rks4('fn1',a,b,Za,m);
U = Z(:,1)';

% resolver v'' = p(t)v'(t) + q(t)v(t)

% con v(a) = 0
% y v'(a) = 1
Za = [0 1];
[T,Z] = rks4('fn2',a,b,Za,m);
V = Z(:,1)';

% Forma la combinación lineal para x(t).

X = U + (beta - U(m+1))*V/V(m+1);
points = [T;X];

clc; figure(1); clf;

%~~~~~~~~~~~~~~~~~~~~~~~h
% inicio sección gráfica
%~~~~~~~~~~~~~~~~~~~~~~~
a = 0;
b = 4;
c = -1.25;
d = 1.75;
plot([a b],[0 0],'b',[0 0],[c d],'b');
axis([a b c d]);
axis(axis);
hold on;
plot(T,X,'-g');
if m<=30,
plot(T,X,'or');
end;
xlabel('t');
ylabel('x');
Mx1 = 'La solución del disparo lineal a x`` = f(t,x,x`).';
title(Mx1);
grid;
hold off;
figure(gcf);

clc;
%............................................
% inicio sección impresión resultados.
%............................................
Mx2 = ' t(k) x(k)';
clc,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points')

Página 105 de 115

 La solución numérica elemental con Octave

function [T,Z]=rks4(F,a,b,Za, M)

%entrada - F es el sistema ingresado como string 'F'

% - a y b son los extremos del intervalo
% - Za=[x(a) y(a)] condiciones iniciales
% - M número de pasos
%salida - T vector de pasos
% - Z=[x1(t) . . . xn(t)]con xk(t) la aproximación a la
% kth variable dependiente

h=(b-a)/M;
T=zeros(1,M+1);
Z=zeros(M+1,length(Za));
T=a:h:b;
Z(1,:)=Za;

for j=1:M
k1=h*feval(F,T(j),Z(j,:));
k2=h*feval(F,T(j)+h/2,Z(j,:)+k1/2);
k3=h*feval(F,T(j)+h/2,Z(j,:)+k2/2);
k4=h*feval(F,T(j)+h,Z(j,:)+k3);
Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;
end

Ejemplo de la ejecución del método del disparo lineal para resolver

el PVF x'' = 2t/(1+t^2) x' - 2/(1+t^2) + 1
con x(0) = 1.25 y x(4) = -0.95.

octave:4> disparo_lineal
function W = fn1(t,Z)
x = Z(1); y = Z(2);
W = [y, (2*t*y/(1 + t^2) - 2*x/(1 + t^2) + 1)];
function W = fn2(t,Z)
x = Z(1); y = Z(2);
W = [y, (2*t*y/(1 + t^2) - 2*x/(1 + t^2))];
La solución del disparo lineal a x`` = f(t,x,x`).
t(k) x(k)
0.00000 1.25000
0.20000 1.31731
0.40000 1.32643
0.60000 1.28165
0.80000 1.18928
1.00000 1.05673
1.20000 0.89191
1.40000 0.70276
1.60000 0.49699
1.80000 0.28198
2.00000 0.06473
2.20000 -0.14818
2.40000 -0.35052
2.60000 -0.53644
2.80000 -0.70043
3.00000 -0.83726
3.20000 -0.94201
3.40000 -1.01000
3.60000 -1.03678
3.80000 -1.01812
4.00000 -0.95000

Página 106 de 115

 La solución numérica elemental con Octave

% Método del Disparo no lineal TECNICA DE DISPARO PARA EL PROBLEMA NO LINEAL

%
% Aproxima la solución del problema nolineal del valor de frontera
%
% Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
%
%
% ENTRADA: Intervalo A,B; condiciones de contorno ALPHA, BETA; número de
% subintervalos N; tolerancia TOL; máximo número de iteraciones M.
% SALIDA: Approximaciones W(1,I) a Y(X(I)); W(2,I) TO Y'(X(I))
% por cada I=0,1,...,N o un mensaje de que el número de iteraciones
% ha sido excedido.
% EJEMPLO
% Resolver por diferencias finitas la
% EDO con valor de frontera
% y ''= (32+2*x^3-y*y')/8, 1<=x<=3,
% y(1) = 17 , y(3) = 43 / 3 .
% Usando h = 0.5; 0.2; 0.05; 0.025.
% Solución: (La solución exacta es y(x) = x2 +16 / x ).
TRUE = 1;
FALSE = 0;
fprintf(1,'este es el método de disparo no lineal.\n');
fprintf(1,'Ingrese la funcion F(X,Y,Z) en términos de x, y, z.\n');
fprintf(1,'seguida de la deriv. parcial de F respecto a y en \n');
fprintf(1,'la siguiente línea seguida de la de F respecto a \n');
fprintf(1,'z o y-prima en la línea siguiente. \n');
fprintf(1,'ejemplp: (32+2*x^3-y*z)/8 \n'); %z=y' {f(x,y,y'}
fprintf(1,' -z/8 \n'); %fy:=-z/8 {df/dy}
fprintf(1,' -y/8 \n'); %fy1:=-y/8 {df/dy'}
s = input('');
F = inline(vectorize(s),'x','y','z');
s = input('');
FY = inline(vectorize(s),'x','y','z');
s = input('');
FYP = inline(vectorize(s),'x','y','z');
OK=0;
while OK == 0
fprintf(1,'Ingrese extremos derecho e izquierdo en líneas separadas.\n');
A = input(' ');
B = input(' ');
if A >= B
fprintf(1,'extremo izquierdo debe ser menor que el derecho.\n');
else OK = 1;
end;
end;
fprintf(1,'Ingrese Y(%.10e).\n', A);
ALPHA = input(' ');
fprintf(1,'Ingrese Y(%.10e).\n', B);
BETA = input(' ');
TK = (BETA-ALPHA)/(B-A);
fprintf(1,'TK = %.8e\n', TK);
fprintf(1,'Ingresa nuevo TK? Entre Y o N.\n');
AA = input(' ','s');
if AA == 'Y' | AA == 'y'
fprintf(1,'ingrese nuevo TK\n');
TK = input(' ');
end;
OK = 0;
while OK == 0
fprintf(1,'Ingrese un entero > 1 par el número de subintervalos.\n');

Página 107 de 115

 La solución numérica elemental con Octave
N = input(' ');
if N <= 1
fprintf(1,'Número debe ser mayor que 1.\n');
else
OK = 1;
end;
end;
OK = 0;
while OK == 0
fprintf(1,'Ingrese Tolerancia.\n');
TOL = input(' ');
if TOL <= 0
fprintf(1,'Tolerancia debe ser positiva.\n');
else
OK = 1;
end;
end;
OK = 0;
while OK == 0
fprintf(1,'Ingrese número máximo de iteraciones.\n');
NN = input(' ');
if NN <= 0
fprintf(1,'debe ser entero positivo .\n');
else
OK = 1;
end;
end;
if OK == 1
FLAG = 1;
if FLAG == 2
fprintf(1,'Ingrese nombrearchivo en la forma - drive:\\[Link]\n');
fprintf(1,'ejemplo A:\\[Link]\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else
OUP = 1;
end;
fprintf(OUP, 'METODO DISPARO NO LINEAL\n\n');
fprintf(OUP, ' I X(I) W1(I) W2(I)\n');
% STEP 1
W1 = zeros(1,N+1);
W2 = zeros(1,N+1);
H = (B-A)/N;
K = 1;
% TK CALCULADA
OK = 0;
% etapa 2
while K <= NN & OK == 0
% etapa 3
W1(1) = ALPHA;
W2(1) = TK;
U1 = 0 ;
U2 = 1;
% etapa 4
% Método Runge-Kutta para sistemas se usa en etapas 5 y 6
for I = 1 : N
% etapa 5
X = A+(I-1)*H;
T = X+0.5*H;
% etapa 6
K11 = H*W2(I);
K12 = H*F(X,W1(I),W2(I));
Página 108 de 115
 La solución numérica elemental con Octave
K21 = H*(W2(I)+0.5*K12);
K22 = H*F(T,W1(I)+0.5*K11,W2(I)+0.5*K12);
K31 = H*(W2(I)+0.5*K22);
K32 = H*F(T,W1(I)+0.5*K21,W2(I)+0.5*K22);
K41 = H*(W2(I)+K32);
K42 = H*F(X+H,W1(I)+K31,W2(I)+K32);
W1(I+1) = W1(I)+(K11+2*(K21+K31)+K41)/6;
W2(I+1) = W2(I)+(K12+2*(K22+K32)+K42)/6;
K11 = H*U2;
K12 = H*(FY(X,W1(I),W2(I))*U1+FYP(X,W1(I),W2(I))*U2);
K21 = H*(U2+0.5*K12);
K22 = H*(FY(T,W1(I),W2(I))*(U1+0.5*K11)+FYP(T,W1(I),W2(I))*(U2+0.5*K21));
K31 = H*(U2+0.5*K22);
K32 = H*(FY(T,W1(I),W2(I))*(U1+0.5*K21)+FYP(T,W1(I),W2(I))*(U2+0.5*K22));
K41 = H*(U2+K32);
K42 = H*(FY(X+H,W1(I),W2(I))*(U1+K31)+FYP(X+H,W1(I),W2(I))*(U2+K32));
U1 = U1+(K11+2*(K21+K31)+K41)/6;
U2 = U2+(K12+2*(K22+K32)+K42)/6;
end;
% etapa 7
% prueba de exactitud
if abs(W1(N+1)-BETA) < TOL
% etapa 8
I = 0;
fprintf(OUP, '%3d %13.8f %13.8f %13.8f\n', I, A, ALPHA, TK);
for I = 1 : N
J = I+1;
X = A+I*H;
fprintf(OUP, '%3d %13.8f %13.8f %13.8f\n', I, X, W1(J), W2(J));
end;
fprintf(OUP, 'Convergencia en %d iteraciones\n', K);
fprintf(OUP, ' t = %14.7e\n', TK);
% etapa 9
OK = TRUE;
else
% etapa 10
% SE aplica método de Newton para mejorar TK
TK = TK-(W1(N+1)-BETA)/U1;
K = K+1;
end;
end;
% etapa 11
% método falló
if OK == 0
fprintf(OUP, 'Método falló luego de %d iteraciones\n', NN);
end;
end;
if OUP ~= 1
fclose(OUP);
fprintf(1,'archivo de salida %s creado con éxito \n',NAME);
end;

function U=dirich(f1,f2,f3,f4,a,b,h,tol,max1) ECUACIONES ELÍPTICAS

Página 109 de 115

 La solución numérica elemental con Octave

%Entrada - f1,f2,f3,f4 son las funciones de contorno ingresadas como string

% - a y b extremos derechos como [0,a] y [0,b]
% - h tamaño de paso
% - tol es la tolerancia
%SALIDA - U la matríz solución del problema

% EJEMPLO DE CARGA
% f1=@(x)0*x+10
% f2=@(x)0*x+100
% f3=@(x)0*x+50
% f4=@(x)0*x
% U=dirich(f1,f2,f3,f4,4,4,0.5,0.001,20)

%Inicializa parámetros y U

n=fix(a/h)+1;
m=fix(b/h)+1;
ave=(a*(feval(f1,0)+feval(f2,0)) ...
+b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);
U=ave*ones(n,m);

%Condiciones de borde

U(1,1:m)=feval(f3,0:h:(m-1)*h)';
U(n,1:m)=feval(f4,0:h:(m-1)*h)';
U(1:n,1)=feval(f1,0:h:(n-1)*h);
U(1:n,m)=feval(f2,0:h:(n-1)*h);
U(1,1)=(U(1,2)+U(2,1))/2;
U(1,m)=(U(1,m-1)+U(2,m))/2;
U(n,1)=(U(n-1,1)+U(n,2))/2;
U(n,m)=(U(n-1,m)+U(n,m-1))/2;

%parametros SOR

w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));

%Refinando aproximaciones
err=1;
cnt=0;
while((err>tol)&(cnt<=max1))
err=0;
for j=2:m-1
for i=2:n-1
relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4;
U(i,j)=U(i,j)+relx;
if (err<=abs(relx))
err=abs(relx);
end
end
end
cnt=cnt+1;
end

U=flipud(U');

Elíptica
Página 110 de 115
 La solución numérica elemental con Octave
Se plantea la ecuación de Laplace, div(div)=0 en R={(x,y):0≤x≤4,0≤y≤4}con
las CF dadas por:
u(x,0)=10 y u(x,4)=100 para 0<x<4 y
u(0,y)=50 y u(4,y)=0 para 0<y<4
se toma h=0.5 para x e y, o sea 64 cuadrados en la malla
octave:13> f1=@(x)0*x+10;
octave:14> f2=@(x)0*x+100;
octave:15> f3=@(x)0*x+50;
octave:16> f4=@(x)0*x;
octave:17> U=dirich(f1,f2,f3,f4,4,4,0.5,0.001,20)
U =
Columns 1 through 6:
75.00000 100.00000 100.00000 100.00000 100.00000 100.00000
50.00000 72.56216 79.88936 81.64786 80.52445 76.68815
50.00000 60.35926 65.34734 66.17761 63.76179 57.91969
50.00000 53.52760 54.96319 53.95345 50.42541 43.98667
50.00000 48.78785 47.02429 44.24762 39.99980 33.76675
50.00000 44.59979 40.09850 36.01288 31.55942 26.04620
50.00000 39.51325 32.75749 28.14616 24.17901 19.88825
50.00000 30.69626 23.27224 19.63521 17.12241 14.67547
30.00000 10.00000 10.00000 10.00000 10.00000 10.00000
Columns 7 through 9:
100.00000 100.00000 50.00000
68.30843 49.30347 0.00000
47.24210 28.90547 0.00000
33.83481 19.07627 0.00000
25.03420 13.56482 0.00000
18.97039 10.14874 0.00000
14.65245 8.05973 0.00000
11.69142 7.43778 0.00000
10.00000 10.00000 5.00000

function U=crnich(f,c1,c2,a,b,c,n,m) ECUACIÓN PARABÓLICA

%Entrada - f=u(x,0) 'f' como un string

% - c1=u(0,t) y c2=u(a,t)
% - a y b extremos derechos como [0,a] y [0,b]
% - c es la constante en la ecuación de calor

% - n y m numeros de puntos de malla [0,a] y [0,b]

%SALIDA - U la matríz solución del problema

% Ejemplo de carga
% f=@(x)sin(pi*x)+sin(3*pi*x)
% crnich(f,0,0,1,0.5,1,11,11)

%Inicializa parámetros y U

h=a/(n-1);
k=b/(m-1);
r=c^2*k/h^2;
s1=2+2/r;
s2=2/r-2;

Página 111 de 115

 La solución numérica elemental con Octave
U=zeros(n,m);

%Condiciones de borde

U(1,1:m)=c1;
U(n,1:m)=c2;

%Generación de la primera fila

U(2:n-1,1)=feval(f,h:h:(n-2)*h)';

%Se forma la diagonal y los elementos no diagonales de A

%el vector constante B y resuelve el sistema tridiagonal AX=B

Vd(1,1:n)=s1*ones(1,n);
Vd(1)=1;
Vd(n)=1;
Va=-ones(1,n-1);
Va(n-1)=0;
Vc=-ones(1,n-1);
Vc(1)=0;
Vb(1)=c1;
Vb(n)=c2;
for j=2:m
for i=2:n-1
Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1);
end
X=trisys(Va,Vd,Vc,Vb);
U(1:n,j)=X';
end

U=U'

function X=trisys(A,D,C,B)
%Entradas – A es la subdiagonal de la matríz de coeficientes
% - D es la diagonal principal de la matríz de coeficientes
% - C es la diagonal superior de la matríz de coeficientes
% - B es el vector constante del sistema lineal
%SALIDA - X es el vector solución

N=length(B);

for k=2:N
mult=A(k-1)/D(k-1);
D(k)=D(k)-mult*C(k-1);
B(k)=B(k)-mult*B(k-1);
end

X(N)=B(N)/D(N);
for k= N-1:-1:1
X(k)=(B(k)-C(k)*X(k+1))/D(k);
end

Página 112 de 115

 La solución numérica elemental con Octave

Ecuación de calor

(∂u (x,t)/∂t)=c2(∂2u (x,t)/∂x2) en 0<x<1 y 0<t<0.1

u(0,t)=0,u(1,t)=0 son las CF para 0≤t≤0.1
u(x,0)=sin(πx)+sin(3πx) a t=0 y 0≤x≤1
h=0.1,k=0.01,r=1, generando n=11 y m=11
octave:25> f=@(x)1.5-1.5*x;
octave:26> g=@(x)0*x;a=1;b=0.5;c=2;n=10;m=10
m = 10
octave:27> f=@(x)sin(pi*x)+sin(3*pi*x);c1=0;c2=0;c=1;n=11;m=11;
octave:28> forwdif(f,c1,c2,a,b,c,n,m)
ans =
Columns 1 through 6:
0.0000e+000 1.1180e+000 1.5388e+000 1.1180e+000 3.6327e-001 0.0000e+000
0.0000e+000 -2.3681e+000 -2.6692e+000 -5.5174e-001 2.3207e+000 3.6327e+000
0.0000e+000 7.9667e+000 9.4239e+000 3.2231e+000 -5.4817e+000 -9.4871e+000
0.0000e+000 -2.4581e+001 -2.8866e+001 -9.2970e+000 1.8015e+001 3.0567e+001
0.0000e+000 7.6894e+001 9.0409e+001 2.9418e+001 -5.5787e+001 -9.4952e+001
0.0000e+000 -2.4000e+002 -2.8213e+002 -9.1647e+001 1.7441e+002 2.9670e+002
0.0000e+000 7.4934e+002 8.8091e+002 2.8624e+002 -5.4441e+002 -9.2622e+002
0.0000e+000 -2.3395e+003 -2.7503e+003 -8.9362e+002 1.6998e+003 2.8918e+003
0.0000e+000 7.3044e+003 8.5868e+003 2.7900e+003 -5.3070e+003 -9.0287e+003
0.0000e+000 -2.2805e+004 -2.6809e+004 -8.7109e+003 1.6569e+004 2.8189e+004
0.0000e+000 7.1202e+004 8.3703e+004 2.7197e+004 -5.1731e+004 -8.8010e+004
Columns 7 through 11:
3.6327e-001 1.1180e+000 1.5388e+000 1.1180e+000 0.0000e+000
2.3207e+000 -5.5174e-001 -2.6692e+000 -2.3681e+000 0.0000e+000
-5.4817e+000 3.2231e+000 9.4239e+000 7.9667e+000 0.0000e+000
1.8015e+001 -9.2970e+000 -2.8866e+001 -2.4581e+001 0.0000e+000
-5.5787e+001 2.9418e+001 9.0409e+001 7.6894e+001 0.0000e+000
1.7441e+002 -9.1647e+001 -2.8213e+002 -2.4000e+002 0.0000e+000
-5.4441e+002 2.8624e+002 8.8091e+002 7.4934e+002 0.0000e+000
1.6998e+003 -8.9362e+002 -2.7503e+003 -2.3395e+003 0.0000e+000
-5.3070e+003 2.7900e+003 8.5868e+003 7.3044e+003 0.0000e+000
1.6569e+004 -8.7109e+003 -2.6809e+004 -2.2805e+004 0.0000e+000
-5.1731e+004 2.7197e+004 8.3703e+004 7.1202e+004 0.0000e+000

Otro ejemplo:
octave:20> f=@(x)sin(pi*x)+sin(3*pi*x);
octave:21> U=crnich(f,0,0,1,0.5,1,11,11)
U =
Columns 1 through 7:
0.00000 1.11803 1.53884 1.11803 0.36327 0.00000 0.36327
0.00000 -0.09292 0.02699 0.38379 0.78084 0.95342 0.78084
0.00000 0.21099 0.33069 0.33501 0.27955 0.24804 0.27955
0.00000 0.03534 0.09171 0.16788 0.23696 0.26507 0.23696
0.00000 0.05357 0.09342 0.11414 0.12045 0.12113 0.12045
0.00000 0.02137 0.04359 0.06500 0.08118 0.08727 0.08118
0.00000 0.01683 0.03099 0.04092 0.04645 0.04818 0.04645
0.00000 0.00887 0.01723 0.02432 0.02916 0.03089 0.02916
0.00000 0.00585 0.01100 0.01493 0.01736 0.01817 0.01736
0.00000 0.00339 0.00649 0.00900 0.01065 0.01122 0.01065
0.00000 0.00211 0.00400 0.00548 0.00642 0.00674 0.00642
Columns 8 through 11:
1.11803 1.53884 1.11803 0.00000
0.38379 0.02699 -0.09292 0.00000
0.33501 0.33069 0.21099 0.00000
0.16788 0.09171 0.03534 0.00000
0.11414 0.09342 0.05357 0.00000
0.06500 0.04359 0.02137 0.00000
Página 113 de 115
 La solución numérica elemental con Octave
0.04092 0.03099 0.01683 0.00000
0.02432 0.01723 0.00887 0.00000
0.01493 0.01100 0.00585 0.00000
0.00900 0.00649 0.00339 0.00000
0.00548 0.00400 0.00211 0.00000
U =
Columns 1 through 7:
0.00000 1.11803 1.53884 1.11803 0.36327 0.00000 0.36327
0.00000 -0.09292 0.02699 0.38379 0.78084 0.95342 0.78084
0.00000 0.21099 0.33069 0.33501 0.27955 0.24804 0.27955
0.00000 0.03534 0.09171 0.16788 0.23696 0.26507 0.23696
0.00000 0.05357 0.09342 0.11414 0.12045 0.12113 0.12045
0.00000 0.02137 0.04359 0.06500 0.08118 0.08727 0.08118
0.00000 0.01683 0.03099 0.04092 0.04645 0.04818 0.04645
0.00000 0.00887 0.01723 0.02432 0.02916 0.03089 0.02916
0.00000 0.00585 0.01100 0.01493 0.01736 0.01817 0.01736
0.00000 0.00339 0.00649 0.00900 0.01065 0.01122 0.01065
0.00000 0.00211 0.00400 0.00548 0.00642 0.00674 0.00642
Columns 8 through 11:
1.11803 1.53884 1.11803 0.00000
0.38379 0.02699 -0.09292 0.00000
0.33501 0.33069 0.21099 0.00000
0.16788 0.09171 0.03534 0.00000
0.11414 0.09342 0.05357 0.00000
0.06500 0.04359 0.02137 0.00000
0.04092 0.03099 0.01683 0.00000
0.02432 0.01723 0.00887 0.00000
0.01493 0.01100 0.00585 0.00000
0.00900 0.00649 0.00339 0.00000
0.00548 0.00400 0.00211 0.00000

function F=finedif(f,g,a,b,c,n,m) ECUACIÓN HIPERBÓLICA

%Entrada – f*u(x,0) como un string

% g*ut(x,0) como un string
% - a y b extremos derechos como [0,a] y [0,b]
% - c es la constant de la ecuación de onda
% - n y m numeros de puntos de malla [0,a] y [0,b]
%SALIDA - U la matríz solución del problema

%Inicializa parámetros y U
h=a/(n-1);
k=b/(m-1);
r=c+k/h;
r2=r^2;
r22=r^2/2;
s1=1-r^2;
s2=2-2*r^2;
U=zeros(n,m);

%Calcular la primera y segunda fila

for i=2:(n-1)
U(i,1) = feval(f,h*(i-1));
Página 114 de 115
 La solución numérica elemental con Octave
U(i,2) = s1*feval(f,h*(i-1)) + k*feval(g,h*(i-1)) ...
+ r22*(feval(f,h*(i)) + feval(f,h*(i-2)));
end

%Calcular las filas restantes de U

for j=3:m
for i=2:(n-1),
U(i,j) = s2*U(i,j-1) + r2*(U(i-1,j-1) + U(i+1,j-1)) - U(i,j-2);
end
end
U=U';

Ecuación de onda

(∂2u (x,t)/∂t2)=c2(∂2u (x,t)/∂x2) en 0<x<1 y 0<t<0.5

u(0,t)=0,u(1,t)=0 son las CF
u(x,0)= u(x,0)=1.5- 1.5x en 0≤ x≤1
ut(x,0)=0 en 0<x<1
octave:22> f=@(x)1.5-1.5*x;
octave:23> g=@(x)0*x;a=1;b=0.5;c=2;n=10;m=10;
octave:24> finedif(f,g,a,b,c,n,m)
ans =
Columns 1 through 7:
0.00000 1.33333 1.16667 1.00000 0.83333 0.66667 0.50000
0.00000 1.33333 1.16667 1.00000 0.83333 0.66667 0.50000
0.00000 -0.16667 1.16667 1.00000 0.83333 0.66667 0.50000
0.00000 -0.16667 -0.33333 1.00000 0.83333 0.66667 0.50000
0.00000 -0.16667 -0.33333 -0.50000 0.83333 0.66667 0.50000
0.00000 -0.16667 -0.33333 -0.50000 -0.66667 0.66667 0.50000
0.00000 -0.16667 -0.33333 -0.50000 -0.66667 -0.83333 0.50000
0.00000 -0.16667 -0.33333 -0.50000 -0.66667 -0.83333 -1.00000
0.00000 -0.16667 -0.33333 -0.50000 -0.66667 -0.83333 -1.00000
0.00000 -0.16667 -0.33333 -0.50000 -0.66667 -0.83333 -1.00000
Columns 8 through 10:
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
0.33333 0.16667 0.00000
-1.16667 0.16667 0.00000
-1.16667 -1.33333 0.00000

Página 115 de 115