1
Profesor Leopoldo Silva Bijit 20-01-2010
Capítulo 23
Algoritmos numéricos.
23.1. Solución de sistema simultáneo de ecuaciones lineales.
Si se aplica método nodal con modificaciones, para tratar fuentes de voltajes controladas e
independientes, se obtiene un sistema de ecuaciones, del tipo:
A x b
Donde A es la matriz nodal aumentada, x es el vector de incógnitas y b el vector de excitaciones.
Existen dos esquemas generales para resolver sistemas lineales de ecuaciones: Métodos de
eliminación directa y Métodos Iterativos. Los métodos directos, están basados en la técnica de
eliminación de Gauss, que mediante la aplicación sistemática de operaciones sobre los
renglones transforma el problema original de ecuaciones en uno más simple de resolver.
De entre los variados esquemas, basados en la eliminación de Gauss, el método de
descomposición en submatrices triangulares (LU, de Lower y Upper) es preferentemente
empleado en implementaciones computacionales, para sistemas de menos de 300 ecuaciones.
Para sistemas de un mayor número de ecuaciones se emplean métodos iterativos.
La mayoría de estos procedimientos están basados en el método de Gauss Seidel, con
aceleraciones para la convergencia.
23.1.1. Descomposición LU.
Está basado en descomponer la matriz de coeficientes en dos matrices triangulares L y U, según:
A L U
Donde L es una matriz triangular inferior (lower), y U es una matriz triangular superior (upper).
El sistema original de ecuaciones, queda:
L U x b
Que puede ser interpretado como dos sistemas de ecuaciones:
2 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
L d b
U x d
Los dos sistemas anteriores son sencillos de resolver, como se verá más adelante. El sistema con
matriz L, puede ser resuelto por substituciones hacia adelante; el sistema con matriz U se
resuelve por substituciones hacia atrás.
El procedimiento está basado en obtener las matrices L y U, a partir de A; luego en obtener el
vector d; y finalmente en calcular la solución en el vector x.
Existen varias formas de efectuar la descomposición, el método de Doolittle asigna unos a los
elementos de la diagonal principal de L.
Veremos a través de un ejemplo, las principales ideas, intentando obtener un algoritmo para el
cálculo.
Se tiene la matriz A de 4x4 y se desea obtener L y U.
11 12 13 14 11 12 13 14
21 22 23 24 21 22 23 24
31 32 33 34 31 32 33 34
41 42 43 44 41 42 43 44
1 0 0 0
1 0 0 0
1 0 0 0
1 0 0 0
a a a a u u u u
a a a a l u u uA
a a a a l l u u
a a a a l l l u
Efectuando la multiplicación de las matrices L y U, se obtiene:
11 12 13 14
21 11 21 12 22 21 13 23 21 14 24
31 11 31 12 32 22 31 13 32 23 33 31 14 32 24 34
41 11 41 12 42 22 41 13 42 23 43 33 41 14 42 24 43 34 44
u u u u
l u l u u l u u l u uA
l u l u l u l u l u u l u l u u
l u l u l u l u l u l u l u l u l u u
El primer renglón de A permite, por comparación, determinar el primer renglón de U.
11 11 12 12 13 13 14 14 u ; u ; u ; ua a a a
Una vez conocido u11, la primera columna de A permite determinar el primer renglón de L, se
obtienen:
21 21 11 31 31 11 41 41 11/ ; / ; /l a u l a u l a u
Algoritmos numéricos 3
Profesor Leopoldo Silva Bijit 20-01-2010
El segundo renglón de A, permite calcular el segundo renglón de U, una vez conocidos los
elementos del primer renglón de U, se tienen:
21 12 22 22 21 13 23 23 21 14 24 24; ; l u u a l u u a l u u a
Despejando los elementos del segundo renglón de U, se obtienen:
22 22 21 12
23 23 21 13
24 24 21 14
u a l u
u a l u
u a l u
La segunda columna de A, permite calcular la segunda columna de L.
31 12 32 22 32 41 12 42 22 42; l u l u a l u l u a
Despejando los elementos de la segunda columna de L. se obtienen:
32 32 31 12 22
42 42 41 12 22
( ) /
( ) /
l a l u u
l a l u u
Del tercer renglón de A, resultan:
31 13 32 23 33 33 31 14 32 24 34 34; l u l u u a l u l u u a
Las que permiten despejar los elementos del tercer renglón de U:
33 33 31 13 32 23
34 34 31 14 32 24
u a l u l u
u a l u l u
De la tercera columna de A, se puede calcular la tercera columna de L:
43 43 41 13 42 23 33( ) /l a l u l u u
Finalmente, el cuarto renglón de A, permite calcular el cuarto renglón de U.
44 44 41 14 42 24 43 34u a l u l u l u
Si bien se ha desarrollado para una matriz de 4x4, de las expresiones obtenidas puede inducirse
relaciones generales, como veremos a continuación:
Con N es el número de renglones y columnas de A.
4 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Para: 1,..., ;n N
, 1n nl
La diagonal de L, se obtiene por definición de la descomposición de Doolittle.
El n-avo renglón de U se obtiene según: 1
, , , ,1
n
n i n i n k k ik
u a l u
Para: ,..., ;i n N
Y la n-ava columna de L con:
1
, , , , ,1
/n
j n j n j k k n n nk
l a l u u
Para: 1,...,j n N
A continuación se obtiene el algoritmo para la substitución hacia delante:
De la relación:
L d b
Se obtiene:
11 1 1
21 22 2 2
31 32 33 3 3
41 42 43 44 4 4
0 0 0 d b
0 0 d b
0 d b
d b
l
l l
l l l
l l l l
Efectuando las multiplicaciones, en el lado derecho, se tienen:
11 1 1
21 1 22 2 2
31 1 32 2 33 3 3
41 1 42 2 43 3 44 4 4
b
b
b
b
l d
l d l d
l d l d l d
l d l d l d l d
Las componentes del vector d, se obtienen según:
Algoritmos numéricos 5
Profesor Leopoldo Silva Bijit 20-01-2010
1 1 11
2 2 21 1 22
3 3 31 1 32 2 33
4 4 41 1 42 2 43 3 44
/
( ) /
( ) /
( ) /
d b l
d b l d l
d b l d l d l
d b l d l d l d l
Una vez obtenido d1, se substituye en la expresión siguiente para calcular d2; con d1 y d2, se
puede calcular d3; y así sucesivamente. Por esta razón, al procedimiento se lo denomina
substitución hacia adelante (forward).
El vector d, puede recalcularse para diferentes valores del vector b, que es la situación que se
produce en un barrido DC. Debido a que en el método de Gauss se ocupa, desde el inicio de las
operaciones, los valores de b; el efectuar cálculos con b variable lo realiza con ventajas el
método de descomposición triangular.
La relación anterior, permite deducir una expresión para calcular los di, en una matriz de orden
N.
1
( ) /i l
i i ij j ii
j
d b l d l
Para: 1,2, ,i N
En la descomposición de Doolittle, los iil son unos.
El algoritmo para la substitución hacia atrás se obtiene de manera similar a las anteriores.
Para la triangular superior:
U x d
Se tiene:
11 12 13 14 1 1
22 23 24 2 2
33 34 3 3
44 4 4
x d
0 x d
0 0 x d
0 0 0 x d
u u u u
u u u
u u
u
Efectuando las multiplicaciones, se obtiene:
11 1 12 2 13 3 14 4 1
22 2 23 3 24 4 2
33 3 34 4 3
44 4 4
d
d
d
d
u x u x u x u x
u x u x u x
u x u x
u x
Despejando los xi, se obtienen:
6 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
4 4 44
3 3 34 4 33
2 2 23 3 24 4 22
1 1 12 2 13 3 14 4 11
/
( ) /
( ) /
( ) /
x d u
x d u x u
x d u x u x u
x d u x u x u x u
Que entrega la solución del sistema de ecuaciones. Nótese que primero se obtiene x4; y luego x3,
que se calcula en términos de x4; y así sucesivamente. Por esta razón a este algoritmo se lo
denomina substitución hacia atrás (back).
En general:
1
/N N NN
N
i ij j
j i
i
ii
x d u
d u x
xu
Para: ( 1), ( 2), ,3,2,1i N N
La observación cuidadosa de la generación de las matrices L y U, muestra que no es necesario
emplear espacio adicional para éstas. Los valores que se van calculando pueden almacenarse en
la matriz A. Tampoco es necesario almacenar los elementos de la diagonal principal de L, ya que
son unos por definición. Entonces luego de la descomposición, la matriz original queda, en el
caso del ejemplo de 4x4:
11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44
u u u u
l u u uA
l l u u
l l l u
void ludcmp(float **a, int N)
/*Dada a[1..N][1..N], la reemplaza por la descomposición LU Doolittle.*/
{ int n,i,j,k;
float sum;
for (n=1; n<=N; n++)
{
for (i=n; i<=N; i++)
{ sum=0;
for (k=1; k<=(n-1); k++) sum += a[n][k]*a[k][i];
a[n][i]-=sum;
}
1
, , , ,1
n
n i n i n k k ik
u a l u
Algoritmos numéricos 7
Profesor Leopoldo Silva Bijit 20-01-2010
for (j=n+1; j<=N; j++)
{ sum=0;
for (k=1; k<=(n-1); k++) sum += a[j][k]*a[k][n];
a[j][n]=(a[j][n]-sum)/a[n][n];
}
}
}
/*Realiza substituciones hacia adelante y hacia atrás.
a[1..N][1..N] es de entrada y debe contener la descomposición L y U.
b[1..N] el vector de excitación.
Retorna en b la solución.
a no es modificada y pueden realizarse sucesivos llamados
con diferentes valores de b.
*/
void lufwbksb(float **a, int N, float b[])
{ int i,j;
float sum;
for (i=1; i<=N; i++)
{ sum=0;
for (j=1; j<=(i-1); j++) sum += a[i][j]*b[j];
b[i]-=sum; //l[i][i]=1.
} //se genera vector d en el espacio ocupado por b.
b[N]/=a[N][N];
for (i=(N-1); i>=1; i--)
{ sum=0;
for (j=(i+1); j<=N; j++) sum += a[i][j]*b[j];
b[i]=(b[i]-sum)/a[i][i];
} //se almacena solución x en el espacio ocupado por b.
}
Ejemplo de uso.
//Resuelve el sistema de ecuaciones lineales a·X = b.
ludcmp(a, n);
lufwbksb(a, n, b);
23.1.2. Métodos iterativos.
Para deducir expresiones generales que permitan escribir algoritmos iterativos, consideremos el
sistema lineal de tres ecuaciones:
11 12 13 1 1
21 22 23 2 2
31 32 33 3 3
a a a x b
a a a x b
a a a x b
1
, , , , ,1
/n
j n j n j k k n n nk
l a l u u
1
( ) /i l
i i ij j ii
j
d b l d l
1
/N N NN
N
i ij j
j i
i
ii
x d u
d u x
xu
8 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Despejando de la primera ecuación, la variable 1x ; de la segunda
2x ; y de la tercera 3x ,
obtenemos:
1 1 12 2 13 3 11
2 2 21 1 23 3 22
3 3 31 1 32 2 33
( ) /
( ) /
( ) /
x b a x a x a
x b a x a x a
x b a x a x a
Si consideramos conocidos los valores de las variables del lado derecho, podremos estimar un
nuevo valor para las variables del lado izquierdo de las ecuaciones. Podemos anotar lo anterior,
mediante:
1 1 12 2 13 3 11
2 2 21 1 23 3 22
3 3 31 1 32 2 33
[ 1] ( [ ] [ ]) /
[ 1] ( [ ] [ ]) /
[ 1] ( [ ] [ ]) /
x n b a x n a x n a
x n b a x n a x n a
x n b a x n a x n a
Donde [ 1]ix n es el valor de ix en la iteración (i+1); y [ ]ix n es el valor obtenido en la
iteración anterior.
Durante el proceso iterativo se verifica la convergencia calculando el mayor cambio relativo
entre una iteración y la siguiente, y comparando el valor absoluto de esta diferencia con la
tolerancia deseada.
| [ 1] [ ] | i ix n x n tolerancia
Si el error es menor que la exactitud requerida el proceso termina; en caso contrario se realiza
una nueva iteración.
Si se tienen N variables, pueden generalizarse las iteraciones según:
1
1 1
[ 1] ( [ ] [ ]) /j i j N
i i ij j ij j ii
j j i
x n b a x n a x n a
El esquema anterior se reconoce como método de Jacobi.
Si el cálculo de las variables se realiza en orden, desde 1x hasta Nx , puede observarse que una
vez obtenido 1x puede usarse este valor para calcular 2x ; y así sucesivamente. Entonces en el
cálculo ix se pueden emplear los nuevos valores de las variables desde 1x hasta 1ix .
Entonces el esquema iterativo puede plantearse:
Algoritmos numéricos 9
Profesor Leopoldo Silva Bijit 20-01-2010
1
1 1
[ 1] ( [ 1] [ ]) /j i j N
i i ij j ij j ii
j j i
x n b a x n a x n a
El que se denomina método de Gauss Seidel.
Mejores resultados se logran calculando las variables en orden decreciente de los valores de la
diagonal principal.
Una mejora notable de la convergencia se logra empleando un promedio ponderado de los
resultados de las dos últimas iteraciones para obtener el nuevo valor. Esto se denomina método
de sucesivas sobre relajaciones (SOR Successive Over-Relaxation).
[ 1] [ 1] (1 ) [ ]i i ix n ax n a x n
Con: 0 2a
Si a es 1, se tiene la fórmula de Gauss Seidel. Con a>1, el nuevo valor, en la iteración (n+1),
tiene mayor importancia. Con a<1, se tiene subrelajación. La elección de este valor, y su
influencia en la convergencia debería aclararse en un curso de análisis numérico.
La recurrencia para encontrar el nuevo valor por el método SOR:
1
1 1
[ 1] (1 ) [ ] ( [ 1] [ ]) /j i j N
i i i ij j ij j ii
j j i
x n a x n a b a x n a x n a
El algoritmo descrito en pseudo código:
SOR. Dados: N, a, b, x[0], tol, nmax.
for n = 1, … nmax do Comienza iteraciones
for i = 1,.... N do
1
1 1
[ ] (1 ) [ ] ( [ ] [ ]) /j i j N
i i i ij j ij j ii
j j i
y n a x n a b a y n a x n a
end for i
i iy x
for j = 1,….. N do
j jx y
end for i
if ( tol ) stop
end for n
10 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
/*Dados a[1..N][1..N], b[1..N] calcula x[1..N]
Con nmax iteraciones y factor de relajación alfa. Retorna solución en x[1..N]*/
int sor(float **a, int N, float *x, float *b, int nmax, float alfa)
{ int n, i, j, and;
float sumi, sums;
float *y;
y=vector(1, N); //vector de flotantes
for (n=1; n<=nmax; n++)
{
for(i=1; <=N; i++)
{
sumi=0; for (j=1; j<=(i-1); j++) sumi += a[i][j]*y[j];
sums=0; for (j=i+1; j<=N; j++) sums += a[i][j]*x[j];
y[i]=(1.-alfa)*x[i] + alfa*(b[i]-sumi-sums)/a[i][i];
}
and=1;
for(j=1; j<=N; j++)
{ and = and && (fabs(y[j]-x[j]) < tol*(fabs(x[j])); //error relativo
if( and==0) break;
}
for(j=1; j<=N; j++) x[j]=y[j];
if (and) break;
//putchar('.');
}
free_vector(y, 1, N);
return (n);
}
En lugar de emplear errores absolutos: (fabs(y[j]-x[j]) < tol), es preferible utilizar errores
relativos: (fabs(y[j]-x[j]) < tol*(fabs(x[j]) );
23.2. Solución numérica de sistemas de ecuaciones diferenciales.
Una ecuación diferencial de primer orden puede resolverse numéricamente mediante
integración.
Si se tiene:
( )( )
dr tF t
dt
Entonces:
0
( ) (0) ( )
t
r t r F d
( )F t considera la variación de r(t) y de las excitaciones que producen la respuesta r(t).
Algoritmos numéricos 11
Profesor Leopoldo Silva Bijit 20-01-2010
Una manera simple y aproximada de realizar la integración es calcular el área mediante la suma
de rectángulos, que estudiaremos como el método de Euler.
Una mejor aproximación se logra sumando trapecios, si se desea mayor precisión se emplea
aproximación por segmentos parabólicos, con la regla de Simpson. Para disminuir la
acumulación de errores se emplea el método de Runge-Kutta.
23.2.1. Formulación de ecuaciones de estado.
La formulación de las ecuaciones de una red eléctrica en términos de las variables de estado
permite encontrar la solución de un sistema de ecuaciones diferenciales de primer orden en el
dominio del tiempo. La solución numérica, que veremos a continuación, puede extenderse a
sistemas no lineales.
La representación se logra con un sistema de ecuaciones diferenciales de primer orden:
dxAx Bu
dt
Donde x es el vector de estado, u es el vector de entrada o de excitaciones.
El resto de las variables del sistema puede expresarse en términos del estado, según:
y Cx Du
Donde y es el vector de salida.
A se denomina matriz de estado del sistema, B es la matriz de entrada, C es la matriz de salida, y
D se denomina matriz de alimentaciones directas (feedforward).
23.2.2. Método de Euler.
A partir de la expansión en serie de Taylor, para una variable escalar x, se tiene:
2
2
2
( ) 1 ( )( ) ( ) ....
2
dx t dx tx t t x t t t
dt dt
La relación anterior, puede generalizarse considerando a x como el vector de estado. Pueden
calcularse, aproximadamente, los valores de las variables de estado en el instante siguiente
(k+1), a partir de los valores en el instante k-ésimo, mediante:
( )[ 1] [ ] i k
i i
dx tx k x k t
dt
Este procedimiento iterativo se denomina esquema simple de Euler.
12 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Los valores de las derivadas, en un instante determinado, se obtienen mediante la matriz de
estado.
1 1
[ 1] [ ] ( [ ] [ ])j n j n
i i ij j ij j
j j
x k x k a x k b u k t
Para: 1,2, ,i n
A partir de la ecuación de estado se determina el valor de las derivadas en un punto.
La siguiente función calcula en la matriz x, los valores de las variables de estado en npuntos
separados en intervalos de tiempo Delta.
No se considera la matriz b, ni el vector u de excitaciones. Esto equivale a resolver un sistema
de ecuaciones diferenciales lineales homogéneas y de primer orden.
void euler(float **a, int N, float **x, float *ic, int npuntos, float Delta)
/*Dados a[1..N][1..N], ic[1..N] calcula x[1..N][1..npuntos]*/
{ int i, j, k;
float sum, t=0.;
for(i=1; i<=N; i++) x[i][1]=ic[i]; //condiciones iniciales.
for (k=1; k<npuntos; k++)
{ t= t+Delta;
for(i=1; i<=N; i++)
{
sum=0; for (j=1; j<=N; j++) sum += a[i][j]*x[j][k];
x[i][k+1]= x[i][k]+sum*Delta;
}
}
}
Una alternativa de diseño es generar un archivo de datos en lugar de almacenar los puntos en
una matriz. Con el archivo de datos se pueden generar formas de ondas.
La siguiente función ilustra la generación de un archivo de datos compatible con el comando
pointplot de Maple. Generando los puntos (t[k], x[i][k]) para la variable ix .
void genseq(float **a, int N, float **x, int npuntos, float Delta, int i)
/*Dados a[1..N][1..N], ic[1..N] y x[1..N][1..npuntos]
genera seq compatible para gráficos de tipo pointplot en Maple.*/
{ int k;
float t=0.;
printf("Seq:=[");
for (k=1; k<=npuntos; k++, t=t+Delta)
{ printf("[%g,%g]\n", t, x[i][k]);
if (k<npuntos) putchar(',');
}
putchar(']'); putchar('\n');
Algoritmos numéricos 13
Profesor Leopoldo Silva Bijit 20-01-2010
}
Ejemplo de uso:
n=2;npuntos=20;
a=matrix(1, n, 1, n); //pide espacio
a[1][1]=0.; a[1][2]=1.;
a[2][1]=-3.; a[2][2]=-2.;
ic=vector(1, n);
ic[1]=1.;ic[2]=0.;
x=matrix(1, n, 1, npuntos); //pide espacio
//Resuelve sistema ecuaciones diferenciales.
euler(a,n,x,ic,npuntos,0.1);
genseq(a,n,x,npuntos,0.1,1);
La última invocación genera para la variable 1x :
Seq:=[[0,1],[0.1,1],[0.2,0.97],[0.3,0.916],[0.4,0.8437],[0.5,0.75838],[0.6,0.664813]
,[0.7,0.567208],[0.8,0.46918],[0.9,0.373741],[1,0.283314],[1.1,0.199761]
,[1.2,0.124418],[1.3,0.0581519],[1.4,0.00140604],[1.5,-0.0457352]
,[1.6,-0.0834903],[1.7,-0.112322],[1.8,-0.132883],[1.9,-0.145962]
]
Debido a la acumulación de errores no suele emplearse el algoritmo de Euler.
23.2.3. Algoritmo de Euler.
Para el caso de una variable, tenemos:
( )( , ( ))
dy tf t y t
dt
Se pueden obtener los valores sucesivos de y, mediante:
1 ( , ( ))n n n ny y hf t y t
Si definimos:
1 ( , )n nk hf t y
Se realiza la integración mediante:
1 1n ny y k
Se tiene que ( , )n nf t y es la derivada de la función ( )y t , y 1k es el área del rectángulo bajo la
curva de ( , )f t y , y también el incremento de la ordenada. Las relaciones se muestran en la
Figura 1, para la variable ( )y t .
14 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
1 1
1
( ) ( , )n nn n n
n n
y y ktg f t y
t t h
tn tn+1
k1 yn
yn+1
h tn tn+1 h
f(tn,yn) n
f(tn+1,yn+1)
Figura 23.1. Método de Euler.
23.2.4. Algoritmo trapezoidal.
Una mejor aproximación para el área bajo la curva de ( , )f t y es mediante el trapecio entre las
paralelas nt y 1nt .
tn tn+1
(k1 + k2)/2
yn
yn+1
h tn tn+1 h
f(tn,yn) n
f(tn+1,yn+1)
Figura 23.2. Método trapezoidal.
Entonces:
1 1 1( ( , ) ( , ))2
n n n n n n
hy y f t y f t y
Con:
1
2 1 1
( , )
( , )
n n
n n
k hf t y
k hf t y
Se realiza la integración trapezoidal mediante:
Algoritmos numéricos 15
Profesor Leopoldo Silva Bijit 20-01-2010
1 1 2
1( )
2n ny y k k
En cada paso de integración es preciso evaluar dos veces la función: ( , )f t y . La determinación
de 1ny que se emplea en el cálculo de 2k , puede determinarse, aplicando Euler,
según: 1 1n ny y k
23.2.5. Algoritmo de Simpson.
Para refinar el cálculo del área, se define un punto dentro del intervalo, de este modo puede
calcularse el área bajo una parábola que pasa por los tres puntos:
0 0( , )f t y , 1 1( , )f t y , 2 2( , )f t y , con 1 0 / 2t t h y 2 0t t h .
tn tn+1
(k0 +4 k1+ k2)/6
yn
yn+1
h tn tn+1 h
f(tn,yn) n
f(tn+1,yn+1)
fa(t)
f(tn+h/2, y(tn+h/2))
Figura 23.3. Método de Simpson.
Si la ecuación de la parábola que aproxima el área bajo la curva de ( , )f t y es:
2( )af t at bt c
Se tienen: 2
0 0 0
2
1 1 1
2
2 2 2
f at bt c
f at bt c
f at bt c
Que permiten calcular , ,a b c , conociendo: 0 1 2 0 1 2, , , , ,f f f t t t . Luego se realiza la integral:
2
0
( )
t t
a
t t
A f t dt
16 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Aplicando que 1 0 / 2t t h y
2 0t t h , se obtiene:
0 1 2( 4 )6
hA f f f
Con:
0
1
2
( , )
( / 2, ( / 2))
( , ( ))
n n
n n
n n
k hf t y
k hf t h y t h
k hf t h y t h
El paso de integración puede realizarse según:
1 0 1 2
1( 4 )
6n ny y k k k
En cada paso de integración es preciso evaluar tres veces la función: ( , )f t y .
Los valores de las ordenadas intermedias, pueden calcularse, según:
0
0
( / 2) / 2
( )
n n
n n
y t h y k h
y t h y k h
La fórmula anterior permite deducir la conocida fórmula para la aproximación de Simpson.
Si se tienen cinco valores de tiempos para los cuales se calcula la integral, se tiene el área
acumulada, según:
0 1 2 2 3 4
1 1( 4 ) ( 4 )
6 6Area k k k k k k
La cual puede simplificarse a:
0 1 2 3 4
1( 4 2 4 )
6Area k k k k k
Con 2n+1 puntos en total, se tiene en general:
0 1 2 3 2 2 2 1 2( 4 2 4 ...... 2 4 )6
n n n
hArea f f f f f f f
double SimpsonIntegral(double a, double b, int n)
{ double h, suma, sumap, sumai;
int i;
if (n%2==1) n++; //deja a n como número par.
h=(b-a)/n;
suma=f(a)+f(b); //suma los extremos
for(i=1, sumap=0; i<n; i+=2){
sumap+=f(a+i*h); //acumula impares
Algoritmos numéricos 17
Profesor Leopoldo Silva Bijit 20-01-2010
}
for(int i=2, sumai=0; i<n; i+=2){
sumai+=f(a+i*h); //acumula pares.
}
return ((suma+4*sumap+2*sumai)*h/3);
}
23.4.6. Métodos multietapas.
Existen métodos más elaborados para efectuar un paso de integración, En éstos la función se
evalúa en varias etapas entre dos puntos. Los puntos de las etapas sólo se emplean para el
cálculo del nuevo valor.
Consisten en definir una función en la cual se escogen los parámetros de tal modo de minimizar
los errores.
Sea la función , ;n ny t h , entonces, la integración se realiza mediante:
1 , ;n n n ny y h y t h
El método de Euler es: , ; ( , )n n n ny t h f y t
El siguiente esquema es de Runge-Kutta, de segundo orden, con cuatro parámetros.
1
2 1
1 1 2
,
,
n n
n n
n n
k f y t
k f y hk t h
y y h ak bk
Con:
, ( ), ( ) ( ), ,n n n n n n n ny t t af y t t bf y t hf y t t t h
Donde deben determinarse: , , ,a b por consideraciones de exactitud.
Si se define el error de la aproximación por:
,n n n n nT y t h y t h y t t
Para el primer término de Tn, por expansión de Taylor de y(t), resulta:
2 32 3
4
2 3
( ) ( ) ( )( )
2! 3!
n n nn n
dy t d y t d y th hy t h y t h O h
dt dt dt
Como se tiene:
18 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
,dy
f y t tdt
La segunda derivada de y, se obtiene según:
2
2, t y
d y d f f dy f ff y t t f f f f
dt dt t y dt t y
Las derivadas parciales de f se han representado, en forma abreviada, mediante subíndices.
Para la tercera derivada, se tiene:
3
3( ) ( )
d y d f f d f d f f dff f
dt dt t y dt t dt y y dt
La cual permite obtener:
3 2 2 2 2
3 2 2
d y f f f f f f ff f f f
dt t t y y t y y t y
Finalmente para la tercera derivada:
3
3( ) ( )tt yt yt yy y t y
d yf f f f f f f f f f f
dt
Reemplazando las derivadas obtenidas, en la expansión de Taylor, se logra:
2 3
2 2 4( ) ( 2 ) ( )2! 3!
n n t y tt ty yy y t y
h hy t h y t hf f f f f f f f f f f f f O h
Donde se ha empleado: ,n nf f y t t
Para el segundo término de Tn, se requiere expandir el término que está multiplicado por b. Para
esto se emplea la expansión de Taylor de dos variables de ,n nf y y t t ,
con: , y hf t h . Como está multiplicado por h, sólo es necesario considerar hasta
los términos de segundo orden.
Entonces;
Algoritmos numéricos 19
Profesor Leopoldo Silva Bijit 20-01-2010
2 2 22 2
2 2
, ,
( , ) ( , ) 1 ( , ) ( , ) 1 ( , )
2 2
f y y t t f y t
f y t f y t f y t f y t f y ty t y y t t
y t y y t t
Reemplazando en la anterior, los valores de los incrementos y empleando notación abreviada
para las derivadas parciales de f, se logra:
2 2 3
, ,
1 1( ) ( ) ( ) ( )( ) ( ) ( )
2 2
n n n n
y t yy yt tt
f y y t t f y t
f hf f h f hf f hf h f h O h
Reemplazando la expresión anterior en:
, ( ), ( ) ( ), ,n n n n n n n ny t t af y t t bf y t hf y t t t h
Se obtiene:
2 2 2 31 1, ( ( ) ( ) )
2 2n n y t yy yt tty t t af b f f hf f h f hf f hf h f h O h
Donde las funciones y las derivadas parciales están evaluadas en tn.
Se tiene finalmente para el error:
2
2
22
2 2 2 4
22 3!
( ) ( )2 2
n t y tt ty yy y t y
yy tty t yt
h hT h f f ff f f f f f f f f f
f f fh af b f f f f h h f fh h O h
Los parámetros: , , ,a b , se escogen del tal modo de minimizar el error. En este caso pueden
hacerse cero los coeficientes de h y 2h , pero no es posible eliminar la parte que depende de
3h .
Para eliminar la parte proporcional a h , se debe cumplir: 0f af bf
Para eliminar la parte proporcional a2h , se debe cumplir: ( ) / 2 0t y y tf ff bf f bf
De la primera se obtiene: 1a b
De la segunda, debe cumplirse: 1 1
( ) ( ) 02 2
t yf b ff b , de la que se desprenden:
1
2b b
20 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Con 1, 0a b se tiene el método de Euler. No pueden ajustarse , . Por esta razón puede
decirse que el algoritmo de Euler pertenece a la familia Runge-Kutta de segundo orden.
23.4.6.1. Método de Heun.
Considera: 1 , 0 , 1
2a b que satisfacen las dos relaciones anteriores.
Se realiza la integración mediante:
1 1 22
n n
hy y k k
Con:
1 ,n nk f y t
2 1 ,n nk f y k h t h
23.4.6.2. Método del punto medio.
Considera: 1
2 , 0, 1a b que satisfacen las dos relaciones anteriores.
1 2n ny y hk
Con:
1 ,n nk f y t
12 ,
2 2n n
k h hk f y t
La pendiente de y está dada ahora por el valor de f en el punto medio del intervalo.
tn tn+1
hk2 yn
yn+1
h tn tn+1 h
f(tn,yn) n
f(tn+1,yn+1)
f(tn+h/2,yn+k1h/2)
Figura 23.4. Método del punto medio.
23.4.6.3. Método de Ralston.
Considera: 3
4
1 2,
3 3a b que satisfacen las dos relaciones anteriores.
Algoritmos numéricos 21
Profesor Leopoldo Silva Bijit 20-01-2010
Se integra mediante:
1 1 2( 2 )3
n n
hy y k k
Con:
1 ,n nk f y t
12
3 3,
4 4n n
k h hk f y t
En los algoritmos del método de Runge-Kutta de segundo orden, debe evaluarse dos veces la
función: ( ),f y t t , para obtener el siguiente punto de la función.
23.4.7. Métodos de Runge-Kutta de cuarto orden.
Para derivar el algoritmo de cuarto orden de Runge Kutta, se definen los 10 parámetros:
1 1 2 2 3 3 1 2 3 4, , , , , , , , ,a b a b a b w w w w . Los cuales deben elegirse de tal modo de eliminar los errores
proporcionales hasta la cuarta potencia del intervalo temporal h.
1
2 1 1 1
3 2 2 2
4 3 3 3
1 1 1 2 2 3 3 4 4
,
,
,
,
n n
n n
n n
n n
n n
k hf y t
k hf y b k t a h
k hf y b k t a h
k hf y b k t a h
y y w k w k w k w k
Con:
1 2 1 1 1
3 2 2 2 4 3 3 3
, , ,
, ,
n n n n n n
n n n n
y t t w f y t w f y b k t a h
w f y b k t a h w f y b k t a h
La expresión para el error es:
1 ,n n n n nT y y h y t t
Antes se obtuvieron las tres primeras derivadas de y. Se requiere calcular la cuarta derivada de
y, se tiene: 4
2 2
4( 2 )tt ty yy y t y
d y df f f f f f f f f
dt dt
Derivando, se obtiene:
22 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
4 2
4 2
2 3 2
2 3 2
(( ) )
2( ) ( )
yyy yyt yy yyt ytt
yy yt y yyt ytt yt ytt ttt
d y dy dy d y dy dyf f f f f
dt dt dt dt dt dt
dy d y d y dy dy d y dyf f f f f f f f
dt dt dt dt dt dt dt
Reemplazando la primera derivada de y:
4 2
4 2
2 3 2
2 3 2
(( ) )
2( ) ( )
yyy yyt yy yyt ytt
yy yt y yyt ytt yt ytt ttt
d y d yf f f f f f f f f
dt dt
d y d y d yf f f f f f f f f f f f
dt dt dt
Reemplazando la segunda derivada de y, se tiene:
4
4
3
3
(( ) ( ) )
2( )( ) ( ) ( )
yyy yyt yy t y yyt ytt
yy yt t y y yyt ytt yt t y ytt ttt
d yf f f f f f f f f f f f
dt
d yf f f f f f f f f f f f f f f f f f
dt
Reemplazando la tercera derivada de y, se obtiene:
4
4(( ) ( ) )
2( )( ) ( ( )
( )) ( ) ( )
yyy yyt yy t y yyt ytt
yy yt t y y tt ty yt yy
y t y yyt ytt yt t y ytt ttt
d yf f f f f f f f f f f f
dt
f f f f f f f f f f f f f f
f f f f f f f f f f f f f f f
Reemplazando las derivadas obtenidas, en la expansión de Taylor, se logra la serie: 2 3
2 2
43 2 2
3 2 5
( ) ( 2 )2! 3!
( 3 4 3 3 54!
3 ) ( )
n n t y tt ty yy y t y
yyy yyt yy y yy t yyt yt y
yt t y y t y tt ttt
h hy t h y t hf f f f f f f f f f f f f
hf f f f f f f f f f f f f f f
f f f f f f f f f O h
Donde se ha empleado: ,n nf f y t t
La segunda componente del error, requiere calcular:
Algoritmos numéricos 23
Profesor Leopoldo Silva Bijit 20-01-2010
1 1 2 2 3 3 4 4,n nh y t t w k w k w k w k
Se emplea la expansión de Taylor de dos variables de ,n nf y y t t , con los
diferentes , y t . Como está multiplicado por h, sólo es necesario considerar hasta los
términos de tercer orden.
Entonces;
2 2
3 2 2 3
, ,
1 1
2 2
1 1 1 1
6 2 2 6
n n n n
y t yy yt tt
yyy yyt ytt ttt
f y y t t f y t
f y f t f y f y t f t
f y f y t f y t f t
Reemplazando en la anterior, los valores de los incrementos y empleando notación abreviada
para las derivadas parciales de f, se logra:
1k hf
k2 =
1
6b1 3 f3 fyyy
1
2b1 2 f2 fyyt a1
1
2b1 f a12 fytt
1
6a13 fttt h4
1
2fyy b1 2 f2 b1 f fyt a1
1
2ftt a12 h3 ( )fy b1 f ft a1 h2 h f
k3 =
fy b2fyy b1 2 f2
2b1 f fyt a1
ftt a12
2fyy b2 2 f ( )fy b1 f ft a1
b2 ( )fy b1 f ft a1 fyt a2b2 2 f2 fyyt a2
2
b2 f a22 fytt
2
b2 3 f3 fyyy
6
a23 fttt
6
h4 fy b2 ( )fy b1 f ft a1ftt a22
2
fyy b2 2 f2
2b2 f fyt a2 h3
( )fy b2 f ft a2 h2 h f
k4 =
h fa33 fttt
6fyy b3 2 f ( )fy b2 f ft a2
b3 3 f3 fyyy
6
b3 f a32 fytt
2
fy b3 fy b2 ( )fy b1 f ft a1ftt a22
2
fyy b2 2 f2
2b2 f fyt a2
24 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
b3 ( )fy b2 f ft a2 fyt a3b3 2 f2 fyyt a3
2h4
fyy b3 2 f2
2
ftt a32
2fy b3 ( )fy b2 f ft a2 b3 f fyt a3 h3
( )fy b3 f ft a3 h2
Reemplazando las expresiones anteriores en: ,n ny t t
Se obtiene:
,n nh y t t
w2b1 3 f3 fyyy
6
b1 2 f2 fyyt a1
2
b1 f a12 fytt
2
a13 fttt
6w4
a33 fttt
6
fyy b3 2 f ( )fy b2 f ft a2b3 3 f3 fyyy
6
b3 f a32 fytt
2
fy b3 fy b2 ( )fy b1 f ft a1ftt a22
2
fyy b2 2 f2
2b2 f fyt a2
b3 ( )fy b2 f ft a2 fyt a3b3 2 f2 fyyt a3
2w3
fy b2fyy b1 2 f2
2b1 f fyt a1
ftt a12
2fyy b2 2 f ( )fy b1 f ft a1
b2 ( )fy b1 f ft a1 fyt a2b2 2 f2 fyyt a2
2
b2 f a22 fytt
2
b2 3 f3 fyyy
6
a23 fttt
6
h4 w4fyy b3 2 f2
2
ftt a32
2fy b3 ( )fy b2 f ft a2 b3 f fyt a3
w2fyy b1 2 f2
2b1 f fyt a1
ftt a12
2
w3 fy b2 ( )fy b1 f ft a1ftt a22
2
fyy b2 2 f2
2b2 f fyt a2 h3
( )w4 ( )fy b3 f ft a3 w2 ( )fy b1 f ft a1 w3 ( )fy b2 f ft a2 h2
( )w1 f w2 f w4 f w3 f h
Donde las funciones y las derivadas parciales están evaluadas en tn.
Deducir los valores de los parámetros para eliminar términos en la expresión para el error de la
aproximación resulta complejo.
Algoritmos numéricos 25
Profesor Leopoldo Silva Bijit 20-01-2010
La generación automática de la función , se logra con el segmento: > restart;
> tmv:=mtaylor(f(y,t),[y=yn,t=tn],4):
> sus1:={D[1,1](f)(yn,tn)=fyy,D[1,2](f)(yn,tn)=fyt,D[2](f)(yn,tn)=ft,
D[1](f)(yn,tn)=fy,D[2,2](f)(yn,tn)=ftt,f(yn,tn)=f,
D[1,1,2](f)(yn,tn)=fyyt,D[1,2,2](f)(yn,tn)=fytt,
D[2,2,2](f)(yn,tn)=fttt,D[1,1,1](f)(yn,tn)=fyyy}:
> tmv:=subs(sus1,tmv):
> k1:=h*f:
> k2:=h*eval(tmv,{y-yn=b1*k1,t-tn=a1*h}):collect(k2,h):
> k3:=h*eval(tmv,{y-yn=b2*k2,t-tn=a2*h}):k3:=collect(k3,h):
> temp3:=k3:
for i from 5 to 13 do
temp3:=eval(temp3,h^i=0)
od:
collect(temp3,h):
> k4:=h*eval(tmv,{y-yn=b3*k3,t-tn=a3*h}):
k4:=collect(k4,h):
temp4:=k4:
for i from 5 to 40 do
temp4:=eval(temp4,h^i=0)
od:
> temp:=collect(w1*k1+w2*k2+w3*k3+w4*k4,h):
> for i from 5 to 40 do
temp:=eval(temp,h^i=0)
od:
Phi:=collect(temp,h):
La generación del primer término del error, puede generarse automáticamente con: > restart;
> y1:=taylor(y(t),t=tn,5):y1:=subs({t-tn=h},y1):
> sus0:={y(tn)=0,D(y)(tn)=f,`@@`(D,2)(y)(tn)=d2,`@@`(D,3)(y)(tn)=d3,`@@`
(D,4)(y)(tn)=d4}:
F:=subs(sus0,y1):
> d2:=diff(f(y(t),t),t):
> d3:=diff(d2,t):
> d4:=diff(d3,t):
> sus1:={D[1,1](f)(y(t),t)=fyy,D[1,2](f)(y(t),t)=fyt,D[2](f)(y(t),t)=ft,
D[1](f)(y(t),t)=fy,
D[2,2](f)(y(t),t)=ftt,f(y(t),t)=f,D[1,1,1](f)(y(t),t)=fyyy,D[1,1,2](f)
(y(t),t)=fyyt,
D[1,2,2](f)(y(t),t)=fytt,D[2,2,2](f)(y(t),t)=fttt}:
sus2:={diff(y(t),t)=f}:
sus3:={diff(y(t),`$`(t,2))=fy*f+ft}:
sus4:={diff(y(t),`$`(t,3))=(fyy*f+fyt)*f+fy*(fy*f+ft)+fyt*f+ftt}:
> d2:=subs(sus1,d2):d2:=subs(sus2,d2):
d3:=subs(sus1,d3):d3:=subs(sus3,d3):d3:=subs(sus2,d3):d3:=expand(d3):
26 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
d4:=subs(sus1,d4):d4:=subs(sus4,d4):d4:=subs(sus3,d4):d4:=subs(sus2,d4
):d4:=expand(d4):
> F;
Efectuando la comparación de los coeficientes del error, de tal forma de anular hasta la potencia
cuarta de h, se obtienen 19 ecuaciones:
Para eliminar el término del error proporcional a h:
:= e1 w1 w2 w3 w4 1
Para eliminar el término del error proporcional a h2:
:= e2 w4 b3 w2 b1 w3 b21
2
:= e3 w4 a3 w2 a1 w3 a21
2
Para eliminar el término del error proporcional a h3:
:= e4w2 b1 2
2
w4 b3 2
2
w3 b2 2
2
1
6
:= e5 w4 b3 a3 w2 b1 a1 w3 b2 a21
3
:= e6 w4 b3 b2 w3 b2 b11
6
:= e7 w4 b3 a2 w3 b2 a11
6
:= e8w4 a32
2
w2 a12
2
w3 a22
2
1
6
Para eliminar el término del error proporcional a h4:
:= ec9w2 b1 3
6
w4 b33
6
w3 b2 3
6
1
24
:= ec10w2 b1 2 a1
2
w4 b32 a3
2
w3 b22 a2
2
1
8
:= ec11 w4 b32 b21
2b3 b22 w3
1
2b2 b12 b22 b1
1
6
:= ec12 w4 b32 a2 w3 b22 a11
8
:= ec13w2 b1 a12
2
w4 b3 a32
2
w3 b2 a22
2
1
8
Algoritmos numéricos 27
Profesor Leopoldo Silva Bijit 20-01-2010
:= ec14 w4 ( )b3 b2 a2 b3 b2 a3 w3 ( )b2 b1 a1 b2 b1 a25
24
:= ec15 w4 b3 a2 a3 w3 b2 a1 a21
8
:= ec16 w4 b3 b2 b11
24
:= ec17 w4 b3 b2 a11
24
:= ec18w4 b3 a22
2
w3 b2 a12
2
1
24
:= ec19w2 a13
6
w4 a33
6
w3 a23
6
1
24 Se tienen 19 ecuaciones para 10 incógnitas, sin embargo considerando que:
Las ecuaciones 16 y 17, implican que a1 = b1.
Las ecuaciones 6 y 7, con a1 = b1, implican: a2 = b2
Las ecuaciones 2 y 3, con a1 = b1, y a2 = b2 implican: a3 = b3
Las ecuaciones 12 y 15 implican, con a1 = b1, a2 = b2 que a3 = b3. Por lo tanto no se requieren
las ecuaciones 12 y 15.
Las ecuaciones 9, 10, 13 y 19, con a1 = b1, a2 = b2 y a3 = b3 son idénticas. Sólo se requiere
considerar una de ellas.
Las ecuaciones 4, 5 y 8, con a1 = b1, a2 = b2 y a3 = b3 son idénticas. Sólo se requiere considerar
una de ellas.
Las ecuaciones 12 y 11 implican la ecuación 18. Por lo tanto no se requiere emplear la ecuación
11.
Las ecuaciones 12 y 14 implican la ecuación 18. Por lo tanto no se requiere emplear la ecuación
14.
Lo anterior reduce el sistema a 10 ecuaciones en 10 incógnitas (a1, a2, a3, b1, b2, b3, w1, w2,
w3, w4).
Las 10 ecuaciones son las numeradas: 1, 2, 3, 6, 7, 16, 17, 4, 9, 18.
> ecc:={e1,e2,e3,e6,e7,ec16,ec17,e4,ec9,ec18}:
> sol:=[solve(ecc)]:sol[1];
{ }, , , , , , , , ,b11
2a2
1
2w1
1
6a3 1 b3 1 a1
1
2w3
1
3w2
1
3b2
1
2w4
1
6
Puede comprobarse que la solución satisface las 19 ecuaciones, y que el error queda:
28 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
5( )nT O h
Entonces el algoritmo de Runge Kutta de cuarto orden resulta:
1
12
23
4 3
1 1 2 3 4
,
,2 2
,2 2
,
12 2
6
n n
n n
n n
n n
n n
k hf y t
k hk hf y t
k hk hf y t
k hf y k t h
y y k k k k
Se calculan cuatro valores de la función, para obtener el siguiente punto.
Existen numerosos algoritmos en los que se varía el intervalo de tiempo entre puntos.
23.3. Solución de ecuación no lineal.
El problema consiste en encontrar la raíz de la ecuación no lineal: ( ) 0f x
Normalmente la solución de ( ) 0f x , puede ser difícil de encontrar analíticamente, pero como
veremos es sencilla de resolver iterativamente.
23.3.1. Método de Newton-Raphson.
Para resolver ( ) 0f x , se parte de un valor 0x y se genera una serie de iteraciones ix que se
acerquen a la solución sx , donde ( ) 0sf x .
En cursos de análisis numérico se responden las preguntas: ¿Cuándo la secuencia ix converge a
la solución correcta? ¿Cuán rápido se converge? ¿La convergencia depende del intento inicial
0x ? ¿Cuándo detener las iteraciones?.
El método de Newton-Raphson consiste en reemplazar, mediante la expansión de Taylor, la
función por su versión lineal, en torno a la solución:
( ) ( ) ( )( )s s s
dff x f x x x x
dx
Para un punto cualquiera se obtiene:
Algoritmos numéricos 29
Profesor Leopoldo Silva Bijit 20-01-2010
1 1( ) ( ) ( )( )k k k k k
dff x f x x x x
dx
Efectuando: 1( ) 0kf x , se obtiene la fórmula de la iteración de Newton-Raphson,
despejando 1kx :
1
( )
( )
kk k
k
f xx x
dfx
dx
Podemos interpretar la fórmula de la iteración, planteando la relación anterior en 0x , y
calculando 1x . Situación que se ilustra en la Figura 23.5.
x0 x1
f(x)
x2
f(x0)
xs 0
f(x1)
Figura 23.5. Iteración Newton-Raphson.
Resulta, de la interpretación gráfica de la derivada en 0x :
00 0
0 1
( )( ) ( )
f xdftg x
dx x x
Despejando 1x , se obtiene el primer valor de aproximación del método de Newton-Raphson:
Nótese que 1( )f x no es cero, lo cual implica que 1x es una aproximación de sx . También debe
notarse que para calcular la siguiente aproximación deben calcularse la función y la derivada en
el punto anterior.
El proceso debe repetirse hasta que: 1k kx x tolerancia
Donde el valor de tolerancia debe ser un valor lo suficientemente pequeño, para que la solución
se considere aceptable. Con números reales de precisión simple (float en C), un valor razonable
1
1 0 0 0( ) ( )df
x x x f xdx
30 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
de tolerancia es 10-6
, que es el valor del número real más pequeño representable, en el formato
interno normalizado IEEE754.
Si el valor inicial es adecuado conviene limitar el número máximo de iteraciones, de este modo
si no existe convergencia se asegura que el algoritmo termine.
También puede verificarse que la ordenada en los puntos sucesivos esté dentro de cierto rango:
1( )kf x tolerancia
Para evitar oscilaciones o ciclos no convergentes se limita el número de iteraciones. Pueden
producirse en un caso como el que se ilustra en la Figura 23.6.
Figura 23.6. Oscilación.
Se emplea un intervalo [x1..x2] en el cual se busca la raíz, para evitar la no convergencia debido
a máximos o mínimos dentro del intervalo. Si en una iteración se encuentra un punto con
derivada casi horizontal en el intervalo, el valor para el nuevo x se aleja de la raíz, como se
ilustra en la Figura 23.7. Para prevenir esta situación se verifica que la nueva aproximación de la
raíz permanezca dentro del intervalo.
Figura 23.7. Divergencia.
Se detienen las iteraciones si los dos últimos valores obtenidos difieren en determinada
tolerancia.
La siguiente función intenta encontrar una raíz en el intervalo [x1, x2].
El primer argumento es un puntero a función con un argumento flotante y que retorna un
flotante. Cuando se invoca a NewtonRaphson el argumento actual es el nombre de la función
que calcula el cuociente de la función con su derivada evaluada en el x actual.
Algoritmos numéricos 31
Profesor Leopoldo Silva Bijit 20-01-2010
#define nmax 20 //Máximo número de Iteraciones.
float NewtonRaphson(float (*pfuncion)(float), float x1, float x2, float tolerancia)
{
int k;
float dx, x;
x=0.5*(x1+x2); //Intento inicial.
for (k=1; k<=nmax; k++)
{
dx=(*pfuncion)(x); //se invoca a la función cuyo nombre se pasa como primer argumento.
x -= dx; //nuevo valor
if ((x<x1) || (x>x2))
{printf("Salta fuera del intervalo");
exit(1);
}
if (fabs(dx) < tolerancia) return (x); //Converge.
}
printf("No converge en %d iteraciones.\n", nmax);
exit(1);
return (0.0); //Para evitar warning.
}
La siguiente función, para un polinomio, define la derivada y retorna el cuociente entre la
función y su derivada. Nótese que la función recibe un argumento flotante y retorna un flotante.
De este modo, su nombre fin2 es un puntero constante compatible con el primer argumento de
NewtonRaphson.
float fun2(float x)
{ float funcion, derivada;
funcion=.2*pow(x,4)-2*pow(x,3)+x+7; //función de x, evaluada en x;
derivada=.8*pow(x,3)-6*pow(x,2)+1; //derivada de f evaluada en x;
return (funcion/derivada);
}
La función: 4 3( ) 0.2 2 7f x x x x , tiene dos raíces reales, como se ilustra en la Figura
23.8.
La solución con fsolve de Maple, entrega: 1.742869875, 9.913192964
32 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Figura 23.8.
Empleando la función NewtonRaphson, pueden obtenerse las raíces mediante:
printf("Valor de la raíz = %f\n", NewtonRaphson(fun2,1,2,1e-6 ));
printf("Valor de la raíz = %f\n", NewtonRaphson(fun2,9,11,1e-6 ));
23.3.2. Generalización para sistemas de ecuaciones no lineales.
Para un sistema de ecuaciones no lineales, se emplea la expansión de Taylor para varias
variables.
La expansión es una linealización en torno a la solución:
( ) ( ) ( )( )s s sF x F x J x x x
Las cantidades ( )F x y ( )sx x se expresan como vectores, y ( )sJ x como una matriz,
denominada Jacobiano.
Para un punto cualquiera, con aproximación de primer orden, se tiene:
1 1( ) ( ) ( )( )k k k k kF x F x J x x x
Para entender la relación anterior se ilustra la forma que ella toma para dos funciones de dos
variables x1 y x2, se tiene:
1
2
1 1
1 1 1 1
1 1 2 12 2
( 1 , 2 ) ( 1 , 2 )
( 1 , 2 ) ( 1 , 2 ) 1 11 2
( 1 , 2 ) ( 1 , 2 ) 2 2( 1 , 2 ) ( 1 , 2 )
1 2
k k k k
k k k k k k
k k k k k kk k k k
F x x F x x
F x x F x x x xx x
F x x F x x x xF x x F x x
x x
Una explicación del cambio de la función de dos variables, puede efectuarse considerando el
plano tangente a la superficie, en el punto (x10, x20) que pasa también por el punto (x11, x21).
Algoritmos numéricos 33
Profesor Leopoldo Silva Bijit 20-01-2010
Donde el punto 0 es el inicial, y el punto 1, se obtiene pasando un plano tangente a la superficie
en el punto 0.
x10
x11
x20
x21 x1
x2
F1x1
F1x2
Figura 23.9. Interpretación del Jacobiano de dos variables.
Aplicando la interpretación geométrica de las derivadas parciales, se tienen:
1 0 0 1 11
0 1
1 0 0 1 22
0 1
( 1 , 2 )( )
1 1 1
( 1 , 2 )( )
2 2 2
xx
xx
F x x Ftg
x x x
F x x Ftg
x x x
El cambio total de la función, resulta:
1 0 0 1 0 01 1 1 2 0 1 0 1
( 1 , 2 ) ( 1 , 2 )( 1 1 ) ( 2 2 )
1 2x x
F x x F x xF F x x x x
x x
Aplicando el método de Newton-Raphson, que consiste en asumir que el plano tangente pasa
por el punto que es una aproximación a la solución. Esto equivale a efectuar:
1 1 1
2 1 1
( 1 , 2 )
( 1 , 2 )0k k
k k
F x x
F x x
Entonces la fórmula de iteración, resulta:
1 1
1 1
1 22 2
( 1 , 2 ) ( 1 , 2 )
1 1 ( 1 , 2 )1 2
2 2 ( 1 , 2 )( 1 , 2 ) ( 1 , 2 )
1 2
k k k k
k k k k
k k k kk k k k
F x x F x x
x x F x xx x
x x F x xF x x F x x
x x
Finalmente, despejando el nuevo punto:
34 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
1 1
1 1
1 22 2
1( 1 , 2 ) ( 1 , 2 )
1 1 ( 1 , 2 )1 2
2 2 ( 1 , 2 )( 1 , 2 ) ( 1 , 2 )
1 2
k k k k
k k k k
k k k kk k k k
F x x F x x
x x F x xx x
x x F x xF x x F x x
x x
La que expresada en términos de vectores y la matriz inversa del Jacobiano, resulta en general,
para n variables:
1
1 ( ) ( ) k k k kx x J x F x
Una mejor visualización de la suma de los incrementos, se logra observando los triángulos
semejantes en la Figura 23.10.
Por el punto inicial (2, 2, 10) se pasa el plano z=2x+3y que también pasa por el punto (0, 0, 0).
Se han dibujado además los planos de z constante, z=4 y z=6.
2, 3z z
x y , 64
z zx y
x y
Figura 23.10. Variación total de función de dos variables.
Volviendo al caso de dos variables, considerando el álgebra de matrices, se tiene:
11a b x by dx
c d y cx ayad bc
Entonces las fórmulas de iteración de Newton-Raphson para un sistema de ecuaciones no
lineales de dos variables, resultan:
1 22 1
11 2 1 2
( ) ( )( ( ) ( ))
2 21 1( ) ( ) ( ) ( )
1 2 2 1
k k
F k F kF k F k
x xx xF k F k F k F k
x x x x
Algoritmos numéricos 35
Profesor Leopoldo Silva Bijit 20-01-2010
2 11 2
11 2 1 2
( ) ( )( ( ) ( ))
2 12 2( ) ( ) ( ) ( )
1 2 2 1
k k
F k F kF k F k
x xx xF k F k F k F k
x x x x
En caso de mayores órdenes debe invertirse el Jacobiano, o alternativamente resolverse el
sistema lineal de ecuaciones, para las incógnitas 1kx :
1( )( ) ( ) k k k kJ x x x F x
Donde el vector de las funciones y el Jacobiano deben evaluarse para cada kx
Referencias.
Numerical Recipes In C: The Art of Scientific Computing. Cambridge University Press. 1992.
36 Estructuras de Datos y Algoritmos
Profesor Leopoldo Silva Bijit 20-01-2010
Índice general.
CAPÍTULO 23 ............................................................................................................................................ 1
ALGORITMOS NUMÉRICOS. ................................................................................................................ 1
23.1. SOLUCIÓN DE SISTEMA SIMULTÁNEO DE ECUACIONES LINEALES. ..................................................... 1 23.1.1. Descomposición LU. ................................................................................................................ 1 23.1.2. Métodos iterativos. ................................................................................................................... 7
23.2. SOLUCIÓN NUMÉRICA DE SISTEMAS DE ECUACIONES DIFERENCIALES. ........................................... 10 23.2.1. Formulación de ecuaciones de estado. .................................................................................. 11 23.2.2. Método de Euler. .................................................................................................................... 11 23.2.3. Algoritmo de Euler. ................................................................................................................ 13 23.2.4. Algoritmo trapezoidal. ........................................................................................................... 14 23.2.5. Algoritmo de Simpson. ........................................................................................................... 15 23.4.6. Métodos multietapas. ............................................................................................................. 17
23.4.6.1. Método de Heun. ............................................................................................................................. 20 23.4.6.2. Método del punto medio. ................................................................................................................ 20 23.4.6.3. Método de Ralston. ......................................................................................................................... 20
23.4.7. Métodos de Runge-Kutta de cuarto orden. ............................................................................ 21 23.3. SOLUCIÓN DE ECUACIÓN NO LINEAL. .............................................................................................. 28
23.3.1. Método de Newton-Raphson. ................................................................................................. 28 23.3.2. Generalización para sistemas de ecuaciones no lineales. ..................................................... 32
REFERENCIAS. ......................................................................................................................................... 35 ÍNDICE GENERAL. .................................................................................................................................... 36 ÍNDICE DE FIGURAS. ................................................................................................................................ 37
Algoritmos numéricos 37
Profesor Leopoldo Silva Bijit 20-01-2010
Índice de figuras.
FIGURA 23.1. MÉTODO DE EULER. .............................................................................................................. 14 FIGURA 23.2. MÉTODO TRAPEZOIDAL. ........................................................................................................ 14 FIGURA 23.3. MÉTODO DE SIMPSON. .......................................................................................................... 15 FIGURA 23.4. MÉTODO DEL PUNTO MEDIO. ................................................................................................. 20 FIGURA 23.5. ITERACIÓN NEWTON-RAPHSON. ........................................................................................... 29 FIGURA 23.6. OSCILACIÓN. ......................................................................................................................... 30 FIGURA 23.7. DIVERGENCIA....................................................................................................................... 30 FIGURA 23.8. .............................................................................................................................................. 32 FIGURA 23.9. INTERPRETACIÓN DEL JACOBIANO DE DOS VARIABLES. ........................................................ 33 FIGURA 23.10. VARIACIÓN TOTAL DE FUNCIÓN DE DOS VARIABLES. .......................................................... 34
Top Related