Curso Propedeutico

Maestrı́a en Ingenierı́a Eléctrica

Métodos Numéricos

Dr. Félix Calderon Solorio

7 de noviembre de 2016
Índice general

Introducción 1
1.1. Introducción (Variables y Operadores) . . . . . . . . . . . . . . . . . . . . . 1
1.1.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2. Ayuda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.3. Tipos de Datos y Variables . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.4. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.5. Operadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2. Instrucciones Secuenciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3. Instrucciones Condicionales . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4. Instrucciones de Repetición . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.1. Ciclos for/end . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.3. Ciclos while/end . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.4. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.5. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.6. Ejemplo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.7. Ejemplo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5. Manejo de Matrices y Vectores . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.5.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.5.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.5.3. Arreglos Bidimensionales . . . . . . . . . . . . . . . . . . . . . . . . 16
1.6. Estructuras de Programa y Funciones . . . . . . . . . . . . . . . . . . . . . 19
1.6.1. Funciones que devuelven una sola variable . . . . . . . . . . . . . . . 19
1.6.2. Funciones que devuelven mas de una variable . . . . . . . . . . . . . 19
1.6.3. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6.4. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.6.5. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Sist. Lineales 25
2.1. Suma, resta, multiplicación y división de matrices . . . . . . . . . . . . . . . 25
 ÍNDICE GENERAL

2.1.1. Suma de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.1.2. Resta de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.3. Multiplicación de matrices . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2. Factorización triangular LU . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1. Sustitución hacia adelante . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.2. Sustitución hacia atrás . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.3. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.4. Factorización de Crout . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.5. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.6. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.2.7. Cálculo de la Matriz inversa utilizando descomposición LU . . . . . 36
2.2.8. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3. Eliminación Gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.4. Gauss-Jordan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.4.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.4.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.4.3. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.5. Inversa de Shiplay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
2.5.1. Ejemplo: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
2.6. Sistemas dispersos y estrategias para conservarla . . . . . . . . . . . . . . . 62
2.7. Método iterativo de Jacobi y Gauss-Seidel . . . . . . . . . . . . . . . . . . . 64
2.7.1. Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
2.7.2. Algoritmo iterativo de Gauss-Seidel . . . . . . . . . . . . . . . . . . 66
2.7.3. Ejemplo matrices dispersas . . . . . . . . . . . . . . . . . . . . . . . 68

Sist. no lineales 71
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.2. Métodos de punto fijo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.3. El método de Bisecciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.3.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.4. Método de Regula Falsi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.4.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.4.2. Solución de un circuito con un diodo . . . . . . . . . . . . . . . . . . 80
3.5. El método Newton Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.5.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.5.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.5.3. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
3.6. Convergencia del método de Newton Raphson . . . . . . . . . . . . . . . . . 92

Ajuste de curvas 95
4.1. Regresión lineal por el método de mı́nimos cuadrados . . . . . . . . . . . . 95
ÍNDICE GENERAL

4.1.1. Ajuste por mı́nimos cuadrados . . . . . . . . . . . . . . . . . . . . . 96

4.1.2. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.1.3. Regresión polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.1.4. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.2. Interpolación lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.2.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.3. Interpolación cuadrática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.3.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.3.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.4. Formulas de interpolación de Newton . . . . . . . . . . . . . . . . . . . . . . 107
4.4.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4.4.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4.5. Interpolación de Polinomios de Lagrange . . . . . . . . . . . . . . . . . . . . 110
4.5.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Diferenciación e Integración 113

5.1. Derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.2. Integración por el método de barras . . . . . . . . . . . . . . . . . . . . . . 113
5.2.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.3. Integración por el método de trapezoides . . . . . . . . . . . . . . . . . . . . 115
5.3.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.4. Integración por el método de regla Simpson 1/3 . . . . . . . . . . . . . . . . 116
5.4.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.5. Integración por el método de regla Simpson 3/8 . . . . . . . . . . . . . . . 118
5.5.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
5.6. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
5.6.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
5.6.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Ecuaciones diferenciales ordinarias 121

6.1. Integración por el método de Euler . . . . . . . . . . . . . . . . . . . . . . . 121
6.2. Integración por el método de Heun con solo uno y con varios predictores . . 122
6.3. Integración por el método del punto medio . . . . . . . . . . . . . . . . . . 123
6.4. Runge-Kutta 2do orden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
6.5. Runge-Kutta 3er orden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
6.6. Runge-Kutta 4to orden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
6.7. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.7.1. Una ecuación diferencial sencilla . . . . . . . . . . . . . . . . . . . . 127
6.7.2. Circuito RL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
6.7.3. Sistema Masa Resorte . . . . . . . . . . . . . . . . . . . . . . . . . . 133
ÍNDICE GENERAL
Introducción a la Programación en Matlab

1.1. Introducción (Variables y Operadores)

1.1.1. Introducción

Matlab puede considerarse como un lenguaje de programación tal como C, Fortran, Java,
etc. Algunas de las caracterı́sticas de Matlab son:
La programación es mucho mas sencilla.
Hay continuidad entre valores enteros, reales y complejos.
La amplitud del intervalo y la exactitud de los números es mayor.
Cuenta con una biblioteca matemática amplia.
Abundantes herramientas gráficas, incluidas funciones de interfaz gráfica con el usua-
rio.
Capacidad de vincularse con los lenguajes de programación tradicionales.
Transportabilidad de los programas.
Algunas de sus desventajas son:
Necesita de muchos recursos de sistema como son Memoria, tarjeta de videos, etc.
para funcionar correctamente.
El tiempo de ejecución es lento.
No genera código ejecutable.
Es caro.
En una estación de trabajo UNIX, MATLAB puede abrirse tecleando
¿matlab

1
2 INTRODUCCIÓN

En el caso de sistemas LINUX, existe una versión similar al MATLAB , llamado OCTAVE,
el cual tiene las mismas funciones y desempeño que el MATLAB. En este caso para iniciar
la sesión se da:
¿octave
En Macintosh o Windows, haga clic en el icono de MATLAB.

1.1.2. Ayuda

Si no entiende el significado de un comando, teclee help y el nombre del comando que desea
revisar. Este comando desplegara información concisa respecto que será de utilidad para
el usuario. Ası́ por ejemplo cuando se da el comando
)) help help
Se despliega en la pantalla.
HELP On-line help, display text at command line. HELP, by itself, lists all primary help
topics. Each primary topic corresponds to a directory name on the MATLABPATH.
”HELP TOPIC”gives help on the specified topic. The topic can be a command name, a
directory name, or a MATLABPATH relative partial pathname (see HELP PARTIAL-
PATH). If it is a command name, HELP displays information on that command. If it is
a directory name, HELP displays the Table-Of-Contents for the specified directory. For
example, ”help [Link] ”help matlab/general”both list the Table-Of-Contents for the
directory toolbox/matlab/general.
HELP FUN displays the help for the function FUN.
T = HELP(’topic’) returns the help text in a \ n separated string.
LOOKFOR XYZ looks for the string XYZ in the first comment line of the HELP text in
all M-files found on the MATLABPATH. For all files in which a match occurs, LOOKFOR
displays the matching lines.
MORE ON causes HELP to pause between screenfuls if the help text runs to several
screens.
In the online help, keywords are capitalized to make them stand out. Always type com-
mands in lowercase since all command and function names are actually in lowercase.
For tips on creating help for your m-files ’type help.m’.
See also LOOKFOR, WHAT, WHICH, DIR, MORE.
1.1. INTRODUCCIÓN (VARIABLES Y OPERADORES) 3

El comando what (Qué) muestra una lista de archivos M, MAT y MEX presentes en
el directorio de trabajo. Otra manera de realizar esta operación es utilizar el comando
dir.

El comando who (Quien) lista las variables utilizadas en el espacio de trabajo actual.

1.1.3. Tipos de Datos y Variables

No es necesario declarar los nombres de la variables ni sus tipos. Esto se debe a que
los nombres de la variables en Matlab no son diferentes para los números enteros, reales
y complejos. Sin embargo esto da lugar a conflictos cuando se utilizan los nombres de
variables de palabras reservadas para:

1. nombres de variables especiales.

Nombre de variable Significado

eps Épsilon de la máquina (2.2204e-16)
pi pi (3.14159...)
iyj Unidad imaginaria
inf infinito
NaN no es número
date fecha
flops Contador de operaciones de punto flotante
nargin Número de argumentos de entrada a una función
nargout Número de argumentos de salida de una función

2. nombres de funciones.

Trigonométricas
4 INTRODUCCIÓN

Nombre de variable Significado

sin Seno
sinh Seno Hiperbólico
asin Seno Inverso
asinh Seno hiperbólico inverso
cos Coseno
cosh Coseno Hiperbólico
acos Coseno inverso
acosh Coseno hiperbólico inverso
tan Tangente
tanh Tangente Hiperbólica
atan Tangente inversa.
atan2 Tangente inversa en cuatro cuadrantes
atanh Tangente hiperbólica inversa
sec Secante
sech Secante Hiperbólica
asec Secante inversa
asech Secante hiperbólica Inversa
csc Cosecante
csch cosecante hiperbólica
acsc Cosecante Inversa
acsch Inverse hyperbolic cosecant
cot Cotangent
coth Hyperbolic cotangent
acot Inverse cotangent
acoth Inverse hyperbolic cotangent

Exponencial.

Nombre de variable Significado

exp Exponencial
log Logaritmo Natural
log10 logaritmo común (base 10)
log2 Logaritmo base 2
pow2 Potenciacion en Base 2
sqrt Raı́z Cuadrada
nextpow2 Próxima potencia de 2

Complejos.
1.1. INTRODUCCIÓN (VARIABLES Y OPERADORES) 5

Nombre de variable Significado

abs Valor Absoluto
angle Ángulo de Fase
complex Constructor de números complejos
conj Conjugado de un complejo
imag Parte imaginaria
real Parte real
unwrap Desenvolvimiento de fase
isreal Verifica si un arreglo es real
cplxpair Ordena números complejos

Redondeo y residuos.

Nombre de variable Significado

fix Parte entera de un número
floor Parte fraccionaria
ceil Parte entere de un número
round Redondea al siguiente entero
mod Residuo con signo de una división
rem Residuo sin signo de una división
sign Signo

3. nombres de comandos.

Comandos

Nombre de comando Significado

what Lista archivos en el directorio
dir Lista archivos en el directorio
who Variables utilizadas
clear Borra variables utilizadas
etc etc

1.1.4. Ejemplo

Calcular el volumen en una esfera.

clear;
r = 2
vol = (4/3)*pi*r^3
6 INTRODUCCIÓN

1.1.5. Operadores

Los operadores aritméticos como +, -, *, / son los mismos que en los lenguajes tradicionales
como C, Java, Fortran, etc., ası́ como la precedencia de estos.
Operador Simbolo
suma +
resta -
multiplicación *
división /
reciproco \
Un operador no convencional es el reciproco
clear;
c = 3 \ 1
Los operadores condicionales que existen en Matlab son:
Operador Simbolo
Mayor que >
Menor que <
Mayor igual >=
Menor igual <=
Igual a ==
Diferente de ˜=
Los cuales son utilizados para hacer condicionales con la sentencia if
clear;
r = -2
if r > 0, (4/3)*pi*r\^3
end
Note que todas las sentencias if deben ir acompañadas por un end.
Los operadores lógicos and y or se denotan con & y k
g = 5
if g>3 | g <0, a = 6
end
y se puede hacer agrupamiento con los operadores & y k
clear;
a = 2
b = 3
1.2. INSTRUCCIONES SECUENCIALES 7

c = 5
if((a==2 | b==3) & c < 5) g=1; end;

1.2. Instrucciones Secuenciales

1.3. Instrucciones Condicionales

Para hacer condicionales se utiliza la sentencia if, la cual, siempre debe ir terminada con
la sentencia end
r = -2
if r > 0, (4/3)*pi*r^3
end
Los enunciados lógicos and y or se denotan con & y k, pueden ser utilizados para imple-
mentar condicionales de la manera siguiente.
clear;
g = 5
if g>3 | g <0, a = 6
end
pause;
además los operadores & y k se puede agrupar como
clear;
a = 2
b = 3
c = 5

if((a==2 | b==3) & c < 5) g=1; end;

pause;
Las condicionales if se pueden utilizar con else o elseif
clear;

r = 2;

if r > 3 b = 1;
elseif r == 3 b = 2;
8 INTRODUCCIÓN

else b = 0;
end;
Ver código ejem001.m

1.4. Instrucciones de Repetición

En Matlab existen dos maneras de implementar ciclos. La primera con los comandos for/end
y la segundo con los comandos while/end, de manera muy similar a los lenguajes de alto
nivel.

1.4.1. Ciclos for/end

La sintaxis de este comando es

for r=inicio: incremento: fin
instrucciones_a_repetir
instrucciones_a_repetir
instrucciones_a_repetir
instrucciones_a_repetir
instrucciones_a_repetir
end;
imprimir los números del 1 a 100 se hace :
for x=1: 1:100

end;
El siguiente conjunto de instrucciones realiza una cuenta de 100 a 80 con decrementos de
0.5.
for x=100:-0.5: 80

end;
en el caso de decrementos o incrementos unitarios, se puede omitir el valor del incremen-
to.
1.4. INSTRUCCIONES DE REPETICIÓN 9

for x=1: 1:100, x, end

1.4.2. Ejemplo 1

Utilizando el comando for/end, calcular el volumen de cinco esferas de radio 1, 2, 3, 4 y 5

se hace:
for r=1:5
vol = (4/3)*pi*r^3;
disp([r, vol])
end;
Los ciclos pueden hacerse anidados de la siguiente manera.
for r=1:5
for s=1:r
vol = (4/3)*pi*(r^3-s^3);
disp([r, s, vol])
end
end
Podemos utilizar el comando break para detener la ejecución de un ciclo
for i=1:6
for j=1:20
if j>2*i, break, end
disp([i, j])
end
end
Ver código esferas.m

1.4.3. Ciclos while/end

La sintaxis de esta comando es

while condición
instrucciones_a_repetir
end
Ası́ por ejemplo podemos implementar al igual que con los ciclos for/end, un pequeño
programa que imprima los números del 1 al 100.
10 INTRODUCCIÓN

x = 1;
while x <= 100
x
x = x + 1;
end
Ver ejemplo while.m
El ejemplo para desplegar el volumen de una esfera con radios de 1, 2, 3, 4 y 5 queda.
r = 0;
while r<5
r = r+1;
vol = (4/3)*pi*r^3;
disp([r, vol])
end;
otro ejemplo interesante es:
clear;
r = 0
while r<10
r = input(’Teclee el radio (o -1 para terminar): ’);
if r< 0, break, end
vol = (4/3)*pi*r*3;
fprintf(’volumen = \%7.3f\n’, vol)
end
Ver esferas while.m

1.4.4. Ejemplo 2

Hacer un programa que permita imprimir un triángulo rectángulo formado por asteris-
cos.
% Codigo para imprimir un triangulo

fprintf(’\ntriangulo\n\n’)

for k=1:7
for l=1:k
fprintf(’*’)
end;
1.4. INSTRUCCIONES DE REPETICIÓN 11

fprintf(’\n’)
end
Ver triangulo.m

1.4.5. Ejemplo 3

Hacer un programa para desplegar un rectángulo de base 6 y altura 7.

% Codigo para imprimir un rectangulo

fprintf(’\nrectangulo\n\n’)

for k=1:7
for l=1:6
fprintf(’*’)
end;
fprintf(’\n’)
end
Ver rectangulo.m

1.4.6. Ejemplo 4

Hacer un programa para imprimir un pino utilizando un solo carácter.

a=10; %altura del follaje del pino
n=12; %Posición horizontal del vértice.
t=3; %altura del tronco del pino
d=4; %diámetro del tronco del pino

% Dibujar el follaje del pino, de altura ’a’

for i=1:a
clear cad2 cad1
num_ast=2*i-1;
num_esp=n-i;
cad1(1:num_esp)=’ ’;
cad2(1:num_ast)=’*’;
fprintf(’\%s\%s\n’,cad1,cad2)
end
12 INTRODUCCIÓN

% Dibujar el tronco del pino, de altura ’t’

clear cad2 cad1

num_ast=d;
num_esp=n-d/2;
cad1(1:num_esp)=’ ’;
cad2(1:num_ast)=’I’;

for i=1:t
fprintf(’\%s\%s\n’,cad1,cad2)
end

Ver pino.m

1.4.7. Ejemplo 5

Hacer un programa para imprimir el triángulo de Pascal.

nr=8; % Numero de renglones del triangulo de Pascal.

n=15; % Numero de espacios en blanco antes del vértice.
x(1)=1;
cad1(1:n)=’ ’;
fprintf(’%s%3.0f\n\n’,cad1,x(1)); % vertice del triangulo
for k=2:nr-1;
clear cad1 cad2
num_esp=n-2*k+1;
cad1(1:num_esp)=’ ’;
clear x
x(1)=1;
for c=2:k;
x(c)=x(c-1)*(k-c+1)/(c-1);
end
fprintf(’\%s’,cad1)
for c=1:k
fprintf(’ \%3.0f’,x(c))
end
fprintf(’\n\n’)
end

Ver pascal.m
1.5. MANEJO DE MATRICES Y VECTORES 13

1.5. Manejo de Matrices y Vectores

Un vector de datos puede definirse como

x = [0, 0.1, 0.2, 0.3, 0.4, 0.5]
si se desea imprimir un dato en particular se teclea
x(3)
el cual imprimirá el número en la tercer posición del arreglo, el primer elemento se numera
con el uno.
Una forma equivalente de definir la misma x es
clear;
for i=1:6
x(i) = (i-1)*0.1;
end

x
x(3)
Otra forma de escribir un arreglo es
clear;

x = 2: -0.4: -2

pause;
Podemos definir arreglos en fila o columna
clear;

z = [0; 0.1; 0.2; 0.3; 0.4; 0.5]

z = [0, 0.1, 0.2, 0.3, 0.4, 0.5]’

podemos realizar operaciones entre arreglos fila o columna con
clear;

x = [1, 2, 3, 4]
y = [4, 3, 2, 1]
14 INTRODUCCIÓN

suma = x + y
resta = x - y
mult = x .* y
div = x ./ y
Las operaciones de suma y resta son iguales que para el álgebra lineal. En cambio .* y ./
son operadores nombrados para la división de arreglos. Si se omite el punto el significado
es diferente lo cual es equivalente a
clear;

x = [1, 2, 3, 4]
y = [4, 3, 2, 1]

Para la suma hacemos

for i=1:4; suma(i) = x(i) + y(i); end;
En la resta
for i=1:4; resta(i) = x(i) - y(i); end;
Multiplicación
for i=1:4; mult(i) = x(i) * y(i); end;
y División
for i=1:4; div(i) = x(i) / y(i); end;
Ver arreglos.m
Propiedades de los arreglos
Concatenación
clear;

x = [2, 3]
x = [x, 4]
en el caso de arreglos columna
clear;

y = [2, 3]’
y = [y; 4]
1.5. MANEJO DE MATRICES Y VECTORES 15

para extraer la una parte de un vector

clear;

y = [1,2,3,4,5,6,7,8,9]
w = y(3:6)
para obtener la longitud de un arreglo se utiliza
clear;

y = [1,2,3,4,5,6,7,8,9]
length(y)
La variables de cadena también puede tener caracteres
clear;

v = ’Hola Mundo’

w = [’H’, ’o’, ’l’, ’a’, ’ ’, ’M’, ’u’, ’n’, ’d’, ’o’]

y podemos cambiar el orden de impresión haciendo
clear;

v = ’Hola Mundo’
v = v’

w = [’H’, ’o’, ’l’, ’a’, ’ ’, ’M’, ’u’, ’n’, ’d’, ’o’]

w = w’
Ver arreglos propiedades.m

1.5.1. Ejemplo 1

Dado un arreglo de datos calcular el promedio de este y el mayor de los valores.

x = [3 9 5 8 2]
n = length(x);
suma =0;
max = x(1);
for k=1:n
suma = suma + x(k);
16 INTRODUCCIÓN

if(max < x(k)); max = x(k); end;

end;
suma = suma/n;

fprintf(’El promedio es = \%5.2f\n’, suma);

fprintf(’El mayor es = \%5.2f\n’, max);
Ver promedio mayor.m

1.5.2. Ejemplo 2

Escribir un programa que verifique si una cadena de caracteres e un palı́ndromo.

x = ’10011’
y = x(length(x):-1:1)
if(x == y)
fprintf(’Es palı́ndromo \n’);
else
fprintf(’No es palı́ndromo \n’);
end;
Ver palindroma.m

1.5.3. Arreglos Bidimensionales

Para definir arreglos bidimensionales o matrices hacemos

m = [0.1, 0.2, 0.3; 0.4, 0.5, 0.6; 0.7, 0.8, 0.9]
otra manera de definirlos es utilizar
clear;
m(1,1) = 0.1;
m(1,2) = 0.2;
m(1,3) = 0.3;
m(2,1) = 0.4;
m(2,2) = 0.5;
m(2,3) = 0.6;
m(3,1) = 0.7;
m(3,2) = 0.8;
m(3,3) = 0.9;
Podemos listar columnas de la matriz haciendo
1.5. MANEJO DE MATRICES Y VECTORES 17

m(:,1)
m(:,2)
m(:,3)
o también renglones
m(1,:)
m(2,:)
m(3,:)
podemos realizar operaciones de +. -, * y / con matrices
a = [0.1, 0.2, 0.3; 0.4, 0.5, 0.6; 0.7, 0.8, 0.9]
b = [0.3, 0.4, 1.3; 0.6, -0.7, 1.0; -2.0, 1.8, 9]

suma = a + b
resta = a - b
mult = a .* b
div = a ./ b
lo cual es equivalente a
a = [0.1, 0.2, 0.3; 0.4, 0.5, 0.6; 0.7, 0.8, 0.9]
b = [0.3, 0.4, 1.3; 0.6, -0.7, 1.0; -2.0, 1.8, 9]

for i=1:3
for j=1:3
suma(i,j) = a(i,j) + b(i,j);
end
end

%suma

for i=1:3
for j=1:3
resta(i,j) = a(i,j) - b(i,j);
end
end

%resta

for i=1:3
for j=1:3
mult(i,j) = a(i,j) * b(i,j);
18 INTRODUCCIÓN

end
end

%multiplicación

for i=1:3
for j=1:3
div(i,j) = a(i,j) / b(i,j);
end
end

%división

pause;
también podemos utilizar el operador de potenciación en arreglos
clear;

a = [0.1, 0.2, 0.3; 0.4, 0.5, 0.6; 0.7, 0.8, 0.9]

g = a .^2

pause;
el cual es equivalente
a = [0.1, 0.2, 0.3; 0.4, 0.5, 0.6; 0.7, 0.8, 0.9]

for i=1:3
for j=1:3
g(i,j) = a(i,j)^2;
end
end
Ver arreglos 2d.m
Algebra de Matrices
Cuando queremos realizar las operaciones del álgebra lineal de suma, resta, multiplicación
y división hacemos
la suma y resta son iguales pero la multiplicación cambia
A = [0.1, 0.2, 0.3; 0.4, 0.5, 0.6; 0.7, 0.8, 0.9]
1.6. ESTRUCTURAS DE PROGRAMA Y FUNCIONES 19

x = [1, 2, 3]’

b = A*x
para la división tendremos
A = [1, 4; 3, 5]

x = [2, 3]’

b = A\x

1.6. Estructuras de Programa y Funciones

Las funciones en MATLAB, que se guardan como archivos independientes, equivalen a las
subrutinas y funciones de otros lenguajes.

1.6.1. Funciones que devuelven una sola variable

Consideremos el ejemplo de escribir la función

f(x)= 2x3+4x2+ 3
podemos escribir en MATLAB un archivo con el nombre fun00.m como
function y = fun01(x)
y = 2*x^3+4*x^2+ 3;
Para ejecutar esta función, desde el ambiente de MATLAB podemos valuar f(2) como
fun01(2)

1.6.2. Funciones que devuelven mas de una variable

Una función puede devolver más de una variable y la sintaxis para escribir esta función
es
function [Y1, Y2, ..., Yn ] = fun_regresa_varias(X)
...
...
...
20 INTRODUCCIÓN

Y1 = ....
Y2 = ....
...
...

Supongamos que dado un conjunto de datos queremos realizar una función que devuelva la
media y la desviación estándar. Primero escribimos un archivo llamado fun02.m, que tenga
las siguientes instrucciones.

function [media, des] = fun02(x)

n = length(x);
suma = 0;
for k=1:n
suma = suma + x(k);
end;
media = suma/n;
suma = 0;
for k=1:n
suma = (x(k) - media)^2;
end;
des = sqrt(suma/n );

Guardamos el archivo y ejecutamos desde el MATLAB con:

[m d] = fun02(x)

y la ejecución regresa.

m =

d =

0.5774

Note que la función recibe dos argumentos a los que llamamos m y d. De hacer el llamado
de la función sin poner estos dos, no se genera error alguno pero solo se imprime el primer
parámetro que devuelve la función.

fun02(x)

ans =

2
1.6. ESTRUCTURAS DE PROGRAMA Y FUNCIONES 21

1.6.3. Ejemplo 1

Se tiene un circuito eléctrico formado por una fuente de voltaje variable en el tiempo
v(t) = 10cos(20t), una resistencia R = 5 y un diodo conectados en serie tal como se
muestra en la figura 3.4.
Para este circuito hacer
Escribir la función que modela el circuito,
Escribir el código para resolver el problema y
Graficar la corriente como función del tiempo

Figura 1.1: Circuito serie con diodo

La función para modelar el diodo en el circuito es:

function it = diodo(vt)
R = 5;

if vt > 0.7
vd = 0.7;
else
vd = vt;
end

it = (vt - vd)/R;
end
Dado que tenemos un operador condicional, no es posible hacer la sectorización. Es este caso
utilizamos la función arrayfun de matlab tal como se muestra en el siguiente código.
22 INTRODUCCIÓN

t = [0:0.01:2]; % tiempo
vt = 10*cos(20*t); % Voltaje variante en el tiempo
it = arrayfun(@diodo, vt); % Calculo de la corriente

plot(t, it, t, vt); % Corriente Calculada

xlabel(’t (seg.)’);
ylabel(’v(t) y i(t)’);
title (’Solución de un circuito con Diodo’);
La Figura 1.2 muestra la solución encontrada en función del tiempo

Solución de un circuito con Diodo

10

2
v(t) y i(t)

−2

−4

−6

−8

−10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
t (seg.)

Figura 1.2: Solución gráfica del circuito con Diodo

1.6.4. Ejemplo 2

Determinar el cruce por cero de la función f(x) = x – cos(x), utilizando el método de

Newton Raphson. Este algoritmo iterativo se resuelve haciendo

f (xk )
xk+1 = xk −
f 0 (xk )

La implementación en MATLAB es:

function z = cero(x0)

for k = 1: 100
1.6. ESTRUCTURAS DE PROGRAMA Y FUNCIONES 23

xnva = x0 - f(x0)/df(x0);
x0 = xnva;
fprintf(’Iteracion \%d f(\%f) = \%f\n’, k, x0, f(x0));
if(abs(f(x0)) < 0.00001) break; end;
end;

z =x0;

function y = f(x)
y = x-cos(x);

function dy = df(x)
dy = 1+sin(x);

1.6.5. Ejemplo 3

Determinar el cruce por cero de la función f(x) = x – cos(x), utilizando el método de

Bisecciones.

La implementación en MATLAB es:

function b = biseccion(x0,x1)
n=x0:0.1:x1;
plot(n,f(n))
xlabel(’eje x’);
ylabel(’eje y’);
hold on;
fmin = min(f(n));
fmax = max(f(n));

for S = 1:100
mau=((x1+x0)/2);
a=f(x0);
b=f(x1);
fmau=f(mau);
fprintf(’Iteracion %d f(mau)(%f) = %f\n’, S, mau, f(mau));

plot(x0, fmin:0.01:fmax, x1, fmin:0.01:fmax);

pause;
24 INTRODUCCIÓN

if a*fmau > 0 x0=mau;

elseif a*fmau==0 x0=mau;
elseif a*fmau==0 x1=mau;
elseif a*fmau < 0 x1=mau;

if (abs(fmau) < 0.0001), break, end;

end;
end;
disp([’La solucion esta en:’])
disp([mau]);

function y = f(x)
y = x-cos(x);
Solución de Sistemas Lineales de Ecuacio-
nes

2.1. Suma, resta, multiplicación y división de matrices

Una matriz se define como un arreglo bidimensional de datos, que tiene n renglones y m
columnas. Un elemento del arreglo puede ser identificado con aij

 
a1,1 a1,2 a1,3 ... a1,N

 a2,1 a2,2 a2,3 ... a2,N 

A=
 a3,1 a3,2 a3,3 ... a3,N 

 
aM,1 aM,2 aM,2 ... aM,N

Algunas de las operaciones básicas que pueden realizarse con matrices son suma, resta
y multiplicación. La división de matrices como tal no existe y en su lugar se calcula la
inversa.

2.1.1. Suma de matrices

Para que la sumar las matrices AN ×M y BR×S , se requiere que las matrices tengan el
mismo número de renglones N = R y de columnas M = S . Si queremos encontrar la suma
C = A + B, cada elemento de la matriz C lo calculamos de la siguiente forma:

Cn,m = An,m + Bn,m

para todos lo n, m en la matriz C

La implementación en Matlab es:

25
26 SIST. LINEALES

function C = suma(A, B)

[N, M] = size(A);
[R, S] = size(B);

if N~=R || M~=S
fprintf(’Las dimensiones tienen que ser iguales\n’);
return;
end;

C = zeros(N, M);

for n=1:N
for m = 1:M
C(n,m) = A(n,m)+B(n,m);
end;
end
end

2.1.2. Resta de matrices

En este caso, se deben cumplir las mismas propiedades que la rsuma de matrices y el calculo
de los elemento de la matriz C se calculan como:

Cn,m = An,m − Bn,m

para todos lo n, m en la matriz C

La implementación en Matlab es:
function C = resta(A, B)

[N, M] = size(A);
[R, S] = size(B);

if N~=R || M~=S
fprintf(’Las dimensiones tienen que ser iguales\n’);
return;
end;

C = zeros(N, M);
2.1. SUMA, RESTA, MULTIPLICACIÓN Y DIVISIÓN DE MATRICES 27

for n=1:N
for m = 1:M
C(n,m) = A(n,m)-B(n,m);
end;
end
end

2.1.3. Multiplicación de matrices

Para realizar el producto C = A ∗ B tenemos que considerar que el producto existe si

1.- El número de columnas de AN ×M es igual al número de renglones de BR×S , esto es
M = R. y
2.- Las dimensiones de la matriz resultante tendrá el mismo numero de renglones que la
matriz A y el mismo número de columnas que la matriz B. La matriz resultante C será de
tamaño N × S
3.- El cálculo de los elementos de la matriz C se lleva a cabo haciendo :

M
X
Cn,m = an,k ∗ bk,m
k=1

para todos lo n, m en la matriz C

La implementación en Matlab es:
function C = multiplica(A, B)
[N,M] = size(A);
[R,S] = size(B);

if M~=R
fprintf(’El numero de renglones no es igual al numero de columnas\n’);
return;
end;

C = zeros(N, S);

for n=1:N
for m= 1:S
suma = 0;
28 SIST. LINEALES

for k = 1:M
suma = suma + A(n,k)*B(k,m);
end;
C(n,m) = suma;
end;
end;
end

2.2. Factorización triangular LU

Supongamos que la matriz de coeficientes A de un sistema lineal Ax = b admite una

factorización triangular como la siguiente.

     
a1,1 a1,2 a1,3 1 0 0 u1,1 u1,2 u1,3
 a2,1 a2,2 a2,3  =  l2,1 1 0  ∗  0 u2,3 u2,3 
a3,1 a3,2 a3,3 l3,1 l3,2 1 0 0 u3,3

Entonces la solución la podemos plantear como

[LU ]x = b

El cual se puede solucionar resolviendo primero Ly = b y luego el sistema U x = y. Estos

sistemas pueden resolverse de manera muy eficiente utilizando sustitución hacia adelante
y sustitución hacia atrás.

2.2.1. Sustitución hacia adelante

Comencemos por dar un sistema triangular inferior

    
1 0 0 0 y1 b1
 l2,1 1 0 0   y2   b2 
  = 
 l3,1 l3,2 1 0   y3   b3 
l4,1 l4,2 l4,3 1 y4 b4

Si despejamos comenzando con la ecuación I a y1 y ası́ sucesivamente para cada una de

las ecuaciones tenemos:
2.2. FACTORIZACIÓN TRIANGULAR LU 29

y1 = b1
y2 = b2 − l2,1 y1
y3 = b3 − l3,1 y1 − l3,2 y2
y4 = b4 − l4,1 y1 − l4,2 y2 − l4,3 y3

La formula para calcular la sustitución hacia adelante es:

n<k
X
yk = bk − lk,n xn
n=1

Para todos los valores k = 1, 2, 3, · · · , N

La implementación en matlab es:
function y = sustitucion_hacia_adelante(l,b)
N = length(l(:,1));
y = zeros(N,1);
for k = 1: N
suma = 0;
for n = 1: k-1
suma = suma + (l(k,n)*y(n));
end;
y(k) = (b(k) - suma);
end;
end

2.2.2. Sustitución hacia atrás

Consideremos ahora un sistema triangular superior

    
u1,1 u1,2 u1,3 u1,4 x1 y1
 0 u2,2 u2,3 u2,4   x2   y2 
  = 
 0 0 u3,3 u3,4   x3   y3 
0 0 0 u4,4 x4 y4

Para este caso comenzamos despejando x4 de la última ecuación y sustituimos y despejamos

en las anteriores tal como:
30 SIST. LINEALES

y4
x4 =
u4,4
y3 − u3,4 x4
x3 =
u3,3
y2 − u2,3 x3 − u2,4 x4
x2 =
u2,2
y1 − u1,2 x2 − u1,3 x3 − u1,4 x4
x1 =
u1,1

La formula para calcular la sustitución hacia atrás es:

PN
yk − n=k+1 uk,n xn
xk =
uk,k

La implementación en Java es:

function x = sustitucion_hacia_atras(u, y)
N = length(y);
x = zeros(N,1);
for k = N: -1: 1
suma = 0;
for n = k+1: N
suma = suma +(u(k,n)*x(n));
end;
x(k) = (y(k)-suma)/u(k,k);
end;
end

2.2.3. Ejemplo

Resolver el sistema de ecuaciones

x1 + 2x2 + 4x3 + x4 = 21
2x1 + 8x2 + 6x3 + 4x4 = 52
3x1 + 10x2 + 8x3 + 8x4 = 79
4x1 + 12x2 + 10x3 + 6x4 = 82
2.2. FACTORIZACIÓN TRIANGULAR LU 31

Supongamos que tenemos la descomposición triangular siguiente

     
1 2 4 1 1 0 0 0 1 2 4 1
 2 8 6 4   2 1 0 0   0 4 –2 2 
 = × 
 3 10 8 8   3 1 1 0   0 0 –2 3 
4 12 10 6 4 1 2 1 0 0 0 –6

Primero resolvemos el sistema de ecuaciones

y1 = 21
2y1 + y2 = 52
3y1 + y2 + y3 = 79
4y1 + y2 + 2y3 + y4 = 82

El cual corresponde a un sistema triangular inferior

y1 = 21
y2 = 10
y3 = 6
y4 = −24

Posteriormente solucionamos U x = y

x1 + 2x2 + 4x3 + x4 = 21
4x2 − 2x3 + 2x4 = 10
−2x3 + 3x4 = 6
−6x4 = −24

Utilizando sustitución hacia atrás tenemos

X = [1, 2, 3, 4]T
32 SIST. LINEALES

2.2.4. Factorización de Crout

Consideremos un sistema más grande, por ejemplo de 4 ecuaciones con cuatro incógni-
tas.

     
a1,1 a1,2 a1,3 a1,4 1 0 0 0 u1,1 u1,2 u1,3 u1,4
 a2,1 a2,2 a2,3 a2,4   l2,1 1 0 0   0 u2,2 u2,3 u2,4 
 = × 
 a3,1 a3,2 a3,3 a3,4   l3,1 l3,2 1 0   0 0 u3,3 u3,4 
a4,1 a4,2 a4,3 a4,4 l4,1 l4,2 l4,3 1 0 0 0 u4,4

Si hacemos las multiplicaciones indicadas y dejamos los valores en la misma matriz tene-
mos:
Para el primer renglón de la matriz (n = 1)
arriba de la diagonal

m=1 a1,1 = a1,1

m=2 a1,2 = a1,2
m=3 a1,3 = a1,3
m=4 a1,4 = a1,4

Para el segundo renglón (n = 2)

bajo la diagonal

m = 1 a2,1 = a2,1 /a1,1

arriba de la diagonal

m = 2 a2,2 = a2,2 − a2,1 a1,2

m = 3 a2,3 = a2,3 − a2,1 a1,3
m = 4 a2,4 = a2,4 − a2,1 a1,4

Tercer renglón (n = 3)
bajo la diagonal

m=1 a3,1 = a3,1 /a1,1

m = 2 a3,2 = (a3,2 − a3,1 a1,2 )/a2,2
2.2. FACTORIZACIÓN TRIANGULAR LU 33

arriba de la diagonal

m = 3 a3,3 = a3,3 − a3,1 a1,3 − a3,2 a2,3

m = 4 a3,4 = a3,4 − a3,1 a1,4 − a3,2 a2,4

Y finalmente para el cuarto renglón (n = 4)

bajo la diagonal

m=1 a4,1 = a4,1 /a1,1

m=2 a4,2 = (a4,2 − a4,1 a1,2 )/a2,2
m = 3 a4,3 = (a4,3 − a4,1 a1,3 − a4,2 a2,3 )/a3,3

arriba de la diagonal

m = 4 a4,4 = a4,4 − a4,1 a1,4 − a4,2 a2,4 − a4,3 a3,4

Podemos ver, que cualquier elemento bajo la diagonal se calcula como:

Pm−1
an,m − k=1 an,k ak,m
an,m =
am,m
n = 1, · · · , N m = 1, · · · , n − 1

Para los términos arriba de la diagonal

n−1
X
an,m = an,m − ak,m
k=1

n = 1, · · · , N m = n, ..., N.

La implementación en Matlab es:

function a = Factorizacion(A)
N = length(A);
a = A;

for n = 1:N
for m = 1: n-1
suma = 0;
34 SIST. LINEALES

for k = 1:m-1
suma = suma + a(n,k) * a(k,m);
end;

a(n,m) = (a(n,m) - suma)/ a(m,m);

end;

for m = n:N
suma = 0;
for k = 1: n-1
suma = suma + a(n,k) * a(k,m);
end;
a(n,m) = a(n,m) - suma;
end;
end;
end

2.2.5. Ejemplo

Hacer la factorización LU de la siguiente matriz, utilizando eliminación Gaussiana. Corro-

borar utilizando Factorización de Crout.

 
1 2 4 1
 2 8 6 4 
 
 3 10 8 8 
4 12 10 6

Paso 1.

 
1 2 4 1
 2
 4 −2 2 

 3 4 −4 5 
4 4 −6 2

Paso 2.
2.2. FACTORIZACIÓN TRIANGULAR LU 35

 
1 2 4 1
 2
 4 −2 2 

 3 1 −2 3 
4 1 −4 0

Paso 3.

 
1 2 4 1
 2
 4 −2 2 

 3 1 −2 3 
4 1 −2 6

La ejecución en Matlab queda

>> Factorizacion( [1 2 4 1; 2 8 6 4;3 10 8 8;4 12 10 6])

ans =

1 2 4 1
2 4 -2 2
3 1 -2 3
4 1 2 -6

2.2.6. Ejemplo

Dado el sistema lineal de ecuaciones, calcular: a) La factorización de Crout y b) Calcular

la solución utilizando sustitución hacia adelante y atrás

    
10 −1 2 1 x1 1
 −1 15 −3 1   x2
    2 
 = 
 2 −3 6 −3   x3   2 
1 1 −3 7 x4 1

a) Realizamos la factorización de LU o Crout

>> Af = Factorizacion([10 -1 2 1; -1 15 -3 1; 2 -3 6 -3; 1 1 -3 7])

Af =

10.0000 -1.0000 2.0000 1.0000

36 SIST. LINEALES

-0.1000 14.9000 -2.8000 1.1000

0.2000 -0.1879 5.0738 -2.9933
0.1000 0.0738 -0.5899 5.0529
b) Par la sustitución hacia adelante hacemos
>> y = sustitucion_hacia_adelante(Af, [1,2,2,1]’)

y =

1.0000
2.1000
2.1946
2.0397
para luego realizar la sustitución hacia atrás
>> x = sustitucion_hacia_atras(Af, y)

x =

-0.0508
0.2372
0.6707
0.4037

2.2.7. Cálculo de la Matriz inversa utilizando descomposición LU

Si consideramos un sistema lineal Ax = b la solución la podemos llevar a cabo utilizando

el concepto de matriz inversa tal que x = A−1 b, por lo que podemos consideran que son
métodos equivalentes. Consideremos un sistema de 4 ecuaciones con 4 incógnitas como el
siguiente:

Ax = b

    
a1,1 a1,2 a1,3 a1,4 x1 0
 a2,1 a2,2 a2,3 a2,4   x2   0 
  = 
 a3,1 a3,2 a3,3 a3,4   x3   1 
a4,1 a4,2 a4,3 a4,4 x4 0

Si tuviéramos la matriz inversa, podemos calcular los valores de x haciendo:

2.2. FACTORIZACIÓN TRIANGULAR LU 37

x = A−1 b

  0
a01,2 a01,3 a01,4 a01,3
    
x1 a1,1 0
 x2   a02,1 a02,2 a02,3 a02,4   0   a02,3 
 x3  =  a03,1
=
a03,2 a03,3 a03,4   1   a03,3 
    

x4 a04,1 a04,2 a04,3 a04,4 0 a04,3

donde a0n,m son los elementos de la matriz inversa A−1 . Note que poniendo un 1 en b3
obtuvimos la tercer columna de la matriz inversa, por lo tanto si solucionamos el sistema
de ecuaciones utilizando descomposición LU obtendrı́amos la tercer columna.

Podemos generalizar el método si definimos


(k) bn = 1 si n = k
b =
bn = 0 si n <> k

Entonces si resolvemos el sistema Ax = b(k) por cualquier método, calculamos la k-esima

colma de la matriz inversa y la solución será x = A−1
k . Para calcular todas las columnas de
la matriz inversa solamente es necesario calcular este procedimiento para k = 1, 2, · · · , N .
Para hacer eficiente el procedimiento solamente factorizamos las matriz una vez y utilizando
sustitución hacia adelante y atrás resolvemos para cada uno de los valores de k. El código
para hacer esto es:

function Ainv = inversa(A)

N = length(A);

Af = Factorizacion(A);

Ainv = zeros(N,N);

for k = 1: N
b = zeros(N,1);
b(k) = 1;
y = sustitucion_hacia_adelante(Af, b);
x = sustitucion_hacia_atras(Af, y);
Ainv(:,k) = x;
end;
end
38 SIST. LINEALES

2.2.8. Ejemplo

Calcular la inversa de la matriz

 
10 −1 2 1
 −1 15 −3 1 
A=
 2 −3

6 −3 
1 1 −3 7

Para comenzar hacemos la factorización de la matriz utilizando el metido de Crout o

LU

 
10 −1 2 1
 −1 15 −3 1 
A=
 2 −3

6 −3 
1 1 −3 7

El resultado es

  
1.0000 0.0000 0.0000 0.0000 10.0000 −1.0000 2.0000 1.0000
 −0.1000 1.0000 0.0000 0.0000   0.0000 14.9000 −2.8000
  1.1000 
LU = 
 0.2000 −0.1879

1.0000 0.0000   0.0000 0.0000 5.0738 −2.9933 
0.1000 0.0738 −0.5899 1.0000 0.0000 0.0000 0.0000 5.0529

Para calcular la primer columna de la matriz inversa hacemos

  (1) 
y1
  
10.0000 −1.0000 2.0000 1.0000 1
 0.0000 14.9000 −2.8000 (1)
1.1000  y2   0 
   
  = 
 0.0000 0.0000 5.0738 −2.9933   y3(1) 0


0.0000 0.0000 0.0000 5.0529 (1)
y4 0

La solución es y (1) = [1.0000, 0.1000, −0.1812, −0.2143]T

  (1) 
x1
  
1.0000 0.0000 0.0000 0.0000 1.0000
 −0.1000 (1)
1.0000 0.0000 0.0000  x2   0.1000 
  
  = 
 0.2000 −0.1879 1.0000 0.0000   x(1) −0.1812 

3

0.1000 0.0738 −0.5899 1.0000 (1)
x4 −0.2143

La primer columna de la matriz inversa es x(1) = [0.1162, −0.0016, −0.0607, −0.0424]T

2.2. FACTORIZACIÓN TRIANGULAR LU 39

Para la segunda columna de la matriz inversa hacemos

  (2) 
y1
  
10.0000 −1.0000 2.0000 1.0000 0
 0.0000 14.9000 −2.8000 (2)
1.1000  y2   1 
   
  = 
 0.0000 0.0000 5.0738 −2.9933   y3(2) 0


0.0000 0.0000 0.0000 5.0529 y4
(2) 0

La solución es y (2) = [0, 1.0000, 0.1879, 0.0370]T

  (2) 
x1
  
1.0000 0.0000 0.0000 0.0000 0
 −0.1000  (2)
1.0000 0.0000 0.0000   x2   1.0000 
  
 = 
 0.2000 −0.1879 1.0000 0.0000   x(2) 0.1879 

3

0.1000 0.0738 −0.5899 1.0000 (2)
x4 0.0370

La segunda columna de la matriz inversa es x(2) = [−0.0016, 0.0743, 0.0414, 0.0073]T

Para la tercercolumna de la matriz inversa hacemos

  (3) 
y1
  
10.0000 −1.0000 2.0000 1.0000 0
 0.0000 14.9000 −2.8000 (3)
1.1000  y2   0 
   
  = 
 0.0000 0.0000 5.0738 −2.9933   y3(3) 1


0.0000 0.0000 0.0000 5.0529 y4
(3) 0

La solución es y (3) = [0, 0, 1.0000, 0.5899]T

  (3) 
x1
  
1.0000 0.0000 0.0000 0.0000 0
 −0.1000  (3)
1.0000 0.0000 0.0000   x2 0 
  
=
   
 0.2000 −0.1879 1.0000 0.0000   x(3) 1.0000 

3

0.1000 0.0738 −0.5899 1.0000 (3)
x4 0.5899

La tercer columna de la matriz inversa es x(3) = [−0.0607, 0.0414, 0.2660, 0.1168]T

Para la última columna de la matriz inversa hacemos

  (4) 
y1
  
10.0000 −1.0000 2.0000 1.0000 0
 0.0000 14.9000 −2.8000 (4)
1.1000  y2   0 
   
  = 
 0.0000 0.0000 5.0738 −2.9933   y3(4) 0


0.0000 0.0000 0.0000 5.0529 y4
(4) 1
40 SIST. LINEALES

La solución es y (4) = [0, 0, 0, 1]T

  (4) 
x1
 
1.0000 0.0000 0.0000 0.0000 0
 −0.1000 (4)
1.0000 0.0000 0.0000  x2   0 
   
  = 
 0.2000 −0.1879 1.0000 0.0000   x(4) 0

3

0.1000 0.0738 −0.5899 1.0000 (4)
x4 1

La última columna de la matriz inversa es x(4) = [−0.0424, 0.0073, 0.1168, 0.1979]T

Finalmente la matriz inversa es la concatenación de la soluciones

 
0.1162 −0.0016 −0.0607 −0.0424
 −0.0016 0.0743 0.0414 0.0073 
A−1 = [x(1) , x(2) , x(3) , x(4) ] = 
 −0.0607

0.0414 0.2660 0.1168 
−0.0424 0.0073 0.1168 0.1979

2.3. Eliminación Gaussiana

Consideremos que tenemos un sistema lineal Ax = b, donde la matriz A no tiene las

condiciones de ser triangular superior.

    
a1,1 a1,2 a1,3 x1 b1
 a2,1 a2,2 a2,3   x2  =  b2 
a3,1 a3,2 a3,3 x3 b3

Comenzaremos por despejar x1 de la ecuación I

b1 − a1,2 x2 − a1,3 x3
x1 =
a1,1

y sustituimos en las ecuaciones II

 
b1 − a1,2 x2 − a1,3 x3
a2,1 + a2,2 x2 + a2,3 x3 = b2
a1,1

agrupando términos semejantes tenemos:

     
a2,1 ∗ a1,2 a1,3 b1
a2,2 − x2 + a2,3 − a2,1 x3 = b2 − a2,1
a1,1 a1,1 a1,1
2.3. ELIMINACIÓN GAUSSIANA 41

Ahora sustituimos en la ecuación III

 
b1 − a1,2 x2 − a1,3 x3
a3,1 + a3,2 x2 + a3,3 x3 = b3
a1,1
     
a3,1 a1,2 a3,1 a1,3 b1
a3,2 − x2 + a3,3 − x3 = b3 − a3,1
a1,1 a1,1 a1,1

Lo cual nos da un nuevo sistema simplificado dada por

    
a1,1 a1,2 a1,3 x1 b1
 0 a02,2 a02,3   x2  =  b02 
0 a03,2 a03,3 x3 b03

donde los valores de a0n,m los calculamos como:

ak,m
a0n,m = an,m − an,k
ak,k
bk
b0n = bn − an,k
ak,k

Si repetimos el procedimiento, ahora para x1 , tendremos un sistema dado como:

    
a1,1 a1,2 a1,3 x1 b1
 0 a02,2 a02,3   x2  =  b02 
0 0 a003,3 x3 b003

La solución del sistema Triangular superior lo podemos realizar utilizando el método de

sustitución hacia atrás. La implementación en Matlab queda
function x = Eliminacion_Gaussiana(A0, b0)

a = A0;
b = b0;

N = length(b);

for k =1:N-1
for n = k+1:N
b(n) = b(n) - a(n,k)*b(k)/a(k,k);
42 SIST. LINEALES

for m=N:-1:k
a(n,m) = a(n,m) - a(n,k)*a(k,m)/a(k,k);
end;
end;
end;
disp(a);
disp(b);

x = sustitucion_hacia_atras(a, b);
end
que corresponde a un sistema triangular superior que podemos solucionar utilizando susti-
tución hacia atrás.

Ejemplo 1

Dado el sistema lineal de ecuaciones, calcular el sistema Triangular Superior utilizando el

algoritmo de Eliminación Gaussiana

    
10 −1 2 1 x1 1
 −1 15 −3 1   x2
    2 
 = 
 2 −3 6 −3   x3   2 
1 1 −3 7 x4 1

Primer paso k = 1

−1
    
10 2 1 x1 1
    
149
− 14 11 21
    
 0   x2   
 10 5 10    10 
  =  
    
 0
 − 14
5
28
5 − 16
5
  x3  
  
9
5


    
11
0 10 − 16
5
69
10 x4 9
10

    
10.0000 −1.0000 2.0000 1.0000 x1 1.0000

 0 14.9000 −2.8000 1.1000   x2

  2.1000 
= 
 0 −2.8000 5.6000 −3.2000   x3   1.8000 
0 1.1000 −3.2000 6.9000 x4 0.9000

Segundo paso k = 2
2.3. ELIMINACIÓN GAUSSIANA 43

10 −1
    
2 1 x1 1
    
149
− 14 11 21
    
 0   x2   
 10 5 10    10 
  =  
 756
   
 0
 0 149 − 446
149
  x3  
  
327
149


    
446 1016 111
0 0 − 149 149 x4 149

   
10.0000 −1.0000 2.0000 1.0000 x1 1.0000

 0 14.9000 −2.8000 1.1000   x2

  2.1000 
= 
 0 0 5.0738 −2.9933   x3   2.1946 
0 0 −2.9933 6.8188 x4 0.7450

Tercer paso k = 3

10 −1
    
2 1 x1 1
    
149
− 14 11 21
    
 0   x2   
 10 5 10    10 
  =  
 756
   
 0
 0 149 − 446
149
  x3  
  
327
149


    
955 257
0 0 0 189 x4 126

   
10.0000 −1.0000 2.0000 1.0000 x1 1.0000

 0 14.9000 −2.8000 1.1000   x2

  2.1000 
= 
 0 0 5.0738 −2.9933   x3   2.1946 
0 0 0 5.0529 x4 2.0397

Sustitución hacia atras

Paso 1 k = 4

b4
x4 =
a4,4
2.0397
x4 = = 0.4037
5.0529
44 SIST. LINEALES

Paso 2 k = 3

b3 − a3,4 x4
x3 =
a3,3
2.1946 − (−2.9933)(0.4037)
x3 = = 0.6707
5.0738

Paso 3 k = 2

b2 − a2,3 x3 − a2,4 x4
x2 =
a2,2
2.1000 − (−2.8000)(0.6707) − (1.1000)(0.4037)
x2 = = 0.2372
14.9000

Paso 4 k = 1

b1 − a1,2 x2 − a1,3 x3 − a1,4 x4

x1 =
a1,1
1.0000 − (−1.0000)(0.2372) − (2.0000)(0.6707) − (1.0000)(0.4037)
x1 = = −0.0508
10.0000

La solución en Matlab es

>> Eliminacion_Gaussiana([10 -1 2 1; -1 15 -3 1; 2 -3 6 -3; 1 1 -3 7], [1,2,2,1]’)

10.0000 -1.0000 2.0000 1.0000

0 14.9000 -2.8000 1.1000
0 0 5.0738 -2.9933
0 0 0 5.0529

1.0000 2.1000 2.1946 2.0397

ans =

-0.0508
0.2372
0.6707
0.4037
2.3. ELIMINACIÓN GAUSSIANA 45

Ejemplo 2

Calcular el sistema triangular superior utilizando eliminación Gaussiana.

5x + 2y + 1z = 3
2x + 3y − 3z = −10
1x − 3y + 2z = 4

Primer paso

5x + 2y + 1z = 3
0 + (11/5)y − (17/5)z = −(56/5)
0 − (17/5)y + (9/5)z = (17/5)

Segundo paso

5x + 2y + 1z = 3
0x + (11/5)y − (17/5)z = −(56/5)
0x − 0y − (38/11)z = −(153/11)

La solución en Matlab es
>> Eliminacion_Gaussiana([5 2 1; 2 3 -3; 1 -3 2], [3, -10, 4])
5.0000 2.0000 1.0000
0 2.2000 -3.4000
0 0 -3.4545

3.0000
-11.2000
-13.9091

ans =

-0.6579
1.1316
4.0263
46 SIST. LINEALES

Ejemplo 3

Resolver el sistema de ecuaciones

    
3 −1 −1 x 0
 −1 1 0   y = 1 
 
−1 0 1 z 1

Aplicando Eliminación Gaussiana nos queda.

    
3 −1 −1 x 0
 0 2/3 −1/3   y  =  1 
0 0 1/2 z 3/2

La solución en Matlab es

>> Eliminacion_Gaussiana([3 -1 -1; -1 1 0; -1 0 1], [0 1 1])

3.0000 -1.0000 -1.0000
0 0.6667 -0.3333
0 0 0.5000

0
1.0000
1.5000

ans =

2.0000
3.0000
3.0000

Problema con la Eliminación Gaussiana

Un ejemplo de sistema donde es necesario hacer un cambio de renglón por renglón para
que tenga solución es el siguiente sistema.
2.3. ELIMINACIÓN GAUSSIANA 47

 
1 2 6
 4 8 −1 
−2 3 5

Aplicando el primer paso de la eliminación gaussiana tenemos:

 
1 2 6
 0 0 −25 
0 7 17

note, que aparece un cero en el elemento 22, lo cual nos da un sistema sin solución. Per-
mutando los renglones 2 y 3 el sistema tiene solución.

 
1 2 6
 −2 3 5 
4 8 −1

La solución en Matlab tenemos

>> Eliminacion_Gaussiana([1 2 6; 4 8 -1; -2 3 5], [1 1 1])
1 2 6
0 0 -25
0 NaN Inf

1
-3
Inf

ans =

NaN
NaN
NaN
Haciendo el cambio de renglones tenemos
>> Eliminacion_Gaussiana([1 2 6; -2 3 5; 4 8 -1], [1 1 1])
1 2 6
0 7 17
0 0 -25
48 SIST. LINEALES

1
3
-3

ans =

0.0057
0.1371
0.1200

2.4. Gauss-Jordan

El método de Gauss-Jordan es una variación de la eliminación de Gauss. La principal

diferencia consiste em que cuando una incógnita se elimina en el método de Gauss-Jordan,
esta es eliminada de todas las ecuaciones en lugar de hacerlo en la subsecuentes. Además
todos los renglones se normalizan al dividirlos entre su elemento privote. De esta forma,
el paso de eliminación genera una matriz identidad en vez de una matriz triangular. En
consecuencia no es necesario usar la sustitución hacia atrás para obtener la solución.

Comenzaremos planteando un sistema de 3 ecuaciones con 3 incógnitas dado por

    
a1,1 a1,2 a1,3 x1 b1
 a2,1 a2,2 a2,3   x2  =  b2 
a3,1 a3,2 a3,3 x3 b3

Este sistema de ecuaciones lo podemos escribir como:

Ax = Ib

donde I es la matriz identidad

     
a1,1 a1,2 a1,3 x1 1 0 0 b1
 a2,1 a2,2 a2,3   x2  =  0 1 0   b2 
a3,1 a3,2 a3,3 x3 0 0 1 b3
2.4. GAUSS-JORDAN 49

Lo cual da lugar a la matriz aumentada que hemos utilizado en el procedimiento de Gauss-

Jordan

 
a1,1 a1,2 a1,3 1 0 0 b1
 a2,1 a2,2 a2,3 0 1 0 b2  (2.1)
a3,1 a3,2 a3,3 0 0 1 b3

Primer paso
Vamos a despejar la variable x1 del sistema y la sustituimos en las otras dos ecuacio-
nes

b1 − a1,2 x2 − a1,3 x3
= x1 (2.2)
a1,1

Sustituyendo en la ecuación 2

 
b1 − a1,2 x2 − a1,3 x3
a2,1 + a2,2 x2 + a2,3 x3 = b2
a1,1
     
a2,1 a2,1 a1,2 a2,1 a1,3
0x1 + b1 + a2,2 − x2 + a2,3 x3 = b2 (2.3)
a1,1 a1,1 a1,1

Sustituyendo en la ecuación 3

 
b1 − a1,2 x2 − a1,3 x3
a3,1 + a3,2 x2 + a3,3 x3 = b3
a1,1

     
a3,1 a3,1 a1,2 a3,1 a1,3
0x1 + b1 + a3,2 − x2 + a3,3 − x3 = b3 (2.4)
a1,1 a1,1 a1,1

En forma matricial podemos escribir las ecuaciones 2.2, 2.3 y 2.4

 a1,2 a1,3 
1 a1,1 a1,1
   1 
0 0  b1 
    x 1 a1,1
a a1,2 a a1,3 a2,1
a2,2 − 2,1 a2,3 − 2,1
 
 0 a1,1 a1,1
  x2  =  −
 a1,1 1 0 
  b2 
a
     
0
a3,1 a1,2
a3,2 − a1,1
a3,1 a1,3
a3,3 − a1,1 x3 − a3,1 0 1 b3
1,1

La matriz aumentada dada por (2.1 ) tenemos:

50 SIST. LINEALES

 a1,2 a1,3 1 1 
1 a1,1 a1,1 a1,1 0 0 a1,1 b1
   
a a1,2 a a1,3 a a
a2,2 − 2,1 a2,3 − 2,1 − a2,1 b2 − a2,1
 
 0 a1,1 a1,1 1 0 b1 
 1,1 1,1 
   
a3,1 a1,2 a3,1 a1,3 a a
0 a3,2 − a1,1 a3,3 − a1,1 − a3,1
1,1
0 1 b3 − a3,1
1,1
b1

Segundo Paso
Para un sistema equivalente

1 a01,2 a01,3
    0  0 
x1 c1,1 0 0 b1
 0 a02,2 a02,3   x2  =  c02,1 1 0   b02  (2.5)
0 a03,2 a03,3 x3 c03,1 0 1 b03

Despejamos x2 de la ecuación 2

c02,1 b01 + b02 − a02,3 x3

x2 =
a02,2

Sustituyendo en la ecuación 1

!
c02,1 b01 + b02 − a02,3 x3
x1 + a01,2 + a01,3 = c01,1 b01
a02,2

! !
0
a01,2 a02,3 0 0
a01,2 c02,1 a02,1 0
x1 + 0x2 + a1,3 − x3 = c1,1 b1 − − 0 b2
a02,2 a02,2 a2,2

Sustituyendo en la ecuación 3

! !
a03,2 a02,3 a03,2 c02,1 a03,2 0
0x1 + 0x2 + a03,3 − x3 = c03,1 − b01 − b 2 + b0 3
a02,2 a02,2 a02,2

En forma matricial

 a01,2 a02,3   a01,2 c02,1 a01,2 

1 0 a01,3 − a02,2
c01,1 − a02,2 a02,2
0
 b01
   
 x1
a02,3 c02,1
 
1  0 
 0 1  x2  =  0  b2
  
a02,2 a02,2 a02,2
a0 a02,3 x3 a0 c02,1 a0 b03
   
0 0 a03,3 − 3,2 a0
c03,1 − 3,2 a02,2
− a03,2 1
2,2 2,2
2.4. GAUSS-JORDAN 51

Tercer Paso
Dado el sistema equivalente

1 0 a001,3
  00
c1,1 c001,2 0
    00 
x1 b1
 0 1 a002,3   x2  =  c002,1 c002,2 0   b002 
0 0 a003,3 x3 c003,1 c003,2 1 b003

Despejamos x3 de la ecuación 3 y sustituimos en la ecuaciones 1 y 2

 a00 00
1,3 c3,2 a00 
c001,1 c001,2 − 00 − a1,3
00
a3,3
 b001
    
1 0 0 x1 3,3
a00 a00

  00 
 c002,1
 0 1 0   x2  =  00
c2,2 − a2,300 − a2,3
00  b2
3,3 3,3
0 0 1 x3 b003
 00
c00

c3,1 3,2 1
a00
3,3 a00
3,3 a00
3,3

Como resultado tenemos que la matrix c00 es la inversa de nuestro sistema y b00 la solución
del sistema de ecuaciones
Resumen
Dado el sistema

 
a1,1 a1,2 a1,3 · · · a1,M
 a2,1 a2,2 a2,3 · · · a2,M 
a3,1 a3,2 a3,3 · · · a3,M

Paso 1
Dividimos la primer ecuación entre a1,1

1  
a1,1 a1,1 a1,2 a1,3 ··· a1,M
 a2,1 a2,2 a2,3 ··· a2,M 
a3,1 a3,2 a3,3 ··· a3,M

Multiplicamos la ecuación 1 por a2,1 y la restamos a la ecuación 2

a1,2 a1,3 a1,M

1 ···
 
a2,1 a1,1 a1,1 a1,1
 a2,1 a2,2 a2,3 ··· a2,M 
a3,1 a3,2 a3,3 ··· a3,M
52 SIST. LINEALES

a1,2 a1,3 a1,M

···
 
1 a1,1 a1,1 a1,1 
     
a1,2 a1,3 a1,M
 0 a2,2 − a2,1 a2,3 − a2,1 ··· a2,M − a2,1

a1,1 a1,1 a1,1 
a3,1 a3,2 a3,3 ··· a3,M

Multiplicamos la ecuación 1 por a3,1 y la restamos a la ecuación 3

a1,2 a1,3 a1,M

···
 
a3,1 1 a1,1 a1,1 a1,1 
     
a1,2 a1,3 a1,M
 0 a2,2 − a2,1 a2,3 − a2,1 ··· a2,M − a2,1

a1,1 a1,1 a1,1 
a3,1 a3,2 a3,3 ··· a3,M

a1,2 a1,3 a1,M

···
 
1 a1,1 a1,1 a1,1 
     
a
 0 a2,2 − a2,1 aa1,2 a
a2,3 − a2,1 a1,3 ··· − a2,1 a1,M

a2,M 
  1,1   1,1   1,1  
a a a
0 a3,2 − a3,1 a1,2
1,1
a3,3 − a3,1 a1,3
1,1
··· a3,M − a3,1 a1,M
1,1

Paso 2
Dividimos la segunda ecuación entre a2,2

 
a1,1 a1,2 a1,3 ··· a1,M
1
a2,2
 a2,1 a2,2 a2,3 ··· a2,M 
a3,1 a3,2 a3,3 ··· a3,M

a1,1 a1,2 a1,3 ··· a1,M

 
a2,M
 aa2,1 1
a2,3
a2,2 ··· a2,2

2,2
a3,1 a3,2 a3,3 ··· a3,M

Multiplicamos la ecuación 2 por a1,2 y la restamos a la primer ecuación

a1,1 a1,2 a1,3 ··· a1,M

 
a2,M
a1,2  aa2,1 1
a2,3
a2,2 ··· a2,2

2,2
a3,1 a3,2 a3,3 ··· a3,M

      
a2,1 a2,3 a2,M

a1,1 − a1,2 a2,2 0 a1,3 − a1,2 a2,2 ··· a1,M − a1,2 a2,2
 a2,1 a2,3 a2,M 
 a2,2 1 a2,2 ··· a2,2

a3,1 a3,2 a3,3 ··· a3,M
2.4. GAUSS-JORDAN 53

Multiplicamos la ecuación 2 por a3,2 y se la restamos a la ecuación 3

      
a2,1 a2,3 a2,M

a1,1 − a1,2 a2,2 0 a1,3 − a1,2 a2,2 ··· a1,M − a1,2 a2,2
a3,2 
 a2,1 a2,3 a2,M 
a2,2 1 a2,2 ··· a2,2

a3,1 a3,2 a3,3 ··· a3,M
      
a2,1 a2,3 a2,M

a1,1 − a1,2 a2,2 0 a1,3 − a1,2 a2,2 ··· a1,M − a1,2 a2,2
 a2,1 a2,3 a2,M 

 a2,2 1 a2,2 ··· a2,2 

     
a2,1 a2,3 a2,M
a3,1 − a3,2 a2,2 0 a3,3 − a3,2 a2,2 ··· a3,1 − a3,2 a2,2

Paso 3
Dividimos la tercera ecuación entre a3,3

 
a1,1 a1,2 a1,3 ··· a1,M
 a2,1 a2,2 a2,3 ··· a2,M 
1
a3,3 a3,1 a3,2 a3,3 ··· a3,M

a1,1 a1,2 a1,3 ··· a1,M

 
 a2,1 a2,2 a2,3 ··· a2,M 
a3,1 a3,2 a3,M
a3,3 a3,3 1 ··· a3,3

Multiplicamos la ecuación 3 por a1,3 y la restamos a la primer ecuación

a1,1 a1,2 a1,3 ··· a1,M

 
 a2,1 a2,2 a2,3 ··· a2,M 
a3,1 a3,2 a3,M
a1,3 a3,3 a3,3 1 ··· a3,3

      
a3,1 a3,2 a3,M

a1,1 − a1,3 a3,3 a1,2 − a1,3 a3,3 0 ··· a1,M − a1,3 a3,3
 
 a2,1 a2,2 a2,3 ··· a2,M 
a3,1 a3,2 a3,M
a3,3 a3,3 1 ··· a3,3

Multiplicamos la ecuación 3 por a2,3 y la restamos a la segunda ecuación

      
a3,1 a3,2 a3,M

a1,1 − a1,3 a3,3 a1,2 − a1,3 a3,3 0 ··· a1,M − a1,3 a3,3
 
 a2,1 a2,2 a2,3 ··· a2,M 
a2,3 a3,1 a3,2 a3,M
a3,3 a3,3 1 ··· a3,3
54 SIST. LINEALES

      
a3,1 a3,2 a3,M

a1,1 − a1,3 a3,3 a1,2 − a1,3 0 ··· a1,M − a1,3
    a3,3   a3,3  
a3,1 a a
 a −
 2,1 a2,3 a3,3 a2,2 − a2,3 a3,2
3,3
0 ··· a2,M − a2,3 a3,M
3,3


a3,1 a3,2 a3,M
a3,3 a3,3 1 ··· a3,3

Note la similitud del método con el método de sumas y restas. No olvidar que este método
es equivalente al calculo parcial de la matriz inversa. El método se ilustra mejor con un
ejemplo.

2.4.1. Ejemplo 1

Utilice la técnica de Gauss-Jordan para resolver el siguiente sistema de ecuaciones.

30x1 − x2 − 2x3 = 78
x1 + 70x2 − 3x3 = −193
3x1 − 2x2 + 100x3 = 714

Primero exprese los coeficientes y el lado derecho como una matriz aumentada

 
30 −1 −2 1 0 0 78
 1 70 −3 0 1 0 −193 
3 −2 100 0 0 1 714

Primer iteración
Luego se normaliza el primer renglón, al dividirlo entre el elemento pivote, 30, se obtie-
ne

 1 1 1 13 
1 − 30 − 15 30 0 0 5
 
 
 1
 70 −3 0 1 0 −193 

 
3 −2 100 0 0 1 714

El término x1 se puede eliminar del segundo renglón restando el primer renglón multiplicado
por 1 del segundo renglón. En forma similar restando el primer renglón multiplicando por
3 eliminará el término x1 del tercer renglón
2.4. GAUSS-JORDAN 55

 1 1 1 13 
1 − 30 − 15 30 0 0 5
 
 0 2101 − 44 1
− 978
 
 30 15 − 30 1 0 5


 
0 − 19
10
501
5
1
− 10 0 1 3531
5

 
1. −0.0333 −0.0667 0.0333 0. 0. 2.6
 0. 70.0333 −2.9333 −0.0333 1. 0. −195.6 
0. −1.9 100.2 −0.1 0. 1. 706.2

Segunda Iteración
2101
En seguida, se normaliza el segundo renglón dividiéndolo entre 30

 1 1 1 13 
1 − 30 − 15 30 0 0 5
 
8 1 30
− 5868
 
 0
 1 − 191 − 2101 2101 0 2101


 
19 501 1 3531
0 − 10 5 − 10 0 1 5

 
1. −0.0333 −0.0667 0.0333 0. 0. 2.6
 0. 1. −0.0419 −0.0005 0.0143 0. −2.793 
0. −1.9 100.2 −0.1 0. 1. 706.2

Al reducir los términos x2 de las ecuaciones 1 y 2 tenemos

 13 70 1 5267 
1 0 − 191 2101 2101 0 2101
 
 0 1 − 8 1 30
− 5868
 
 191 − 2101 2101 0 2101


 
0 0 19123
191
212
− 2101 57
2101 1 1472577
2101

 
1. 0. −0.0681 0.0333 0.0005 0. 2.5069
 0. 1. −0.0419 −0.0005 0.0143 0. −2.793 
0. 0. 100.12 −0.1009 0.0271 1. 700.893

Tercer iteración
19123
El tercer renglón se normaliza entonces al dividirlo entre 191
56 SIST. LINEALES

 13 70 1 5267 
1 0 − 191 2101 2101 0 2101
 
 0 1 − 8 1 30
− 5868
 
 191 − 2101 2101 0 2101


 
212 57 191 1472577
0 0 1 − 210353 210353 19123 210353

 
1. 0. −0.0681 0.0333 0.0005 0. 2.5069
 0. 1. −0.0419 −0.0005 0.0143 0. −2.793 
0. 0. 1. −0.001 0.0003 0.01 7.0005

Por último, los términos x3 se pueden reducir de la primera y segunda ecuación para
obtener

 538 8 1 48274 
1 0 0 16181 16181 1471 16181
 
 0 1 0 − 109 3006 8
− 525828
 

 210353 210353 19123 210353 
 
212 57 191 1472577
0 0 1 − 210353 210353 19123 210353

 
1. 0. 0. 0.0332 0.0005 0.0007 2.9834
 0. 1. 0. −0.0005 0.0143 0.0004 −2.4997 
0. 0. 1. −0.001 0.0003 0.01 7.0005

2.4.2. Ejemplo 2

Dado el sistema lineal de ecuaciones calcular la solución utilizando el método de Gauss-

Jordan

    
1 1 2 x1 −1
 1 3 −6   x2  =  7 
2 −1 2 x3 0

La matriz aumentada para este sistema es:

 
1 1 2 −1
 1 3 −6 7 
2 −1 2 0

Primer iteración
2.4. GAUSS-JORDAN 57

 
1 1 2 −1
 0 2 −8 8 
0 −3 −2 2

Segunda iteración

 
1 0 6 −5
 0 1 −4 4 
0 0 −14 14

Tercer iteración

 
1 0 0 1
 0 1 0 0 
0 0 1 −1

2.4.3. Ejemplo 3

Dado el sistema lineal de ecuaciones, calcular la solución y la matriz inversa utilizando el

método de Gauss-Jordan

    
10 −1 2 1 x1 1
 −1 15 −3 1   x2   2 
  = 
 2 −3 6 −3   x3   2 
1 1 −3 7 x4 1

El sistema inicial aumentado es

 
10 −1 2 1 1 0 0 0 1
 −1 15 −3 1 0 1 0 0 2 
 
 2 −3 6 −3 0 0 1 0 2 
1 1 −3 7 0 0 0 1 1

Primer iteración
58 SIST. LINEALES

1 1 1 1 1
1 − 10
 
5 10 10 0 0 0 10
 
 0 149 − 14 11 1 21
 
 10 5 10 10 1 0 0 10


 
 
 0 − 14 28
− 16 − 51 0 1 0 9 
 5 5 5 5 
 
11
0 10 − 16
5
69
10
1
− 10 0 0 1 9
10

 
1.0000 −0.1000 0.2000 0.1000 0.1000 0 0 0 0.1000

 0 14.9000 −2.8000 1.1000 0.1000 1.0000 0 0 2.1000 

 0 −2.8000 5.6000 −3.2000 −0.2000 0 1.0000 0 1.8000 
0 1.1000 −3.2000 6.9000 −0.1000 0 0 1.0000 0.9000

Segunda iteración

 27 16 15 1 17 
1 0 149 149 149 149 0 0 149
 
 0 1 − 28 11 1 10 21
 
 149 149 149 149 0 0 149


 
 756

 0 0
 149 − 446
149
27
− 149 28
149 1 0 327
149


 
446 1016 16 11 111
0 0 − 149 149 − 149 − 149 0 1 149

 
1.0000 0 0.1812 0.1074 0.1007 0.0067 0 0 0.1141

 0 1.0000 −0.1879 0.0738 0.0067 0.0671 0 0 0.1409 

 0 0 5.0738 −2.9933 −0.1812 0.1879 1.0000 0 2.1946 
0 0 −2.9933 6.8188 −0.1074 −0.0738 0 1.0000 0.7450

Tercer iteración

3 3 1 1
0 − 28
 
1 0 0 14 28 0 28
 
 0 1 0 −1 2 1 2
 
 27 0 27 27 0 9


 
 
 0 0 1 − 223 1
− 28 1 149
0 109 
 378 27 756 252 
 
955 3 1 223 257
0 0 0 189 − 14 27 378 1 126
2.5. INVERSA DE SHIPLAY 59

 
1.0000 0 0 0.2143 0.1071 0.0000 −0.0357 0 0.0357

 0 1.0000 0 −0.0370 0 0.0741 0.0370 0 0.2222 

 0 0 1.0000 −0.5899 −0.0357 0.0370 0.1971 0 0.4325 
0 0 0 5.0529 −0.2143 0.0370 0.5899 1.0000 2.0397

Cuarta iteración

111 3 58 81 97
− 1910 − 955 − 1910 − 1910
 
1 0 0 0 955
 
 0 1 0 0 − 3 71 79 7 453
 

 1910 955 1910 955 1910 
 
 
 0 0 1 0 − 58 79 254 223 1281 
 955 1910 955 1910 1910 
 
81 7 223 189 771
0 0 0 1 − 1910 955 1910 955 1910

 
1.0000 0 0 0 0.1162 −0.0016 −0.0607 −0.0424 −0.0508

 0 1.0000 0 0 −0.0016 0.0743 0.0414 0.0073 0.2372 

 0 0 1.0000 0 −0.0607 0.0414 0.2660 0.1168 0.6707 
0 0 0 1.0000 −0.0424 0.0073 0.1168 0.1979 0.4037

2.5. Inversa de Shiplay

Consideremos el siguiente sistema de ecuaciones

    
a1,1 a1,2 a1,3 x1 b1
 a2,1 a2,2 a2,3   x2  =  b2 
a3,1 a3,2 a3,3 x3 b3

de la ecuación I despejamos el valor de la variable x1

b1 − a1,2 x2 − a1,3 x3
x1 =
a1,1

Reorganizando la ecuación tenemos

b1 a1,2 a1,3
− x2 − x3 = x1
a1,1 a1,1 a1,1
60 SIST. LINEALES

Sustituimos el valor de x1 en la ecuación II:

(b1 − a1,2 x2 − a1,3 x3 )a2,1

+ a2,2 x2 + a2,3 x3 = b2
a1,1
   
a2,1 a2,1 a1,2 a2,1 a1,3
b1 + a2,2 − x2 + a2,3 − x3 = b2
a1,1 a1,1 a1,1

y finalmente en III tenemos :

   
a3,1 a3,1 a1,2 a3,1 a1,3
b1 + a3,2 − x2 + a3,3 − x3 = b3
a1,1 a1,1 a1,1

En forma matricial estas ecuaciones lucen como:

    
1/a1,1 a1,2 /a1,1 a1,3 /a1,1 b1 x1
 a2,1 /a1,1 a2,2 − a2,1 a1,2 /a1,1 a2,3 − a2,1 a1,3 /a1,1   x2  =  b2 
a3,1 /a1,1 a3,2 − a3,1 a1,2 /a1,1 a3,3 − a3,1 a1,3 /a1,1 x3 b3

Lo que nos da un sistema equivalente como el siguiente.

a01,1 a01,2 a01,3

    
b1 x1
 a02,1 a02,2 a02,3   x2  =  b2 
a03,1 a03,2 a03,3 x3 b3

Si repetimos el procedimiento suponiendo que despejamos x2 y sustituimos en las otras dos

ecuaciones tendremos un sistema.

a01,1 − a01,2 a02,1 /a02,2 a01,2 /a02,2 a01,3 − a01,2 a02,3 /a02,2
    
b1 x1
 −a02,1 /a02,2 1/a02,2 −a02,3 /a02,2   b2  =  x2 
0 0 0 0 0 0 0 0
a3,1 − a3,2 a2,1 /a2,2 a3,2 /a2,2 a3,3 − a3,2 a2,3 /a2,20 0 x3 b3

Lo que nos da un sistema equivalente como el siguiente.

a001,1 a001,2 a001,3

    
b1 x1
 a002,1 a002,2 a002,3   b2  =  x2 
a003,1 a003,2 a003,3 x3 b3

Finalmente si despejamos x3 de la ecuación III y sustituimos en las otras tenemos;

2.5. INVERSA DE SHIPLAY 61

a001,1 − a001,3 a003,1 /a003,3 a001,2 − a001,3 a003,2 /a003,3 a001,3 /a003,3
    
b1 x1
 a002,1 − a002,3 a003,1 /a003,3 a002,2 − a002,3 a003,2 /a003,3 a002,3 /a003,3   b2  =  x2 
a003,1 /a003,3 a003,2 /a002,2 /a003,3 1/a003,3 b3 x3

El procedimiento para el calculo de la inversa por el método de Shiplay es:

Dado un valor de k = 1, 2, 3, ..., N al que llamamos pivote hacemos:
an,m = an,m –an,k ak,m /ak,k para todos los valores de n y m diferentes de k
an,k = an,k /ak,k para todos los n diferentes de k
ak,m = ak,m /ak,k para todos los m diferentes de k
ak,k = 1/ak,k
La implementación en Matlab queda
function a = Inversa_Shiplay(A)
a=A;
N=length(a);
for k=1:N
for n=1:N
for m=1:N
if ((m~=k)&(n~=k))
a(n,m)=a(n,m)-(a(n,k)*a(k,m))/a(k,k);
end;
end;
end;
for n=1:N
if (n~=k)
a(n,k)=a(n,k)/a(k,k);
a(k,n)=-(a(k,n)/a(k,k));
end;
end;
a(k,k)=1/a(k,k);
end;
end

2.5.1. Ejemplo:

Determinar la matriz inversa de

62 SIST. LINEALES

 
3 −1 −1
 −1 1 0 
−1 0 1

Primer paso:

 
1/3 1/3 1/3
 −1/3 2/3 −1/3 
−1/3 −1/3 2/3

Segundo paso:

 
1/2 1/2 1/2
 1/2 3/2 1/2 
−1/2 −1/2 1/2

Tercer paso:

 
1 1 1
 1 2 1 
1 1 2

La solución en Matlab es:

>> Inversa_Shiplay([ 3 -1 -1;-1 1 0 ; -1 0 1])

ans =

1.0000 1.0000 1.0000

1.0000 2.0000 1.0000
1.0000 1.0000 2.0000

2.6. Sistemas dispersos y estrategias para conservarla

Consideremos un sistema de ecuaciones Ax = b como la siguiente

2.6. SISTEMAS DISPERSOS Y ESTRATEGIAS PARA CONSERVARLA 63

    
10 −1 −1 −1 1 x1 1

 1 4 0 0 0 
 x2  
  2 

A=
 1 0 4 0 0   x3
 =
  3 

 1 0 0 4 0  x4
   4 
1 0 0 0 3 x5 5

Para ésta matriz el total de elementos es 25 de los cuales 13 son diferentes de cero y 12 son
igual cero. A este tipo de matrices o sistemas de ecuaciones se les conoce como dispersos o
ralos dado que tienen una cantidad de ceros igual o mayor que los diferentes de cero.
Para este tipo de sistema cuando calculamos la solución por lo métodos dados tene-
mos
Descomposición LU
>> Factorizacion(A)

ans =

10.0000 -1.0000 -1.0000 -1.0000 1.0000

0.1000 4.1000 0.1000 0.1000 -0.1000
0.1000 0.0244 4.0976 0.0976 -0.0976
0.1000 0.0244 0.0238 4.0952 -0.0952
0.1000 0.0244 0.0238 0.0233 2.9070
Eliminación Gaussiana
>> Eliminacion_Gaussiana(A,b)
10.0000 -1.0000 -1.0000 -1.0000 1.0000
0 4.1000 0.1000 0.1000 -0.1000
0 0 4.0976 0.0976 -0.0976
0 0 0 4.0952 -0.0952
0 0 0 0 2.9070
Inversa de Shiplay
>> Inversa_Shiplay(A)

ans =

0.0960 0.0240 0.0240 0.0240 -0.0320

-0.0240 0.2440 -0.0060 -0.0060 0.0080
-0.0240 -0.0060 0.2440 -0.0060 0.0080
-0.0240 -0.0060 -0.0060 0.2440 0.0080
64 SIST. LINEALES

-0.0320 -0.0080 -0.0080 -0.0080 0.3440

Note que en todos los casos la dispersidad se pierde y tenemos una matriz con números
diferentes de cero ahı́ donde existı́an ceros.

En este caso una opción es el manejo de Matrices Ralas del Matlab para no almacenar los
valores iguales a cero y utilizar métodos iterativo como los de la siguiente sección.

La forma de representar la matriz dispersa en Matlab es

A = sparse(5,5);
A(1,1) = 10;
A(1,2) = -1;
A(1,3) = -1;
A(1,4) = -1;
A(1,5) = 1;
A(2,1) = 1;
A(3,1) = 1;
A(4,1) = 1;
A(5,1) = 1;
A(2,2) = 4;
A(3,3) = 4;
A(4,4) = 4;
A(5,5) = 3;
disp(A);

2.7. Método iterativo de Jacobi y Gauss-Seidel

2.7.1. Método de Jacobi

Consideremos el siguiente sistema de ecuaciones.

    
a1,1 a1,2 a1,3 x1 b1
 a2,1 a2,2 a2,3   x2  =  b2 
a3,1 a3,2 a3,3 x3 b3

Vamos a representar cada una de las variables en términos de ellas mismas.

2.7. MÉTODO ITERATIVO DE JACOBI Y GAUSS-SEIDEL 65

b1 − a1,2 x2 − a1,3 x3
x1 =
a1,1
b2 − a2,1 x1 − a2,3 x3
x2 =
a2,2
b3 − a3,1 x1 − a3,2 x2
x3 =
a3,3

Lo cual nos sugiere el siguiente esquema iterativo de solución.

(t) (t)
(t+1) b1 − a1,2 x2 − a1,3 x3
x1 =
a1,1
(t) (t)
(t+1) b2 − a2,1 x1 − a2,3 x3
x2 =
a2,2
(t) (t)
(t+1) b3 − a3,1 x1 − a3,2 x2
x3 =
a3,3

En general podemos escribir como

PN (t)
(t+1) bk – l=1,l6=k ak,l xl
xk =
ak,k

La implementación en Matlab es
function y = Jacobi(A, x, b)

N = length(x);
y = zeros(N,1);

for iter=1:100000
for k = 1:N
suma =0;
for l= 1:N
if k ~= l
suma = suma + A(k,l)*x(l);
end;
66 SIST. LINEALES

end;
y(k) = (b(k) - suma)/A(k,k);
end;
if sqrt(norm(x-y)) < 1e-6
break;
else
x = y;
end;
end;

2.7.2. Algoritmo iterativo de Gauss-Seidel

El cambio que debemos hacer respecto al de Jacobi, es que las variables nuevas son uti-
lizadas una vez que se realiza el cálculo de ellas ası́, para un sistema de tres ecuaciones
tendremos:

(t) (t)
(t+1) b1 − a1,2 x2 − a1,3 x3
x1 =
a1,1
(t+1) (t)
(t+1) b2 − a2,1 x1 − a2,3 x3
x2 =
a2,2
(t+1) (t+1)
(t+1) b3 − a3,1 x1 − a3,2 x2
x3 =
a3,3

La implementación en Matlab es
function y = Gauss_Seidel(A, x, b)

N = length(x);
y = zeros(N,1);

for iter=1:100000
for k = 1:N
suma =0;
for l= 1:N
if k ~= l
suma = suma + A(k,l)*x(l);
end;
2.7. MÉTODO ITERATIVO DE JACOBI Y GAUSS-SEIDEL 67

end;
x(k) = (b(k) - suma)/A(k,k);
end;
if sqrt(norm(x-y)) < 1e-6
break;
else
y=x;
end;
end;

Ejemplo

Resolver el siguiente sistema de ecuaciones utilizando el método de Jacobi y comparar con

el método de Gauss-Seidel.

    
4 −1 1 x1 7
 4 −8 1   x2  =  −21 
−2 1 5 x3 15

Las primeras 20 iteraciones del algoritmo de Jacobi son :

k x1 x2 x3
1 1.0000 2.0000 2.0000
2 1.7500 3.3750 3.0000
3 1.8438 3.8750 3.0250
4 1.9625 3.9250 2.9625
5 1.9906 3.9766 3.0000
6 1.9941 3.9953 3.0009
7 1.9986 3.9972 2.9986
8 1.9996 3.9991 3.0000
9 1.9998 3.9998 3.0000
10 1.9999 3.9999 2.9999
11 2.0000 4.0000 3.0000
12 2.0000 4.0000 3.0000
Para correr hacer
>> Jacobi([4 -1 1; 4 -8 1; -2 1 5], [1,2,2]’, [7,-21, 15]’)

ans =
68 SIST. LINEALES

2.0000
4.0000
3.0000
La solución utilizando Gauss-Seidel es :
k x1 x2 x3
1 1.0000 2.0000 2.0000
2 1.7500 3.7500 2.9500
3 1.9500 3.9688 2.9863
4 1.9956 3.9961 2.9990
5 1.9993 3.9995 2.9998
6 1.9999 3.9999 3.0000
7 2.0000 4.0000 3.0000
8 2.0000 4.0000 3.0000
Note que Gauss-Seidel requiere de 7 iteraciones mientras Jacobi de 11, para convergir. Para
correr hacer
>> Gauss_Seidel([4 -1 1; 4 -8 1; -2 1 5], [1,2,2]’, [7,-21, 15]’)

ans =

2.0000
4.0000
3.0000

2.7.3. Ejemplo matrices dispersas

Resolver el sistema de ecuaciones utilizando matrices dispersas y los métodos de Jacobi y

Gauss-Seidel

    
10 −1 −1 −1 1 x1 1

 1 4 0 0 0   x2  
   2 

A=
 1 0 4 0 0   x3 
 
=
 3 

 1 0 0 4 0   x4   4 
1 0 0 0 3 x5 5
A = sparse(5,5);
A(1,1) = 10;
A(1,2) = -1;
A(1,3) = -1;
2.7. MÉTODO ITERATIVO DE JACOBI Y GAUSS-SEIDEL 69

A(1,4) = -1;
A(1,5) = 1;
A(2,1) = 1;
A(3,1) = 1;
A(4,1) = 1;
A(5,1) = 1;
A(2,2) = 4;
A(3,3) = 4;
A(4,4) = 4;
A(5,5) = 3;
disp(A);

size(A)

b = [1,2,3,4,5]’;

disp(b);
Jacobi(A,[0,0,0,0,0]’, b)
Gauss_Seidel(A,[0,0,0,0,0]’, b)
La solución para ambos métodos es:
ans =

0.1520
0.4620
0.7120
0.9620
1.6160

ans =

0.1520
0.4620
0.7120
0.9620
1.6160
70 SIST. LINEALES
Solución de Ecuaciones No lineales

3.1. Introducción

Una ecuación no lineal es aquella que tiene una forma diferente a f (x) = a0 + a1 x en cuyo
caso calcular la solución consiste en resolver despejando de x = −a0 /a1 . Pero el caso es
que queremos resolver un sistema de ecuaciones no lineales de la forma

f1 (x1 , x2 , · · · , xN ) = 0
f2 (x1 , x2 , · · · , xN ) = 0
.. .
. = ..
fn (x1 , x2 , · · · , xN ) = 0

donde N es el número de ecuaciones y se tiene el mismo números de funciones fi que

variables xi . En esta sección veremos dos de estos métodos que son la iteración de punto
fijo y el método de Newton-Raphson, los cuales se pueden utilizar para resolver sistemas
de ecuaciones no lineales

3.2. Métodos de punto fijo

Considerando un sistema no lineal de ecuaciones

f1 (x1 , x2 , · · · , xn ) = 0
f2 (x1 , x2 , · · · , xn ) = 0
.. .
. = ..
fn (x1 , x2 , · · · , xn ) = 0

71
72 SIST. NO LINEALES

El método de iteración de punto fijo intentara despejar de cada una de la las ecuaciones
fi (x1 , x2 , · · · , xn ) la i-esima variable tal que

x1 − g1 (x1 , x2 , · · · , xn ) = f1 (x1 , x2 , · · · , xn )
x2 − g2 (x1 , x2 , · · · , xn ) = f2 (x1 , x2 , · · · , xn )
.. .
. = ..
xn − gn (x1 , x2 , · · · , xn ) = fn (x1 , x2 , · · · , xn )

Para resolver de manera iterativa

(t+1) (t) (t)

x1 = g1 (x1 , x2 , · · · , x(t)
n )
(t+1) (t) (t)
x2 = g2 (x1 , x2 , · · · , x(t)
n )
.. ..
. = .
(t) (t)
x(t+1)
n = gn (x1 , x2 , · · · , x(t)
n )

La implementación en Matlab es
function r = Punto_Fijo(g, x1)

while 1
x2 = g(x1)
error = abs(norm((x2-x1)/x2));
if(error < 0.0001) break
else x1 = x2;
end;
end;
r = x2;

Ejemplo 1

Utilice el método de iteración de punto fijo para determinar las raı́ces del sistema de
(0) (0)
ecuaciones dado. Considere como valores iniciales x1 = 1.5 y x2 = 3.5.

x21 + x1 x2 − 10 = 0
x2 + 3x1 x22 − 57 = 0
3.2. MÉTODOS DE PUNTO FIJO 73

De acuerdo con lo descrito despejamos x1 y x2

√
f1 (x1 , x2 ) = x1 − 10 − x1 x2
r
57 − x2
f2 (x1 , x2 ) = x2 −
3x1

Despejando tenemos el siguiente sistema iterativo de ecuaciones

q
(t+1) (t) (t)
x1 = 10 − x1 x2
v
u 57 − x(t)
u
(t+1) 2
x2 = t (t+1)
3x1

La implementación en Matlab es
function x = g1(x)

x(1) = sqrt(10 - x(1)*x(2));

x(2) = sqrt((57 - x(2))/(3 * x(1)));
El proceso iterativo será:
Iteración 1
Con x(0) = [1.5, 3.5]T tenemos

q
(1) (0) (0)
p
x1 = 10 − x1 x2 = 10 − (1.5)(3.5) = 2.1794
v s
u 57 − x(0)
u
(1) 2 57 − (3.5)
x2 = t (1)
= = 2.8605
3x1 3(2.1794)

Iteración 2
Con x(1) = [2.1794, 2.8605]T tenemos

q
(2) (1) (1)
p
x1 = 10 − x1 x2 = 10 − (2.1794)(2.8605) = 1.9405
v s
u 57 − x(0)
u
(2) 2 57 − (2.8605)
x2 = t (1)
= = 3.0496
3x1 3(1.9405)
74 SIST. NO LINEALES

El resumen el proceso iterativo se muestra en la siguiente tabla, donde podemos ver que la
solución es x∗ = [2, 3]T en 9 iteraciones.
(k) (k)
k x1 x2
0 1.5000 3.5000
1 2.1794 2.8605
2 1.9405 3.0496
3 2.0205 2.9834
4 1.9930 3.0057
5 2.0024 2.9981
6 1.9992 3.0007
7 2.0003 2.9998
8 1.9999 3.0001
9 2.0000 3.0000

Para realizar la ejecución dar en Matlab

Punto_Fijo(@g1, [1.5; 3.5])

Ejemplo 2

Utilice el método de iteración de punto fijo para determinar las raı́ces del sistema de
(0) (0)
ecuaciones dado. Considere como valores iniciales x1 = 0.0 y x2 = 1.0.

x21 − 2x1 − x2 + 0.5 = 0

x21 + 4x22 − 4 = 0

De acuerdo con lo descrito despejamos x1 y x2

√
f1 (x1 , x2 ) = x1 − 2x1 + x2 − 0.5
r
4 − x21
f2 (x1 , x2 ) = x2 −
4

Despejando tenemos el siguiente sistema iterativo de ecuaciones

3.2. MÉTODOS DE PUNTO FIJO 75

q
(t+1) (t) (t)
x1 = 2x1 + x2 − 0.5
v
u 4 − x(t+1) 2
u  
(t+1)
t 1
x2 =
4

La implementación en Matlab es

function x = g2(x)

x(1) = sqrt(2*x(1) + x(2) - 0.5);

x(2) = sqrt(4 -x(1)^2)/2;

El proceso iterativo será:

Iteración 1

Con x(0) = [0, 1]T tenemos

q
(1) (0) (0)
p
x1 = 2x1 + x2 − 0.5 = 2(0) + (1) − 0.5 = 0.7071
v
u 4 − x(1) 2 r
u  
(1)
t 1 4 − (0.7071)2
x2 = = = 0.9354
4 4

Iteración 2

Con x(1) = [0.7071, 0.9354]T tenemos

q
(2) (1) (1)
p
x1 = 2x1 + x2 − 0.5 = 2(0.7071) + (0.9354) − 0.5 = 1.3600
v
u 4 − x(2) 2 r
u  
(2)
t 1 4 − (1.3600)2
x2 = = = 0.7332
4 4

El resumen el proceso iterativo se muestra en la siguiente tabla, donde podemos ver que la
solución es x∗ = [1.9007, 0.3112]T en 8 iteraciones.
76 SIST. NO LINEALES

(k) (k)
k x1 x2
0 0.0000 1.0000
1 0.7071 0.9354
2 1.3600 0.7332
3 1.7185 0.5116
4 1.8570 0.3713
5 1.8935 0.3220
6 1.8997 0.3127
7 1.9006 0.3114
8 1.9007 0.3112
Para realizar la ejecución dar en Matlab
Punto_Fijo(@g2, [0; 1])

3.3. El método de Bisecciones

Para este método debemos considerar una función unidimensional f : R → R continua

dentro de un intervalo [a, b] tal que f (a) tenga diferente signo f (a) ∗ f (b) < 0.
El proceso de decisión para encontrar la raı́z consiste en dividir el intervalo [inicio, f in] a
la mitad mitad = (inicio + f in)/2 y luego analizar las tres posibilidades que se pueden
dar.
1. Si f (inicio) y f (mitad) tienen signos opuestos, entonces hay un cero entre [inicio, mitad].
2. Si f (mitad) y f (f in) tienen signos opuestos, entonces, hay un cero en [mitad, f in].
3. Si f (mitad) es igual a cero, entonces mitad es un cero
La implementación en Matlab es:
function r = Biseccion(f, inicio, fin)
mitad = 0;
while abs((fin - inicio)/fin) > 0.0001
mitad = (fin+inicio)/2.0;

if(f(mitad) == 0)
r = mitad;
return;
end;

if f(inicio)*f(mitad) < 0
3.4. MÉTODO DE REGULA FALSI 77

fin = mitad;
else
inicio = mitad;
end;
end;

r= mitad;

3.3.1. Ejemplo

Calcular los ceros de la función f (x) = x − cos(x) utilizando el algoritmo de Bisección en

el intervalo [0, 1].

iter a c b f(a) f( c) f(b)

0 0.0000 0.5000 1.0000 -1.0000 -0.3776 0.4597
1 0.5000 0.7500 1.0000 -0.3776 0.0183 -0.3776
2 0.5000 0.6250 0.7500 -0.3776 -0.1860 -0.3776
3 0.6250 0.6875 0.7500 -0.1860 -0.0853 -0.1860
4 0.6875 0.7188 0.7500 -0.0853 -0.0339 -0.0853
5 0.7188 0.7344 0.7500 -0.0339 -0.0079 -0.0339
6 0.7344 0.7422 0.7500 -0.0079 0.0052 -0.0079
7 0.7344 0.7383 0.7422 -0.0079 -0.0013 -0.0079
8 0.7383 0.7402 0.7422 -0.0013 0.0019 -0.0013
9 0.7383 0.7393 0.7402 -0.0013 0.0003 -0.0013
10 0.7383 0.7388 0.7393 -0.0013 -0.0005 -0.0013
11 0.7388 0.7390 0.7393 -0.0005 -0.0001 -0.0005
12 0.7390 0.7391 0.7393 -0.0001 0.0001 -0.0001
13 0.7390 0.7391 0.7391 -0.0001 0.0000 -0.0001
Para ver los resultados correr

Biseccion(@f2, 0, 1)

3.4. Método de Regula Falsi

Una de las razones de la introducción de este método es que la velocidad de convergencia

del método de Bisecciones es bastante baja. En el método de Bisección se usa el punto
medio del intervalo [a, b] para llevar a cabo el siguiente paso. Suele conseguirse una mejor
78 SIST. NO LINEALES

aproximación usando el punto [c, 0] en el que la recta secante L, que pasa por los puntos
[a, f (a)] y [b, f (b)].
En la Figura 3.3 se puede ver como funciona el método. En esta figura en azul esta la
función de la cual queremos calcular el cruce por cero y en negro dos lı́neas rectas que
aproximan la solución.

Figura 3.3: Método Regula Falsi

Para calcula la ecuación de la lı́nea secante hacemos

p1 = [a, f (a)]
p2 = [b, f (b)]

y sustituimos en la ecuación de la lı́nea recta.

y–f (a) = (f (b)–f (a)) ∗ (x–a)/(b–a)

El cruce por cero de esta ecuación esta en

c = a–f (a) ∗ (b − a)/(f (b)–f (a))

Entonces el método de Bisecciones puede ser modificado, en lugar de calcular c = (a + b)/2

hacemos c = a–f (a)∗(b−a)/(f (b)–f (a)) y aplicamos las mismas tres reglas de la Bisección.
La implementación en Matlab es:
3.4. MÉTODO DE REGULA FALSI 79

function r = Regula_Falsi(f, a, b)

while abs((b - a)/a) > 0.0001

c = a - f(a)*(b-a)/(f(b) - f(a));

if(f(c) == 0)
r = c;
return;
end;

if f(a)*f(c) < 0
b = c;
else
a = c;
end;
end;

r= c;

Para ejecutar hacer

Regula_Falsi(@f2, 0, 1)

3.4.1. Ejemplo

Calcular los ceros de la función f (x) = x − cos(x) utilizando el algoritmo de regula falsi en
el intervalo [0, 1].

iter. a c b f(a) f(c ) f(b)

0 0.0000 0.6851 1.0000 -1.0000 -0.0893 0.4597
1 0.6851 0.7363 1.0000 -0.0893 -0.0047 -0.0893
2 0.7363 0.7389 1.0000 -0.0047 -0.0002 -0.0047
3 0.7389 0.7391 1.0000 -0.0002 0.0000 -0.0002
4 0.7391 0.7391 1.0000 0.0000 0.0000 0.0000
5 0.7391 0.7391 1.0000 0.0000 0.0000 0.0000
6 0.7391 0.7391 1.0000 0.0000 0.0000 0.0000
80 SIST. NO LINEALES

3.4.2. Solución de un circuito con un diodo

Consideremos un circuito de corriente alterna formado por una fuente, una resistencia y
un diodo tal como se muestra en la figura 3.4.

Figura 3.4: Circuito de corriente alterna con un diodo

La ecuación que modela al diodo esta dada por (3.6) y la curva se muestra en la figura
3.5:

id = Is (evd /Vt − 1) (3.6)

donde Is = 1−12 y Vt = 25.85−3

Para este ejemplo la ecuación que hay que resolver es:

V (t) − Ri(t) − Vd (t) = 0


 log(Id (t)/Is + 1) ∗ Vt Si Id (t) ≥ 0
Vd (t) =
V sino Id(t) < 0


La solución implementada se muestra en el siguiente código y la solución se muestra en la

Fig. 3.6
function Simula_Circuito_Diodo()
% Paramentros del Cicuito

Vmax = 10; % Volts;

R = 1.8; % ohms
3.4. MÉTODO DE REGULA FALSI 81

155
x 10
12

10

6
Id

0
5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10
Vd

Figura 3.5: Comportamiento del diodo

%tiempo de evaluación
t = [0:0.001:0.1];
%Voltaje en función del tiempo
v = Vmax*sin(377*t);

for k = 1:length(t)
I(k) = ReglaFalsa(@Circuito_Diodo, -Vmax/R, Vmax/R, [v(k), R]);
end

plot(t, v, t, I);
title(’Solución de un circuito con Diodo’);
xlabel(’tiempo s’);
legend(’Voltaje’, ’Corriente’);
end

function di = Circuito_Diodo(I, p)
V = p(1);
R = p(2);
di = V - R*I - Vd(V, I);
end
82 SIST. NO LINEALES

function r = ReglaFalsa(f, a, b, p)
%ReglaFalsa(f, a, b, p)
% f : función a evaluar
% a : inicio del intervalo
% b : fin frl intervalo
% p : parametros

while 1
c = a - f(a,p)*(b-a)/(f(b,p) - f(a,p));
if(abs(f(c,p)) <= 1e-06)
break;
end;
if f(a,p)*f(c,p) < 0
b = c;
else
a = c;
end;
end;
r = c;
end

function V = Vd(Vf, I)
% V = Vd(Vf, I)
% Vf voltaje de la Fuente
% I Corriente en el circuito

Is = 1e-12;
Vt = 25.85e-3;

if I>=0
V = log(I/Is+1)*Vt;
else
V = Vf;
end
end

3.5. El método Newton Raphson

Consideremos el sistema de ecuaciones no lineales

3.5. EL MÉTODO NEWTON RAPHSON 83

Solución de un circuito con Diodo

10
Voltaje
Corriente
8

−2

−4

−6

−8

−10
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
tiempo s

Figura 3.6: Solución al circuito con diodo

f1 (x1 , x2 , · · · , xn ) = 0
f2 (x1 , x2 , · · · , xn ) = 0
.. .
. = ..
fn (x1 , x2 , · · · , xn ) = 0

Utilizando la serie de Taylor podemos hacer una aproximación lineal para una función
(t) (t) (t) (t) (t) (t) (t) (t+1) (t)
fi (x1 +δx1 , x2 +δx2 , · · · , xn +δxn ) en un incremento δxi = xi −xi como:
84 SIST. NO LINEALES

(t) (t) (t) (t)

f1 (x1 + δx1 , x2 + δx2 , · · · , x(t) (t)
n + δxn )
(t) (t) ∂f1 (t) ∂f1 (t) ∂f1 (t)
= f1 (x1 , x2 + · · · , +x(t)
n )+ δx + δx + · · · + δx
∂x1 1 ∂x2 2 ∂xn n
(t) (t) (t) (t)
f2 (x1 + δx1 , x2 + δx2 , · · · , x(t) (t)
n + δxn )
(t) (t) ∂f2 (t) ∂f2 (t) ∂f2 (t)
= f2 (x1 , x2 + · · · , +x(t)
n )+ δx1 + δx2 + · · · + δx
∂x1 ∂x2 ∂xn n
..
.
(t) (t) (t) (t)
fn (x1 + δx1 , x2 + δx2 , · · · , x(t) (t)
n + δxn )
(t) (t) ∂fn (t) ∂fn (t) ∂fn (t)
= fn (x1 , x2 + · · · , +x(t)
n )+ δx1 + δx2 + · · · + δx
∂x1 ∂x2 ∂xn n

Si escribimos el sistema en forma matricial tenemos

(t) (t) (t) (t) (t) (t)

 
f1 (x1 + δx1 , x2 + δx2 , · · · , xn + δxn )
(t) (t) (t) (t) (t) (t)
f2 (x1 + δx1 , x2 + δx2 , · · · , xn + δxn )
 
 
 .. =
.
 
 
(t) (t) (t) (t) (t) (t)
fn (x1 + δx1 , x2 + δx2 , · · · , xn + δxn )
  ∂f (x) ∂f (x)
(t) (t) (t)
· · · ∂f∂x 1 (x) (t)
  
1 1
f1 (x1 , x2 , · · · , xn ) ∂x1 ∂x2 n
δx1
∂f2 (x) ∂f2 (x)
(t) (t) (t)
f2 (x1 , x2 , · · · , xn )   ∂x · · · ∂f∂x 2 (x)
  δx(t)
     
∂x2   .2
 1 n

 ..  + 
.. .. ..  .

. . . ··· .  .
   
   
(t) (t) (t) ∂fn (x) ∂fn (x) ∂fn (x) (t)
fn (x1 , x2 , · · · , xn ) ∂x1 ∂x2 ··· ∂xn δxn

En forma compacta

f (x(t) + δx(t) ) = f (x(t) ) + J(x(t) )δx(t)

f (x(t) + δx(t) ) = f (x(t) ) + J(x(t) )(x(t+1) − x(t) )

donde J(x) es la matriz de primeras derivada o Jacobiano.

Como queremos encontrar el cero de la función, la aproximación lineal que debemos resolver
es:
3.5. EL MÉTODO NEWTON RAPHSON 85

0 = f (x(t) ) + J(x(t) )(x(t+1) − x(t) )

donde las actualizaciones las hacemos como

h i−1
x(t+1) = x(t) − J(x(t) ) f (x(t) )

La implementación del algoritmo en Matlab es:

function r = Newton_Raphson(f, J, x1)

while 1
x2 = x1 - inv(J(x1))* f(x1);
if norm((x1-x2)./x2) < 0.0001 break;
else x1 = x2;
end;
end;
r = x2;

3.5.1. Ejemplo 1

Resolver el siguiente sistema de ecuaciones dado por y cuya gráfica se muestra en la Figura
3.7

f1 (x1 , x2 ) = x21 + x1 x2 − 10
f2 (x1 , x2 ) = x2 + 3x1 x22 − 57

El Jacobiano es  
2x1 + x2 x1
J(x) = 2
3x2 1 + 6x1 x2
y el arreglo de funciones

x21 + x1 x2 − 10
 
f (x) =
x2 + 3x1 x22 − 57

Considerando como valores iniciales x(0) = [1.5, 3.5]T tenemos:

Primer iteración
86 SIST. NO LINEALES

10

x2
0

-5

-10
-10 -5 0 5 10
x1

Figura 3.7: Curvas con la solución del sistema del ejemplo 1

" # −1 
(1)      
x1 1.5000 6.5000 1.5000 −2.5000 2.0360
(1) = − =
x2 3.5000 36.7500 32.5000 1.6250 2.8439

Segunda iteración

" # −1 
(2)      
x1 2.0360 6.9159 2.0360 −0.0644 1.9987
(2) = − =
x2 2.8439 24.2629 35.7413 −4.7562 3.0023

Tercer iteración

" # −1 
(3)      
x1 1.9987 6.9997 1.9987 −0.0045 2.0000
(3) = − =
x2 3.0023 27.0412 37.0041 0.0496 3.0000

Cuarta iteración

" # −1 
(4)      
x1 2.0000 7.0000 2.0000 −0.00000129 2.0000
(4) = − =
x2 3.0000 27.0000 37.0000 −0.00002214 3.0000
3.5. EL MÉTODO NEWTON RAPHSON 87

Note que solamente 4 iteraciones son necesarias para llegar a la solución x = [2, 3]T

La implementación en Matlab para este ejemplo son:

Para la función tenemos

function f = f1(x)

n = length(x);
f = zeros(n,1);

f(1) = x(1)^2 + x(1)*x(2) - 10;

f(2) = x(2) + 3*x(1)*x(2)^2 - 57;
end

El Jacobiano

function J = J1(x)
n = length(x);
J = zeros(n, n);

J(1,1) = 2*x(1) + x(2);

J(1,2) = x(1);
J(2,1) = 3*x(2)^2;
J(2,2) = 1 + 6*x(1)*x(2);
end

y la ejecución

>> Newton_Raphson(@f1, @J1, [1.5,3.5]’)

ans =

2.0000
3.0000

3.5.2. Ejemplo 2

Resolver el siguiente sistema de ecuaciones dado por las ecuaciones y cuya solución gráfica
se muestra en la figura 3.8
88 SIST. NO LINEALES

x2
0

-1

-2
-2 -1 0 1 2
x1

Figura 3.8: Curvas con la solución del sistema del ejemplo 2

f1 (x1 , x2 ) = x21 –2x1 –x2 + 0.5

f2 (x1 , x2 ) = x21 + 4x22 –4

El Jacobiano es  
2x1 − 2 −1
J(x) =
2x1 8x2
y el arreglo de funciones

x21 –2x1 –x2 + 0.5

 
f (x) =
x21 + 4x22 –4

Considerando como valores iniciales x(0) = [0, 1]T tenemos:

Primer iteración

" # −1 
(1)      
x1 0 −2 −1 −0.5 −0.25
(1) = − =
x2 1 0 8 0.0 1.00
3.5. EL MÉTODO NEWTON RAPHSON 89

Segunda iteración
" # −1 
(2)      
x1 −0.2500 2.5000 −1.0000 0.0625 −0.2226
(2) = − =
x2 1.0000 −0.5000 8.0000 0.0625 0.9939

Tercer iteración

" # −1 
(3)      
x1 −0.2226 −2.4451 −1.0000 −0.0008 −0.2222
(3) = − =
x2 0.9939 −0.4451 7.9512 0.0009 0.9938

Cuarta iteración

" # −1 
(4)      
x1 −0.2222 −2.4444 −1.0000 −0.00002300 −0.2222
(4) = − =
x2 0.9938 −0.4444 7.9505 0.00000038 0.9938

La implementación en Matlab para este ejemplo son:

Para la función tenemos
function f = f2(x)

n = length(x);
f = zeros(n,1);

f(1) = x(1)*x(1) - 2*x(1) - x(2) + 0.5;

f(2) = x(1)*x(1) + 4*x(2)*x(2) - 4;
end
El Jacobiano
function y = J2(x)
n = length(x);
y = zeros(n, n);
y(1,1) = 2*x(1) -2;
y(1,2) = -2;
y(2,1) = 2*x(1);
y(2,2) = 8*x(2);
end
y la ejecución
90 SIST. NO LINEALES

>> Newton_Raphson(@f2, @J2, x)

ans =

-0.2222
0.9938

3.5.3. Ejemplo

Para el circuito que se muestra en la figura 3.9, esta constituido por dos mallas y dos
cargas. Se sea calcular la corriente que circula por cada uno de los elementos y la potencia
que debe suministrar la fuente de voltaje.

Figura 3.9: Circuito eléctrico de dos mallas

La ecuaciones de Voltaje para cada una de las mallas es:

V1 − R1 (i1 + i2 ) − L1 /i1 = 0
V1 − R1 (i1 + i2 ) − R2 i2 − L2 /i2 = 0

Reorganizando las funciones tenemos

 
V i1 − R1 (i1 + i2 )i1 − L1
F = =0
V i2 − R1 (i1 + i2 )i2 − R2 i22 − L2
function f = Fun(i)
N = length(i);
f = zeros(N,1);
3.5. EL MÉTODO NEWTON RAPHSON 91

V = 100;
R1 = 2;
R2 = 3;
L1 = 100;
L2 = 200;

f(1) = V*i(1) - R1*(i(1)+i(2))*i(1) - L1;

f(2) = V*i(2) - R1*(i(1)+i(2))*i(2) - R2*i(2)^2 - L2;
end

El Jacobiano para este sistema es:

 
V − R1 (2i1 + i2 ) −i1 R1
J=
−i2 R1 V1 − R1 (i1 + 2i2 ) − 2R2 i2
function J = JacSistema(i)
V = 100;
R1 = 2;
R2 = 3;
%L1 = 100;
%L2 = 200;

J(1,1) = V - R1*(2*i(1) + i(2)) ;

J(1,2) = -i(1)*R1;
J(2,1) = -i(2)*R1;
J(2,2) = V - (i(1) + 2*i(2))*R1 - 2*i(2)*R2 ;
end

La solución considerando un valor inicial de corrientes I (0) = [0, 0]T :

>> I = NewtonRaphson(@FunSistema, @JacSistema, [0,0]’)

1.- 1 2

2.- 1.0720 2.3114

3.- 1.0728 2.3185

4.- 1.0728 2.3185

I =
92 SIST. NO LINEALES

1.0728
2.3185

Potencia de la fuente PV = V (i1 + i2 ) = 100(1.0728 + 2.3185) = +339.1282

Potencia consumida por R1 es PR1 = −R1 (i1 + i2 )2 = −2(1.0728 + 2.3185)2 = −23.0016
Potencia consumida por R2 es PR1 = −R2 i22 = −3(2.3185)2 = −16.1266
La potencia consumida por las cargas es L1 + L2 = −300.
La suma de la potencia es cero.

3.6. Convergencia del método de Newton Raphson

Aunque el método de Newton-Raphson en general es muy eficiente, hay situaciones en que

se comporta en forma deficiente. Un caso especial, raı́ces múltiples.

Ejemplo

Determinar la raı́z de la función f (x) = x10 –1.

La solución utilizando el método de Newton-Raphson queda:

f (x(k) )
x(k+1) = x(k) −
f 0 (x(k+1) )

Sustituyendo valores tenemos

10
(k+1) (k) x(k) −1
x =x −
9
10 x(k)

Y la solución numérica es:

3.6. CONVERGENCIA DEL MÉTODO DE NEWTON RAPHSON 93

x(k) f (x(k) ) f 0 (x(k) )

0.5000 -0.9990 0.0195
51.6500 135114904483914000.0000 26159710451871000.0000
46.4850 47111654129711500.0000 10134807815362300.0000
41.8365 16426818072478500.0000 3926432199748670.0000
37.6529 5727677301318310.0000 1521180282851980.0000
33.8876 1997117586819850.0000 589336409039672.0000
30.4988 696351844868619.0000 228320999775654.0000
27.4489 242802875029547.0000 88456233382052.8000
24.7040 84660127717097.5000 34269757191973.2000
22.2336 29519161271064.1000 13276806089225.7000
20.0103 10292695105054.7000 5143706707446.1600
18.0092 3588840873655.1100 1992777367871.5700
16.2083 1251351437592.9200 772042782329.1500
14.5875 436319267276.5290 299105192259.1190
13.1287 152135121499.2910 115879479847.7330
11.8159 53046236848.5329 44894084747.9692
10.6343 18496079117.2577 17392888266.5936
9.5708 6449184014.3077 6738361277.7304
8.6138 2248691421.7628 2610579221.6818
7.7524 784070216.9426 1011391879.0870
94 SIST. NO LINEALES
Ajuste de curvas

4.1. Regresión lineal por el método de mı́nimos cuadrados

El ejemplo más simple de una aproximación por mı́nimos cuadrados es mediante el ajuste
de un conjunto de pares de observaciones: [x1 , y1 ], [x2 , y2 ], [x3 , y3 ],..., [xn , yn ] a una lı́nea
recta. La expresión matemática de la lı́nea recta es:

y = a0 + a1 x + e

donde a0 y a1 , son los coeficientes que representan el cruce con el eje y y la pendiente de
la lı́nea y e representa el error de nuestra aproximación. En la figura ??, se muestra un
ejemplo de puntos a los que se les desea ajustar a una linea recta.

y
6

x
2 3 4 5 6

Figura 4.10: Un conjunto de puntos

95
96 AJUSTE DE CURVAS

Una estrategia, para ajustar a la mejor lı́nea, es minimizar la suma al cuadrado de los
errores para todos los datos disponibles

n
X n
X
E(a) = e2i = (yi − a0 − a1 xi )2
i=1 i=1

Esta misma ecuación la podemos escribir en forma matricial como

 
e1
 e2 
E(a) = [e1 , e2 , · · · , en ] 
 
.. 
 . 
en

Los elementos del vector de error ei en forma matricial queda como

 
a0
ei = yi − [1, xi ]
a1

y en general

     
e1 y1 1 x1
 e2   y2   1 x2  
a0
    
 e3 =
  y3 −
  1 x3 
.. ..  a1
 
 ..   ..  
 .   .   . . 
en yn 1 xn

En forma compacta tenemos que e = y–M a, donde M es la matriz de coordenadas en x y

a el vector de parámetros.

4.1.1. Ajuste por mı́nimos cuadrados

Si queremos encontrar el vector de parámetros a que minimiza nuestra suma de cuadrados,

tenemos que calcular la derivada de la función de error respecto al vector de parámetros e
igualar a cero.

∂E(a) ∂
= [y − M a]T [y − M a] = 0
∂a ∂a
∂E(a)
= −2M T [y − M a] = 0
∂a
4.1. REGRESIÓN LINEAL POR EL MÉTODO DE MÍNIMOS CUADRADOS 97

El valor del vector de parámetros a los calculamos resolviendo el siguiente sistema de

ecuaciones

[M T M ]a = M T y

Es común encontrar la solución de este sistema como:

 Pn Pn    Pn 
1 x a0 yi
Pni=1 Pni=1 2i = Pn
i=1
i=1 xi i=1 xi a1 i=1 xi yi

La implementación en Matlab es:

function a = AjustePolinomial (x, y, grado)

N = length(x);
M = ones(N, 1);

for k=1:grado
M = [M,x.^k];
end

a=inv(M’*M)*M’*y;
end

4.1.2. Ejemplo 1

Hacer el ajuste a una lı́nea recta de los siguientes valores

x y
1.00 0.50
2.00 2.50
3.00 2.00
4.00 4.00
5.00 3.50
6.00 6.00
7.00 5.50

La matriz M queda como

98 AJUSTE DE CURVAS

 
1 1.00

 1 2.00 


 1 3.00 

M =
 1 4.00 


 1 5.00 

 1 6.00 
1 7.00

 
1 1.00

 1 2.00 

  1 3.00   
1 1 1 1 1 1 1 7.00 28.00
MT M =
 
 1 4.00 =
1.00 2.00 3.00 4.00 5.00 6.00 7.00   28.00 140.00

 1 5.00 

 1 6.00 
1 7.00

 
0.50

 2.50 

  2.00   
1 1 1 1 1 1 1 24.00
MT y =
 
 4.00 =
1.00 2.00 3.00 4.00 5.00 6.00 7.00   119.50

 3.50 

 6.00 
5.50

Aplicando las formulas anteriores, tenemos que el sistema de ecuaciones a resolver es

    
7.00 28.00 a0 24.00
=
28.00 140.00 a1 119.50

La solución es a = [0.07142, 0.8392] y en la figura 4.11 se muestra el ajuste encontrado

Para ejecutar dar

AjustePolinomial([1,2,3,4,5,6,7]’, [0.5, 2.5, 2, 4, 3.5, 6, 5.5]’,1)

4.1. REGRESIÓN LINEAL POR EL MÉTODO DE MÍNIMOS CUADRADOS 99

y
6

x
1 2 3 4 5 6 7

Figura 4.11: Un conjunto de puntos y la linea recta ajustada

4.1.3. Regresión polinomial

Podemos generalizar el caso de la regresión lineal y extenderla a cualquier polinomio de

orden m, hacemos la siguiente representación.
 
 2 3 m
 a0  
1 x1 x1 x1 · · · x1 y0
 1 x2 x2 x3 · · · xm  
 a1 
2 2 2
  y1 

 1 x3 x2 x3 · · · xm  
 a2   
3 3 3  =
  y2 

 .. .. .. .. .. ..   a3   ..


 . . . . . .  ..   . 
 . 
1 xn xn xn · · · xm
2 3
n yn
am

La cual podemos escribir de manera compacta como

Ma = y

Si pre-multiplicamos ambos lados de la ecuación tenemos

[M T M ]a = M T y (4.7)

El sistema se resuelve para a o simplemente se calcula haciendo

a = [M T M ]−1 [M T y]
100 AJUSTE DE CURVAS

4.1.4. Ejemplo 2

Ajustar a un polinomio de segundo orden los datos en la siguiente tabla.

x y
0.00 2.10
1.00 7.70
2.00 13.60
3.00 27.20
4.00 40.90
5.00 61.10
En este caso por ser un polinomio de segundo grado tenemos una matriz M dada por
 
1 0 0

 1 1 1 

 1 2 4 
M = 

 1 3 9 

 1 4 16 
1 5 25

El sistema de ecuaciones que debemos resolver es:

    
6.00 15.00 55.00 a0 152.60
 15.00 55.00 225.00   a1  =  585.60 
55.00 225.00 979.00 a2 2488.80

La solución del sistema es a = [2.4785, 2.3592, 1.8607] y el ajuste da como resultado el

polinomio p(x) = 2.4785+2.3592x+1.8607x2. En la figura 4.12 se muestra la aproximación
a un polinomio de segundo orden.
Para correr en Matlab hacer
AjustePolinomial([0,1,2,3,4,5]’, [2.1, 7.70, 13.60, 27.20, 40.90, 61.10]’,2)

4.2. Interpolación lineal

El modo más simple de interpolación es conectar dos puntos con una lı́nea recta. Esta técni-
ca, llamada interpolación lineal, la podemos representar por la siguiente formulación, la cual
se calcula considerando que se tienen dos puntos P1 = [x1 , f (x1 )]T y P0 = [x0 , f (x0 )]T
4.2. INTERPOLACIÓN LINEAL 101

60

50

40

30

20

x
1 2 3 4 5

Figura 4.12: Un conjunto de puntos y su aproximación cuadrática

f (x) − f (x0 ) f (x1 ) − f (x0 ))

=
x − x0 x1 − x0

La cual da lugar a la siguiente ecuación

f (x1 ) − f (x0 )
f (x) = f (x0 ) + (x − x0 )
x1 − x0

La implementación para una función de interpolación lineal para un conjunto de puntos

es:

function y = Interpolacion_Lineal(x0, x, y)

N =length(x);

for k=1:N-1
if x0 >= x(k) && x0 <= x(k+1) break;
end
end

y = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k))*(x0-x(k));

end
102 AJUSTE DE CURVAS

4.2.1. Ejemplo 1

Dados los puntos

x y
0.00 2.10
1.00 7.70
2.00 13.60
3.00 27.20
4.00 40.90
5.00 61.10

Implementar el código correspondiente para graficar los puntos intermedios a los valores
dados utilizando interpolación lineal

x = [0,1,2,3,4,5];
y = [2.1, 7.70, 13.60, 27.20, 40.90, 61.10];

x_nva = 0:0.01:5;

N = length(x_nva);
y_nva = zeros(N, 1);

for k=1:N
y_nva(k) = Interpolacion_Lineal(x_nva(k), x, y);
end;

plot(x_nva, y_nva, ’r-.’, x, y, ’k*’);

xlabel(’x’);
ylabel(’f(x)’);

En la figura 4.13, se muestran los puntos correspondientes calculados de acuerdo con el

código implementado

4.3. Interpolación cuadrática

El error que se observa en la figura 4.13, se debe a que suponemos que los puntos se unen por
lineas. Una manera resolver este problema es suponer una función polinomial cuadrática
tomando al menos tres puntos para tal propósito. Una forma de introducir tres puntos en
una función cuadrática es mediante la siguiente ecuación
4.3. INTERPOLACIÓN CUADRÁTICA 103

70

60

50

40
f(x)

30

20

10

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x

Figura 4.13: Un conjunto de puntos y su interpolación lineal

q(x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 )

Para calcular los valores b0 , b1 y b2 necesitamos introducir información de tres puntos

[x0 , f (x0 )], [x1 , f (x1 )] y [x2 , f (x2 )]. Si sustituimos en primer x0 tenemos que

f (x0 ) = b0 + b1 (x0 − x0 ) + b2 (x0 − x0 )(x − x1 )

por lo tanto
b0 = f (x0 )

Para calcular b1 utilizamos el segundo punto ası́ tenemos

f (x0 ) = b0 + b1 (x1 − x0 ) + b2 (x1 − x0 )(x1 − x1 )

despejando tenemos

f (x1 ) − f (x0 )
b1 =
x1 − x0
104 AJUSTE DE CURVAS

Para el calculo de b2 sustituimos los valores de b0 , b1 y [x2 , f (x2 )] de la siguiente mane-

ra

f (x1 ) − f (x0 )
f (x0 ) = f (x0 ) + (x2 − x0 ) + b2 (x2 − x0 )(x2 − x1 )
x1 − x0

despejando tenemos

f (x2 )−f (x1 ) f (x1 )−f (x0 )

x2 −x1 − x1 −x0
b2 =
x2 − x0

4.3.1. Ejemplo 1

Calcular el logaritmo de x = 2 dados tres puntos

x0 = 1 f (x0 ) = 0
x1 = 4 f (x1 ) = 1.386294
x2 = 6 f (x2 ) = 1.791759
Para resolver tenemos

b0 = 0

f (x1 ) − f (x0 ) 1.386294 − 0

b1 = = = 0.4620981
x1 − x0 4−1

El coeficiente b2 tenemos

f (x2 )−f (x1 ) f (x1 )−f (x0 )

x2 −x1 − x1 −x0
b2 =
x2 − x0

1.791759−1.386294
6−4 − 0.4620981
b2 = = −0.0518731
6−1

Finalmente para calcular el valor hacemos

q(x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 )
4.3. INTERPOLACIÓN CUADRÁTICA 105

q(x) = 0 + 0.4620981(x − 1) − 0.0518731(x − 1)(x − 4)

q(2) = 0 + 0.4620981(2 − 1) − 0.0518731(2 − 1)(2 − 4) = 0.5658444

En Matlab dar
Interpolacion_Newton(2, [1, 4, 6] , [0, 1.386294, 1.791759])

ans =

0.5658

4.3.2. Ejemplo 2

Dados los puntos

x y
0.00 2.10
1.00 7.70
2.00 13.60
3.00 27.20
4.00 40.90
5.00 61.10
Hacer una interpolación cuadrática de los datos y comparar con la interpolación lineal
x = [0,1,2,3,4,5];
y = [2.1, 7.70, 13.60, 27.20, 40.90, 61.10];

x_nva = 0:0.01:5;

grado = 2;

N = length(x_nva);
M = length(x);
y_nva = zeros(N, 1);

for n=1:N
for m=1:M-1
if x_nva(n) >= x(m) && x_nva(n) <= x(m+1) break;
106 AJUSTE DE CURVAS

end;
end;
l = m+grado;
if l > M
l = M;
end;

y_nva(n) = Interpolacion_Newton(x_nva(n), x(m:l), y(m:l));

end;

plot(x_nva, y_nva, ’r-.’, x, y, ’k*’);

xlabel(’x’);
ylabel(’f(x)’);

En la figura 4.14, se muestran los puntos correspondientes calculados de acuerdo con el

código implementado

70

60

50

40
f(x)

30

20

10

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x

Figura 4.14: Un conjunto de puntos y su interpolación cuadrática

4.4. FORMULAS DE INTERPOLACIÓN DE NEWTON 107

4.4. Formulas de interpolación de Newton

Las formulaciones anteriores pueden ser generalizadas para ajustar un polinomio de n-ésimo
orden a n + 1 datos. El polinomio de n-ésimo orden es

fn (x) = b0 + b1 (x − x0 ) + b2 (x − x0)(x − x1 ) + · · · + bn (x − x0 )(x − x1 ) · · · (x − xn−1 ) (4.8)

Como se hizo antes con las interpolaciones lineales y cuadráticas, los puntos de los datos
evaluaban los coeficientes b0 , b1 , · · · , bn . Para un polinomio de n-ésimo orden se requiere
n + 1 puntos: [x0 , f (x0 )], [x1 , f (x1 )], ..., [xn , f (xn )]. Usamos estos datos y las siguientes
ecuaciones para evaluar los coeficientes

b0 = f (x0 )
b1 = f [x1, x0 ]
b2 = f [x2 , x1 , x0 ]
..
.
bn = f [xn , xn−1 , · · · , x1 , x0 ]

donde las evaluaciones de la función puestas entre paréntesis son diferencias divididas
finitas. Por ejemplo, la primera diferencia dividida finita se representa por lo general co-
mo

f (xi ) − f (xj )
f [xi , xj ] =
xi − xj

La segunda diferencia dividida finita, la cual representa la diferencia de las dos primeras
diferencias divididas, se expresa por lo general como

f [xi , xj ] − f [xj , xk ]
f [xi , xj , xk ] =
xi − xk

En forma similar, la n-ésima diferencia dividida finja es

f [xn , xn−1 , · · · , x1 ] − f [xn−1 , xn−2 , · · · , x1 , x0 ]

f [xn , xn−1 , · · · , x1 , x0 ] =
xn − x0
108 AJUSTE DE CURVAS

Estas diferencias pueden usarse para evaluar los coeficientes en la ecuación 4.8, las cuales
entonces se sustituirán en la siguiente ecuación, para ası́ obtener el polinomio de interpo-
lación

fn (x) = f (x0 ) + f [x1 , x0 ](x − x0 ) + f [x2 , x1 , x0 ](x − x0)(x − x1 ) +

· · · +f [xn , xn−1 , · · · , x1 , x0 ](x − x0 )(x − x1 ) · · · (x − xn−1 )

En la siguiente tabla se muestra el procedimiento para el calculo de las diferencias finitas

divididas.
i xi f (xi ) Primero Segundo Tercero
0 x0 f (x0 ) f [x1 , x0 ] f [x2 , x1 , x0 ] f [x3 , x2 , x1 , x0 ]
1 x1 f (x1 ) f [x2 , x1 ] f [x3 , x2 , x1 ]
2 x2 f (x2 ) f [x3 , x2 ]
3 x3 f (x3 )
La implementación de este algoritmo en Matlab es:
function yint = Interpolacion_Newton(xi, x, y)
N = length(x);

f = zeros(N, N);

for n=1:N
f(n, 1) = y(n);
end;

for m = 2:N
for n = 1:N+1-m
f(n,m)=(f(n+1,m-1)-f(n, m-1))/(x(n+m-1) - x(n));
end;
end;

yint = f(1,1);
dx = 1;

for n=2:N
dx = dx*(xi-x(n-1));
yint = yint + f(1,n)*dx;
end;
4.4. FORMULAS DE INTERPOLACIÓN DE NEWTON 109

end

4.4.1. Ejemplo 1

Calcular el logaritmo de 2 utilizando la interpolación de Newton y los valores [1, 0], [4, 1.386294],
[6, 1.791759] y [5, 1.609438]
La solución la llevaremos a cabo de forma tabular
i xi f (xi ) Primero Segundo Tercero
0 1 0.0000 0.4621 -0.0519 0.0079
1 4 1.3863 0.2027 -0.0204
2 6 1.7918 0.1823
3 5 1.6094
Para 4 puntos el polinomio de Newton de tercer orden es

f3 (x) = 0 + 0.462098(x − 1) − 0.05187311(x − 1)(x − 4) + 0.00786529(x − 1)(x − 4)(x − 6)

y el valor del logaritmo es f3 (2) = 0.6287686

4.4.2. Ejemplo 2

Dados los puntos

x y
0.00 2.10
1.00 7.70
2.00 13.60
3.00 27.20
4.00 40.90
5.00 61.10
Hacer una interpolación utilizando polinomios de Newton de quinto grado de los datos y
comparar con la interpolación lineal y cuadrática
x = [0,1,2,3,4,5];
y = [2.1, 7.70, 13.60, 27.20, 40.90, 61.10];

x_nva = 0:0.01:5;
110 AJUSTE DE CURVAS

grado = 5;

N = length(x_nva);
M = length(x);
y_nva = zeros(N, 1);

for n=1:N
for m=1:M-1
if x_nva(n) >= x(m) && x_nva(n) <= x(m+1) break;
end;
end;
l = m+grado;
if l > M
l = M;
end;

y_nva(n) = Interpolacion_Newton(x_nva(n), x(m:l), y(m:l));

end;

plot(x_nva, y_nva, ’r-.’, x, y, ’k*’);

xlabel(’x’);
ylabel(’f(x)’);

En la figura 4.15, se muestran los puntos correspondientes calculados de acuerdo con el

código implementado

4.5. Interpolación de Polinomios de Lagrange

La interpolación de polinomios de Lagrange es simplemente una reformulación del polino-

mio de Newton que evita el cálculo por diferencias divididas. Se puede expresar de manera
concisa como

n
X
fn (x) = Li (x)f (xi )
i=0

donde
4.5. INTERPOLACIÓN DE POLINOMIOS DE LAGRANGE 111

70

60

50

40
f(x)

30

20

10

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x

Figura 4.15: Un conjunto de puntos y su interpolación con Polinomios de Newton de Quinto

grado

n
Y x − xj
Li (x) =
xi − xj
j=0
j 6= i

Por ejemplo la versión lineal n = 1 es

x − x1 x − x0
f1 (x) = f (x0 ) + f (x1 )
x0 − x1 x1 − x0

y la versión de segundo orden es

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

f2 (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

El código en Matlab para este algoritmo es

function yint = Interpolacion_Lagrange(xint, x, y)

suma = 0;
112 AJUSTE DE CURVAS

N = length(x);
for n=1:N
producto = y(n);
for m = 1:N
if n~=m
producto = producto*(xint -x(m))/(x(n) - x(m));
end;
end;
suma = suma + producto;
end;
yint = suma;
end

4.5.1. Ejemplo 1

Use una interpolación del polinomio de Lagrange de primer y segundo orden para evaluar
ln(2) con base en los siguientes datos [1, 0], [4, 1.386294] y [6, 1.791759]
La solución lineal n = 1 es

2−4 2−1
f1 (x) = 0+ 1.386294 = 0.4620981
1−4 4−1

y la solución de segundo orden es

(2 − 4)(2 − 6) (2 − 1)(2 − 6) (2 − 1)(2 − 4)

f2 (x) = 0+ 1.386294 + 1.791760 = 0.5658444
(1 − 4)(1 − 6) (4 − 1)(4 − 6) (6 − 1)(6 − 4)

Como era de esperar, ambos resultados concuerdan con los que se obtuvieron utilizando
interpolación utilizando polinomios de Newton.
Diferenciación e Integración

5.1. Diferenciación por diferencias divididas finitas atrás,

adelante y centrales de exactitud simple

Considerando la definición de la derivada

df (x) f (x) − f (x − h)
= lı́m
dx h→0 h

entonces para valores pequeños de h podemos hacer una aproximación de la derivada

haciendo

df (x) f (x) − f (x − h)
≈
dx h

En el caso de considerar la diferencia hacia adelante podemos hacer

df (x) f (x + h) − f (x)
≈
dx h
Considerando que los incrementos se dan hacia adelante y hacia atrás, podemos dar una
representación de la derivada centrada como

df (x) f (x + h) − f (x − h)
≈
dx 2h

5.2. Integración por el método de barras

En general todos los métodos de integración numérica hacen uso de la interpretación de la

integral propia, esta interpretación nos dice que la integral en un intervalo a y b es el área

113
114 DIFERENCIACIÓN E INTEGRACIÓN

bajo la curva. Ası́ podemos hacer este cálculo del área suponiendo que nuestros rectángulos
son de altura f (a). Recordemos que el área de un rectángulo esta dado como

A = (b − a)f (a)

por lo que la integral la podemos aproximar por

Z b Z b
I= f (x)dx ≈ f (a)dx = f (a)x|ba = f (a) ∗ (b − a)
a a

5.2.1. Ejemplo

Para la función f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5 , calcular la integral
utilizando el método de barras en el intervalo a = 0 b = 0.8 .

I ≈ f (0) ∗ (0.8 − 0) = 0.2 × 0.8 = 0.16

En Matlab
>> Barras(@f1, 0, 0.8, 1)

ans =

0.1600
Para una mejor aproximación podemos hacer N divisiones de la región de tamaño h y
aproximar la integral como:

Z b N
X −1
f (x)dx ≈ f (a + kh) ∗ h
a k=0

Esta aproximación da lugar al método de barras cuya codificación en Matlab se presenta

a continuación
function I = Barras(f, a, b, N)

h = abs(b-a)/N;

suma = 0;
for k=0:N-1
5.3. INTEGRACIÓN POR EL MÉTODO DE TRAPEZOIDES 115

suma = suma + f(a+h*k);

end;

I = suma*h;

end

5.3. Integración por el método de trapezoides

Suponer que todos los rectángulos son de altura constante es una suposición fuerte. En
lugar de ello supondremos, para este método, que los rectángulos varı́an linealmente de
acuerdo con la siguiente ecuación

(x − a)
f (x) = f (a) + (f (b) − f (a))
b−a
A partir de esto podemos hacer la aproximación de nuestra integral
Z b Z b 
(x − a)
I= f (x)dx ≈ f (a) + (f (b) − f (a)) dx
a a b−a

resolviendo la integral tenemos

b
(x − a)2


I ≈ f (a)x + (f (b) − f (a)) 
2(b − a) a

f (b) + f (a)
I≈ (b − a)
2
La cual es el área de un trapecio de altura menor f (a) y altura mayor f (b).

5.3.1. Ejemplo

Para la función f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5 , calcular la integral
utilizando el método de la Regla Trapezoidal en el intervalo a = 0 b = 0.8 .

f (b) + f (a) (0.2320 + 0.2)(0.8 − 0)

I≈ (b − a) = = 0.1728
2 2
En Matlab la solución es
116 DIFERENCIACIÓN E INTEGRACIÓN

>> Regla_Trapezoidal(@f1, 0, 0.8, 1)

ans =

0.1728
Si aplicamos la formula para N divisiones de tamaño h tenemos que la regla trapezoidal
la podemos calcular como

PN −1
f (a) + 2 k=1 f (a + kh) + f (b)
I ≈ (b − a)
2N

La implementación en Matlab para la regla trapezoidal

function I = Regla_Trapezoidal(f, a, b, N)

h = abs(b-a)/N;

suma = f(a);

for k=1:N-1
suma = suma + 2*f(a+h*k);
end;
suma = suma + f(b);

I = (b-a)*suma/(2*N);

end

5.4. Integración por el método de regla Simpson 1/3

En general cuando se utiliza una función de interpolación de mayor grado la diferencia entre
los valores estimados y los reales se vuelve pequeña. Para este método en lugar se suponer
que la función se aproxima por lineas haremos la aproximación utilizando un polinomio
de segundo orden. Pare este método utilizaremos los polinomios de Lagrange de segundo
orden dado por

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

f (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )
5.4. INTEGRACIÓN POR EL MÉTODO DE REGLA SIMPSON 1/3 117

Aplicando la integración de la función de interpolación tenemos:

Z b
f (x0 ) + 4f (x1 ) + f (x2 )
I≈ f (x)dx = (b − a)
a 6

5.4.1. Ejemplo

Para la función f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5 , calcular la integral
utilizando el método Simpson 1/3 en el intervalo a = 0 b = 0.8 .

f (x0 ) + 4f (x1 ) + f (x2 ) f (0) + 4f (0.4) + f (0.8)

I ≈ (b − a) = (0.8 − 0) =
6 6
(0.2 + 4 × 2.4560 + 0.2320)
I ≈ (0.8 − 0) = 1.3675
6
En Matlab la solución es
>> Simpson13(@f1, 0, 0.8, 2)

ans =

1.3675
La implementación para N divisiones de tamaño h es:

PN −1 PN −1
f (xa ) + 4 k=1,3,5,··· f (a + k ∗ h) + 2 k=2,4,6,··· f (a + k ∗ h) + f (b)
I ≈ (b − a)
3N

y la implementación en Matlab es:

function I = Simpson13(f, a, b, N)

h = abs(b-a)/N;

s_par = 0;
s_impar = 0;

for k=1:N-1
if mod(k,2) == 1
s_impar = s_impar + f(a+h*k);
else
118 DIFERENCIACIÓN E INTEGRACIÓN

s_par = s_par + f(a+h*k);

end;
end;

I = (b-a)*(f(a) + 4*s_impar + 2*s_par + f(b))/(3*N);

end

5.5. Integración por el método de regla Simpson 3/8

Para este método utilizaremos como función de integración un polinomio de Lagrange de

tercer orden dado por

(x − x1 )(x − x2 )(x − x3 )
f (x) = f (x0 )
(x0 − x1 )(x0 − x2 )(x0 − x3 )
(x − x0 )(x − x2 )(x − x3 )
+ f (x1 )
(x1 − x0 )(x1 − x2 )(x1 − x3 )
(x − x0 )(x − x1 )(x − x3 )
+ f (x2 )
(x2 − x0 )(x2 − x1 )(x2 − x3 )
(x − x0 )(x − x1 )(x − x2 )
+ f (x3 )
(x3 − x0 )(x3 − x1 )(x3 − x2 )

Resolviendo para esta función polinomial

Z b
f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )
I≈ f (x)dx = (b − a)
a 8

5.5.1. Ejemplo

Para la función f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5 , calcular la integral
utilizando la regla Simpson 3/8 en el intervalo a = 0 b = 0.8 .

f (0) + 3f (0.2667) + 3f (0.5333) + f (0.8)

I ≈ (0.8 − 0) =
8

(0.2 + 3 × 1.4327 + 3 × 3.4872 + 0.2320)

I ≈ (0.8 − 0) = 1.5192
8
5.6. EJEMPLOS 119

En Matlab la solución es
>> Simpson38(@f1, 0, 0.8, 1)

ans =

1.5192
La implementación en Matlab de la regla de Simpson 3/8 para N intervalos es
function I = Simpson38(f, a, b, N)
h = abs(b-a)/N;

suma =0;

for k=0:N-1
x0 = a + h*k;
x1 = x0 + h/3;
x2 = x0 + 2*h/3;
x3 = x0 + h;
suma = suma + f(x0) + 3*f(x1) + 3*f(x2) + f(x3);
end;

I = h*suma/8;
end

5.6. Ejemplos

5.6.1. Ejemplo 1

Dada la función

f (x) = 0.2 + 25x − 200x2 + 675x3 − 900x4 + 400x5

Hacer la solución de las integral en el intervalo a = 0 b = 0.8 con N = 1 aplicando los

métodos vistos y comparar con la solución real dada por la ecuación.
Comenzamos por calcular la integral de f (x)

Z b
25 2 200 3 675 4 900 5 400 6
I(x) = f (x)dx = 0.2x + x − x + x − x + x
a 2 3 4 5 6
120 DIFERENCIACIÓN E INTEGRACIÓN

La solución para cada uno de los métodos es

Método Solución error
Barras 01600 1.4805
R Trapezoidal 0.1728 1.4677
Simpson 1/3 1.3675 0.2731
Simpson 3/8 1.5192 0.1214
Note que el mejor desempeño lo tiene el método de Simpson 3/8 en virtud de utilizar
polinomios de grado 3 para llevar a cabo la integración.

5.6.2. Ejemplo 2

Repetir el ejemplo anterior para N = 100.

Método Solución error
Barras 1.6401 0.0004
R Trapezoidal 1.6403 0.0003
Simpson 1/3 1.6405 0.0000
Simpson 3/8 1.6405 0.0000
Note que a medida que se utilizan mas intervalos la solución en mucho mejor acercando a
la solución real.
Ecuaciones diferenciales ordinarias

En este capı́tulo se lleva a cabo la solución de ecuaciones diferenciales de la forma

dy
= f (x, y)
dx

En general los métodos que veremos hacen una formulación lineal de la forma

y (k+1) = y (k) + φh

donde φ es una estimación de la pendiente de la lı́nea recta y utilizaremos esta linea para
hacer una predicción de los valores de la integral.

6.1. Integración por el método de Euler

La primera derivada proporciona una estimación directa de la pendiente en x(k) dada

por

φ = f (x(k) , y (k) )

donde f (x(k) , y (k) ) es la ecuación diferencial evaluada en x(k) y y (k) . Utilizando esta esti-
mación tenemos

y (k+1) = y (k) + f (x(k) , y (k) )h

La ecuación anterior es conocida como el método de Euler (o de Euler-Cauchy o de punto

medio). Su implementación en Matlab es:

121
122 ECUACIONES DIFERENCIALES ORDINARIAS

function solucion = Euler(xini, h, N, f)

D = length(xini);
x = xini(1);
y = xini(2:D);

solucion = zeros(N+1,D);
solucion(1,:) = [x,y];

for k = 2:N+1
x = x + h;
y = y + f(x, y)*h;
solucion(k,:)=[x, y];
end;
end
donde x ini el el vector de valores iniciales, h es el incremento, N es el conjunto de valores
donde se integra y f son el conjunto de ecuaciones diferenciales.

6.2. Integración por el método de Heun con solo uno y con

varios predictores

Un método para mejorar la estimación de la pendiente involucra la determinación de dos

derivadas para cada uno de los intervalos (una en el punto inicial y otra en el punto final).
Las dos derivadas se promedian con el fin de obtener una estimación mejorada de las
pendientes para todo el intervalo.
Recordemos que el método de Euler, la pendiente al inicio de un intervalo es

y 0(k) = f (x(k) , y (k) )

y se usa para extraporlar linealmente a y (k+1)

yp(k+1) = y (k) + f (x(k) , y (k) )h

La pendiente en el punto final del intervalo la podemos calcular como:

y 0(k+1) = f (x(k+1) , yp(k+1) )

y la pendiente promedio la podemos calcular como

6.3. INTEGRACIÓN POR EL MÉTODO DEL PUNTO MEDIO 123

y 0(k+1) + y 0(k)
ȳ 0 =
2

De lo anterior podemos decir que el método tiene dos pasos, uno de predicción que es una
iteración de Euler y un Corrector donde se utiliza la estimación de la pendiente al final del
bloque. Estos pasos son:
Predicción
yp(k+1) = y (k) + f (x(k) , y (k) )h

Corrección
(k+1)
f (x(k) , y (k) ) + f (x(k+1) , yp )
y (k+1) = y (k) + h
2

La implementación en Matlab para el Método de Heun es:

function solucion = Heun(xini, h, N, f)
D = length(xini);
x = xini(1);
y = xini(2:D);

solucion = zeros(N+1,D);
solucion(1,:) = [x,y];

for k = 2:N+1
%Predictor
k1 = f(x, y );
k2 = f(x+h, y + k1*h);
%corrector
x = x + h;
y = y + (k1 + k2 )*h/2.0;
solucion(k,:) = [x,y];
end;
end

6.3. Integración por el método del punto medio

Esta técnica usa el método de Euler para predecir un valor de y en el punto medio del
intervalo
124 ECUACIONES DIFERENCIALES ORDINARIAS

y (k+1/2) = y (k) + f (x(k) , y (k) )h/2

Después este valor es utilizado para calcular la pendiente en el punto medio y hacer la
solución en todo el intervalo como

y (k+1) = y (k) + f (x(k+1/2) , y (k+1/2) )h

La implementación en Matlab es

function solucion = Punto_Medio(xini, h, N, f)

D = length(xini);
x = xini(1);
y = xini(2:D);

solucion = zeros(N+1,D);
solucion(1,:) = [x,y];

for k = 2:N+1

xm = x + h/2;
ym = y + f(x, y)*h/2;

x = x + h;
y = y + f(xm, ym)*h;
solucion(k,:) = [x,y];
end;
end

6.4. Fórmulas de integración por algunos métodos de Runge-

Kutta de segundo orden clásicos

Los algoritmos de Runge-Kutta son los métodos de integración con mejor desempeño y
exactitud. En el caso del método de Runge-Kutta tenemos calcular dos valores k1 y k2
utilizando

k1 = f (x(k) , y (k) )
6.5. RUNGE-KUTTA 3ER ORDEN 125

3 3
k2 = f (x(k) + h, k1 h)
4 4

La actualización se lleva a cabo haciendo

 
(k+1) (k) 1 2
y =y + k1 + k2 h
3 3
function solucion = RK(xini, h, N, f)
D = length(xini);
x = xini(1);
y = xini(2:D);

solucion = zeros(N+1,D);
solucion(1,:) = [x,y];

for k = 2:N+1
k1 = f(x, y);
k2 = f(x + 3*h/4, y+3*h*k1/4);

x = x + h;
y = y + (k1/3 + 2*k2/3)*h;
solucion(k,:) = [x,y];
end;
end

6.5. Fórmula de integración por el método de Runge-Kutta

clásico de tercer orden

En el caso del método de Runge-Kutta de tercer orden tenemos que calcular tres valores
k1 , k2 y k3 utilizando las siguientes formulas

k1 = f (x(k) , y (k) )

1 1
k2 = f (x(k) + h, k1 h)
2 2
126 ECUACIONES DIFERENCIALES ORDINARIAS

k3 = f (x(k) + h, y (k) − k1 h + 2k2 h)

La actualización se lleva a cabo haciendo

1
y (k+1) = y (k) +
(k1 + 4k2 + k3 ) h
6
function solucion = RK3(xini, h, N, f)
D = length(xini);
x = xini(1);
y = xini(2:D);

solucion = zeros(N+1,D);
solucion(1,:) = [x,y];

for k = 2:N+1

k1 = f(x, y);
k2 = f(x + h/2, y + k1*h/2);
k3 = f(x + h, y - k1*h + 2*k2*h);

x = x + h;
y = y + (k1 + 4*k2 + k3)*h/6;
solucion(k,:) = [x,y];
end;
end

6.6. Formula de integración por el método de Runge-Kutta

clásico de cuarto orden

El más popular de los métodos RK es el de cuarto orden. En el caso del método de Runge-
Kutta de cuarto orden tenemos que calcular cuatro valores k1 , k2 , k3 y k4 utilizando las
siguientes formulas

k1 = f (x(k) , y (k) )

1 1
k2 = f (x(k) + h, y (k) + k1 h)
2 2
6.7. EJEMPLO 127

1 1
k3 = f (x(k) + h, y (k) + k2 h)
2 2

k4 = f (x(k) + h, y (k) + k3 h)

La actualización se lleva a cabo haciendo

1
y (k+1) = y (k) + (k1 + 2k2 + 2k3 + k4 ) h
6

La implementación en Matlab para el método de Runge-Kutta de cuarto orden es:

function solucion = RK4(xini, h, N, f)
D = length(xini);
x = xini(1);
y = xini(2:D);

solucion = zeros(N+1,D);
solucion(1,:) = [x,y];

for k = 2:N+1

k1 = f(x, y);
k2 = f(x + h/2, y + k1*h/2);
k3 = f(x + h/2, y + k2*h/2);
k4 = f(x + h, y + k3*h);

x = x + h;
y = y + (k1 + 2*k2 + 2*k3 + k4)*h/6;
solucion(k,:) = [x,y];
end;
end

6.7. Ejemplo

6.7.1. Una ecuación diferencial sencilla

Consideremos la ecuación diferencial

128 ECUACIONES DIFERENCIALES ORDINARIAS

f (x, y) = −2x3 + 12x2 − 20x + 8.5

La integral la podemos calcular como

Z Z
dy
dx = f (x, y)dx = −0.5x4 + 4x3 − 10x2 + 8.5x
dx

La solución para y para un valor inicial y(0) = 1 lo cual da la siguiente expresión

y(x) = −0.5x4 + 4x3 − 10x2 + 8.5x + 1

Encontrar la solución numérica para esta ecuación diferencial, suponiendo un valor inicial
x(0) = 0, y (0) = 1 y h = 0.5, utilizando los métodos analizados en este capı́tulo.

Solución

La implementación de la función en Matlab es

function dy = ecuacion(x, y)
dy = -2*x^3 + 12*x^2- 20*x + 8.5;
end

En la siguiente tabla se muestran los resultados al aplicar la integración numérica.

Real Euler Heun P. Medio RK RK3 RK4

(k) (k) (k) (k) (k) (k)
k x(k) y (k) yE yH yP M yRK yRK3 yRK4
0 0.0 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1 0.5 3.2188 5.2500 3.4375 3.1094 3.2773 3.2188 3.2188
2 1.0 3.0000 5.8750 3.3750 2.8125 3.1016 3.0000 3.0000
3 1.5 2.2188 5.1250 2.6875 1.9844 2.3477 2.2188 2.2188
4 2.0 2.0000 4.5000 2.5000 1.7500 2.1406 2.0000 2.0000
5 2.5 2.7188 4.7500 3.1875 2.4844 2.8555 2.7188 2.7188
6 3.0 4.0000 5.8750 4.3750 3.8125 4.1172 4.0000 4.0000
7 3.5 4.7188 7.1250 4.9375 4.6094 4.8008 4.7188 4.7188
8 4.0 3.0000 7.0000 3.0000 3.0000 3.0312 3.0000 3.0000

Para evaluar la exactitud de los métodos calculamos el valor absoluto de la diferencia entre
el valor real y el calculado.
6.7. EJEMPLO 129

(k) (k) (k) (k) (k) (k)

k |y (k) − yE | |y (k) − yH | |y (k) − yP M | |y (k) − yRK | |y (k) − yRK3 | |y (k) − yRK4 |
0 0 0 0 0 0 0
1 2.0312 0.2188 0.1094 0.0586 0 0
2 2.8750 0.3750 0.1875 0.1016 0 0
3 2.9062 0.4688 0.2344 0.1289 0 0
4 2.5000 0.5000 0.2500 0.1406 0 0
5 2.0312 0.4688 0.2344 0.1367 0 0
6 1.8750 0.3750 0.1875 0.1172 0 0
7 2.4062 0.2188 0.1094 0.0820 0 0
8 4.0000 0 0 0.0312 0 0
Para reproducir estos valores utilizar el código
function ecuacion_simple()
% Numero de puntos
N = 8;

% Valores Iniciales
xini = 0;
xfin = 4;
y_ini = 1;
h = (xfin-xini)/N;
x = [xini:h:xfin]’;

x1 = Euler([xini, y_ini], h, N, @ecuacion);

x2 = Heun([xini, y_ini], h, N, @ecuacion);
x3 = Punto_Medio([xini, y_ini], h, N, @ecuacion);
x4 = RK([xini, y_ini], h, N, @ecuacion);
x5 = RK3([xini, y_ini], h, N, @ecuacion);
x6 = RK4([xini, y_ini], h, N, @ecuacion);

[x, solucion_real(x)- x1(:,2), x2(:,2), x3(:,2), x4(:,2), x5(:,2), x6(:,2)]

[x, abs(solucion_real(x)-x1(:,2)),...
abs(solucion_real(x)-x2(:,2)),...
abs(solucion_real(x)-x3(:,2)),...
abs(solucion_real(x)-x4(:,2)),...
abs(solucion_real(x)-x5(:,2)),...
abs(solucion_real(x)-x6(:,2))]

function dy = ecuacion(x, y)
130 ECUACIONES DIFERENCIALES ORDINARIAS

dy = -2*x^3 + 12*x^2- 20*x + 8.5;

end

function y = solucion_real(x)
y = -0.5*x.^4 + 4*x.^3 - 10*x.^2 + 8.5*x +1;
end
end

6.7.2. Circuito RL

Para el circuito RL que se muestra en la figura 6.16 calcular el valor de la corriente i(t) en
función del tiempo t, asumiendo que el interruptor se cierra en el tiempo t = 0

Figura 6.16: Circuito RL

Solución
La ecuación diferencial para este problema es

di V − Ri
=
dt L

La ecuación diferencial en Matlab es:

function di = Circuito_RL(t, i)
%Parametros

R = 2.0; % Ohms
L = 0.3; % Henries
V = 10.0; % Volts
di = (V - R *i)/L;
6.7. EJEMPLO 131

end
La solución real de la ecuación diferencial, en función del tiempo, es:

V  
i(t) = 1 − e(−Rt/L)
R
La solución utilizando el Método de Runge-Kutta 4 para un valor inicial t = 0, i(0) = 0 y
un incremento h = 0.2 es:
t i(t)
0.0 0
0.2 3.5391
0.4 4.5732
0.6 4.8753
0.8 4.9636
1.0 4.9894
En la figura 6.17 se muestra la solución gráfica de los métodos estudiados en este capı́tulo.
En negro se da la solución real de la ecuación diferencial y en amarillo la solución del
Runge-Kutta de cuarto orden que es el que mejor se desempeña para el incremento dado.
La peor solución es la calculada con el método de Euler la cual se muestra en azul en esta
figura.

Solucion Circuito RL
7

5
Corriente (Amp.)

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
tiempo (seg.)

Figura 6.17: Circuito RL

El código implementado para este ejemplo es:

132 ECUACIONES DIFERENCIALES ORDINARIAS

function Circuito_RL()

% Numero de puntos
N = 5;

% Valores Iniciales
tini = 0;
tfin = 1;
i_ini = 0;
h = (tfin-tini)/N;

x0 = Solucion_Real([tini, i_ini], h, N);

x1 = Euler([tini, i_ini], h, N, @Circuito_RL);
x2 = RK([tini, i_ini], h, N, @Circuito_RL);
x3 = RK3([tini, i_ini], h, N, @Circuito_RL);
x4 = RK4([tini, i_ini], h, N, @Circuito_RL);

x5 = Heun([tini, i_ini], h, N, @Circuito_RL);

plot(x0(:,1), x0(:,2), ’ko-’, x1(:,1), x1(:,2), ’b .-’, ...

x2(:,1), x2(:,2), ’r*-’, x3(:,1), x3(:,2), ’g--’, ...
x4(:,1), x4(:,2), ’y+-’, x5(:,1), x5(:,2), ’c’);

xlabel(’tiempo (seg.)’);
ylabel(’Corriente (Amp.)’);
title(’Solucion Circuito RL’);

% Solucion exacta de la ecuacion Diferencial

function x =Solucion_Real(xini, h, N)
%Parametros

R = 2.0; % Ohms
L = 0.3; % Henries
V = 10.0; % Volts

t = [xini:h:N*h]’;
Tau = L/R;
i = 5 *(1 - exp(-t/Tau));
6.7. EJEMPLO 133

x = [t, i];
end

% Ecuación diferencial del circuito RL

function di = Circuito_RL(t, i)
%Parametros

R = 2.0; % Ohms
L = 0.3; % Henries
V = 10.0; % Volts
di = (V - R *i)/L;
end
end

6.7.3. Sistema Masa Resorte

El sistema masa resorte como el que se muestra en la figura 6.18, esta constituido por una
masa M , sujeta a un resorte de constante K, el cual esta sujeto a fuerzas de amortigua-
miento C y se aplica una fuerza F . El modelo dinámico para este sistema es:

Figura 6.18: Sistema Masa Resorte

dx2 dx
M 2
+C + Kx = F
dt dt

Note que tenemos una ecuación diferencial de segundo orden la cual plantearemos como
un par de ecuaciones de primer orden como:
134 ECUACIONES DIFERENCIALES ORDINARIAS

dx
=v
dt
dv F − Cv − Kx
=
dt M

La implementación de estas funciones en Matlab es:

function D = Sistema_MR(t, r)
M = 10.0; % Masa en Kg
C = 1; % Coeficiente de Amortiguamiento Dinnamica
K = 5; % Constante del Resorte
F = 0; % Fuerza Aplicada
d = r(1);
v = r(2);

D(1) = v; % derivada del desplazamiento

D(2) = (F - C*v - K*d)/M; % derivada de la velocidad
end
En la figura 6.19 se muestra la solución de la ecuaciones diferenciales utilizando el método
de Runge-Kutta de cuarto orden, para un desplazamiento inicial x(0) = 2, una velocidad
inicial v(0) = 0, de una Masa M = 10, con un coeficiente de amortiguamiento C = 1, una
constante del resorte K = 5 y una fuerza externa F = 0. En la parte superior izquierda
tenemos el desplazamiento en función del tiempo a la derecha la velocidad en función del
tiempo y en la parte inferior el diagrama de fase que muestra la v(t) vs x(t).
La implementación de este ejemplo en Matlab es:
function Sistema_Masa_Resorte(Metodo)

% Numero de puntos
N = 1000;

% Valores Iniciales
tini = 0;
tfin = 40;
d_ini = 2;
v_ini = 0;

h = (tfin-tini)/N;

x = Metodo([tini, d_ini, v_ini], h, N, @Sistema_MR);

6.7. EJEMPLO 135

2 1.5

1
Desplazamiento (m.)

Velocidad (m/s.)
0.5

0 0

−0.5
−1
−1

−2 −1.5
0 10 20 30 40 0 10 20 30 40
tiempo (seg.) tiempo (seg.)

1.5

1
Velocidad (m/s.)

0.5

−0.5

−1

−1.5
−2 −1 0 1 2
Desplazamient (m.)

Figura 6.19: Sistema Masa Resorte

136 ECUACIONES DIFERENCIALES ORDINARIAS

title(’Sistema Masa Resorte’);

subplot(2,2,1)
plot(x(:,1), x(:,2));
xlabel(’tiempo (seg.)’);
ylabel(’Desplazamiento (m.)’);

subplot(2,2,2)
plot(x(:,1), x(:,3));
xlabel(’tiempo (seg.)’);
ylabel(’Velocidad (m/s.)’);

subplot(2,2,3)
plot(x(:,2), x(:,3));
xlabel(’Desplazamient (m.)’);
ylabel(’Velocidad (m/s.)’);

% Ecuación diferencial del circuito RL

function D = Sistema_MR(t, r)
M = 10.0; % Masa en Kg
C = 1; % Coeficiente de Amortiguamiento Dinnamica
K = 5; % Constante del Resorte
F = 0; % Fuerza Aplicada
d = r(1);
v = r(2);

D(1) = v; % derivada del desplazamiento

D(2) = (F - C*v - K*d)/M; % derivada de la velocidad
end
end