Metodos Numericos
Metodos Numericos
29 de julio de 2025
Í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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.3. Ciclos while/end . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.4. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.5. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.6. Ejemplo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.7. Ejemplo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5. Manejo de Matrices y Vectores . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.5.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.5.3. Arreglos Bidimensionales . . . . . . . . . . . . . . . . . . . . . . . . 16
1.6. Estructuras de Programa y Funciones . . . . . . . . . . . . . . . . . . . . . 18
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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.6.4. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6.5. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Modelos 25
ÍNDICE GENERAL
SE No Lineales 35
3.1. Método Gráfico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.2. Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2. Método de Bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.3. Método de la falsa posición (Regula Falsi) . . . . . . . . . . . . . . . . . . . 37
3.3.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.4. Método de iteración de punto fijo . . . . . . . . . . . . . . . . . . . . . . . . 39
3.4.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4.2. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.5. Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.5.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.6. Aplicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.6.1. Cálculo de la Raı́z Cuadrada . . . . . . . . . . . . . . . . . . . . . . 43
3.6.2. Solución de un circuito con un diodo . . . . . . . . . . . . . . . . . . 45
3.6.3. Solución del problema de flujos de Potencia . . . . . . . . . . . . . . 47
3.7. Evaluación de Polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.7.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.7.2. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.7.3. Cálculo de derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.7.4. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.7.5. Calculo de la integral . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.7.6. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.7.7. Ejemplo cálculo de raı́ces de un polinomio . . . . . . . . . . . . . . . 55
3.8. Deflación de polinomios y División de Polinomios . . . . . . . . . . . . . . . 57
3.8.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.8.2. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.8.3. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.8.4. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.9. Solución de un polinomio de orden 2 en la forma general . . . . . . . . . . . 65
3.10. Método de Bairstow para la solución de polinomios . . . . . . . . . . . . . 65
ÍNDICE GENERAL
3.10.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.10.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.10.3. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Sistemas de Ecuaciones 75
4.1. Solución de Sistemas lineales . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.1.1. Método iterativo de Jacobi . . . . . . . . . . . . . . . . . . . . . . . 75
4.1.2. Algoritmo iterativo de Gauss-Seidel . . . . . . . . . . . . . . . . . . 77
4.1.3. Ejemplo matrices dispersas . . . . . . . . . . . . . . . . . . . . . . . 79
4.1.4. Eliminación Gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.1.5. Gauss-Jordan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.2. Métodos para sistemas no lineales . . . . . . . . . . . . . . . . . . . . . . . . 102
4.2.1. Método Gráfico para sistemas no lineales en dos dimensiones . . . . 102
4.2.2. Método de iteración de punto fijo para sistemas . . . . . . . . . . . . 103
4.2.3. Método de Newton-Raphson para sistemas . . . . . . . . . . . . . . 107
4.2.4. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
Optimización 119
5.1. Optimización no-restringida. Método de búsqueda de la sección dorada . . . 119
5.1.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.1.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.2. Optimización no-restringida. Método de Newton . . . . . . . . . . . . . . . 122
5.2.1. Método de Newton en una dimensión . . . . . . . . . . . . . . . . . 122
5.2.2. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
5.2.3. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.2.4. Método de Newton en N dimensiones . . . . . . . . . . . . . . . . . 125
5.2.5. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.2.6. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.2.7. Propiedades del Método de Newton . . . . . . . . . . . . . . . . . . 131
5.2.8. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
5.2.9. Ejemplo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
5.2.10. Problemas de convergencia del Método de Newton . . . . . . . . . . 133
5.2.11. Ejemplo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
5.2.12. Método de Newton Modificado . . . . . . . . . . . . . . . . . . . . . 135
5.3. Metodo Simplex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5.3.1. Forma estándar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5.3.2. Método Simplex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
5.3.3. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.3.4. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.3.5. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
ÍNDICE GENERAL
Tareas 219
11.1. Tarea 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
11.2. Tarea 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
11.3. Tarea 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
11.4. Tarea 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
11.5. Tarea 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
11.6. Tarea 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
11.7. Tarea 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
11.8. Tarea 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
11.9. Tarea 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
[Link] 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
[Link] 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222
[Link] 12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
Introducción a la Programación en 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.
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
1.1. INTRODUCCIÓN (VARIABLES Y OPERADORES) 3
dir.
El comando who (Quien) lista las variables utilizadas en el espacio de trabajo actual.
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:
2. nombres de funciones.
Trigonométricas
4 INTRODUCCIÓN
Exponencial.
Complejos.
1.1. INTRODUCCIÓN (VARIABLES Y OPERADORES) 5
1.1.4. Ejemplo
Calcular el volumen en una esfera.
clear;
r = 2
vol = (4/3)*pi*r^3
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.
6 INTRODUCCIÓN
Operador Simbolo
suma +
resta -
multiplicación *
división /
reciproco \
clear;
c = 3 \ 1
Operador Simbolo
Mayor que >
Menor que <
Mayor igual >=
Menor igual <=
Igual a ==
Diferente de ∼=
clear;
r = -2
if r > 0, (4/3)*pi*r\^3
end
g = 5
if g>3 | g <0, a = 6
end
clear;
a = 2
b = 3
c = 5
if((a==2 | b==3) & c < 5) g=1; end;
1.2. INSTRUCCIONES SECUENCIALES 7
r = -2
if r > 0, (4/3)*pi*r^3
end
Los enunciados lógicos and y or se denotan con & y ∥, pueden ser utilizados para imple-
mentar condicionales de la manera siguiente.
clear;
g = 5
if g>3 | g <0, a = 6
end
pause;
clear;
a = 2
b = 3
c = 5
pause;
clear;
r = 2;
if r > 3 b = 1;
elseif r == 3 b = 2;
else b = 0;
end;
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.
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:
1.4. INSTRUCCIONES DE REPETICIÓN 9
for r=1:5
vol = (4/3)*pi*r^3;
disp([r, vol])
end;
for r=1:5
for s=1:r
vol = (4/3)*pi*(r^3-s^3);
disp([r, s, vol])
end
end
for i=1:6
for j=1:20
if j>2*i, break, end
disp([i, j])
end
end
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.
x = 1;
while x <= 100
x
x = x + 1;
end
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;
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
1.4. INSTRUCCIONES DE REPETICIÓN 11
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
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
for i=1:t
fprintf(’\%s\%s\n’,cad1,cad2)
end
12 INTRODUCCIÓN
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
for i=1:6
x(i) = (i-1)*0.1;
end
x
x(3)
clear;
x = 2: -0.4: -2
pause;
clear;
clear;
x = [1, 2, 3, 4]
y = [4, 3, 2, 1]
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]
14 INTRODUCCIÓN
x = [2, 3]
x = [x, 4]
en el caso de arreglos columna
clear;
y = [2, 3]’
y = [y; 4]
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
1.5. MANEJO DE MATRICES Y VECTORES 15
clear;
v = ’Hola Mundo’
v = ’Hola Mundo’
v = v’
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);
if(max < x(k)); max = x(k); end;
end;
suma = suma/n;
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)
16 INTRODUCCIÓN
suma = a + b
resta = a - b
mult = a .* b
1.5. MANEJO DE MATRICES Y VECTORES 17
div = a ./ b
lo cual es equivalente a
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);
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;
clear;
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]
x = [1, 2, 3]’
b = A*x
para la división tendremos
A = [1, 4; 3, 5]
x = [2, 3]’
b = A\x
des = sqrt(suma/n );
end
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.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.8.
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
La función para modelar el diodo en el circuito es:
function it = diodo(vt)
R = 5;
if vt > 0.7
vd = 0.7;
1.6. ESTRUCTURAS DE PROGRAMA Y FUNCIONES 21
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.
t = [0:0.01:2]; % tiempo
vt = 10*cos(20*t); % Voltaje variante en el tiempo
it = arrayfun(@diodo, vt); % Calculo de la corriente
xlabel(’t (seg.)’);
ylabel(’v(t) y i(t)’);
title (’Solución de un 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
22 INTRODUCCIÓN
v(t) y i(t)
0
−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.)
f (x(k) )
x(k+1) = x(k) −
f ′ (x(k) )
for k = 1: 100
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
.
1.6. ESTRUCTURAS DE PROGRAMA Y FUNCIONES 23
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));
function y = f(x)
y = x-cos(x);
24 INTRODUCCIÓN
Modelos, computadoras y errores
1 ii 1 1
f (x) = f (a)+f i (a)(x−a)+ f (a)(x−a)2 + f iii (a)(x−a)3 +· · ·+ f n (a)(x−a)n (2.1)
2! 3! n!
N
X f n (a)
f (x) = (x − a)n
n!
n=0
dn f (x)
donde f 0 (a) = f (x)|x=a y f n (a) = dxn x=a
25
26 MODELOS
1 ii (i) 2 1 1
f (x(i+1) ) = f (x(i) ) + f i (x(i) )h + f (x )h + f iii (x(i) )h3 + · · · + f n (x(i) )hn (2.2)
2! 3! n!
2.2.2. Ejemplo 1
Dada la función f (x) = ex hacer la implementación utilizando la serie de Taylor dada por
(2.1) considerando que a = 0.
f (x) = f (x) = ex
df (x)
f i (x) = = ex
dx
d2 f (x)
f ii (x) = = ex
dx2
d3 f (x)
f iii (x) = = ex
dx3
d4 f (x)
f iv (x) = = ex
dx4
1 ii 1 1
ex ≈ f (a) + f i (a)(x − a) +
f (a)(x − a)2 + f iii (a)(x − a)3 + · · · + f n (a)
2! 3! n!
1 1 1
ex ≈ 1 + 1 × (x − 0) + × 1 × (x − 0)2 + × 1 × (x − 0)3 + × 1 × (x − 0)4 + · · ·
2! 3! 4!
x2 x3 x4 x5 x6 x7
ex ≈ 1 + x + + + + + + + ···
2! 3! 4! 5! 6! 7!
N
X xk
ex ≈
k!
k=0
function y = exponencial(x)
y =0;
for k= 0:40
y = y + x.^k/factorial(k);
end;
Ver exponencial.m
2.2.3. Ejemplo 2
Dada la función f (x) = cos(x) hacer la implementación utilizando la serie de Taylor dada
por (2.1) considerando que a = 0.
Sustituyendo para x = a = 0
28 MODELOS
x2 x4 x6 x8 x10
cos(x) ≈ 1 − + − + −
2! 4! 6! 8! 10!
En forma compacta podemos representar la serie como:
N
X (−1)k x2k
cos(x) ≈
(2k)!
k=0
2.2.4. Ejemplo 3
Utilizar los términos de la serie de Taylor para n = 0 hasta 6, para aproximar la función
f (x) = cos(x) en x(i+1) = π/3 y como condición inicial sus derivadas en x(i) = π/4.
2.2. SERIE DE TAYLOR Y ERRORES DE TRUNCAMIENTO 29
Para nuestra implementación tenemos que el incremento h = π/3 − π/4 = π/12. La apro-
ximación en serie de Taylor utilizando (2.2) es:
1 ii (i) 2 1 1
f (x(i+1) ) = f (x(i) ) + f i (x(i) )h + f (x )h + f iii (x(i) )h3 + · · · + f n (x(i) )hn
2! 3! n!
π π π 1 π π 2 1 π π 3 1 π π n
f (x(i+1) ) ≈ f +f i + f ii + f iii +· · ·+ f n
4 4 12 2! 4 12 3! 4 12 n! 4 12
N
X 1 n π π n
f (x(i+1) ) ≈ f
n! 4 12
n=0
N =0
f (x(i+1) ) ≈ 0.7071
N =1
π
f (x(i+1) ) ≈ 0.7071 + (−0.7071) × = 0.5219
12
N =2
π 1 π 2
f (x(i+1) ) ≈ 0.7071 + (−0.7071) × + × (−0.7071) = 0.4977
12 2 12
2.2.5. Ejemplo 4
Hacer la implementación en Serie de Taylor para la función f (x) = sen(x).
De acuerdo con lo anterior la serie de Taylor para el seno esta dada como
x3 x5 x7 x9
sen(x) ≈ x − + − + + ···
3! 5! 7! 9!
N
X (−1)k x2k+1
sen(x) ≈
(2k + 1)!
k=0
for k= 0:N
y = y + (-1)^k*x.^(2*k+1)/factorial(2*k+1);
end
end
ver seno.m
′ 1 ′′ 1
v(x(i+1) ) = v(x(i) )+v (x(i) )(x(i+1) −x(i) )+ v (x(i) )(x(i+1) −x(i) )2 +· · ·+ v n (xi )(x(i+1) −x(i) )n
2! n!
Ahora, truncando la serie después del término con la primera derivada, se obtiene:
′
v(x(i+1) ) = v(x(i) ) + v (x(i) )(x(i+1) − x(i) ) + R1
′ v(x(i+1) ) − v(x(i) ) R1
v (x(i) ) = (i+1) (i)
− (i+1)
(x −x ) (x − x(i) )
N = length(x);
dv = zeros(N-1);
for k=1:N-1
dy(k) = (y(k+1) - y(k))/(x(k+1) - x(k));
32 MODELOS
end;
xnva = x(1:N-1);
end
Ejemplo
Vamos a calcular la derivada de la función y = cos(x)y comparar con el valor exacto de
la derivada y ′ = −sen(x). Para esto consideramos que 0 ≤ x ≤ 2π y que el incremento
es h = 0.01. El siguiente código muestra la manera de utilizar la función derivada y en
la Figura 2.4 se muestra la gráfica de la derivada real y la aproximada con la serie de
Taylor.
h = 0.01;
x = [0:h:2*pi];
y = cos(x);
dy = sin(x);
[dx, dy] = derivada(x, y);
plot(x, -sin(x), dx, dy)
grid on
ver derivada.m
2.2. SERIE DE TAYLOR Y ERRORES DE TRUNCAMIENTO 33
Figura 2.4: Gráfica mostrando la derivada real del coseno contra la estimación utilizando
serie de Taylor
34 MODELOS
Solución de Ecuaciones no-lineales
3.1.1. Ejemplo
Utilice el método gráfico para observar algunas de las raı́ces de la función f (x) = sen10x +
cos3x, en el intervalo [0, 5]
1.5
0.5
f(x)
−0.5
−1
−1.5
−2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
x
En la Figura 3.5 se muestran los cruces por cero de la función en el intervalo dado.
35
36 SE NO LINEALES
3.1.2. Implementación
La implementación en Matlab es
function Metodo_Grafico(funcion, x)
f = funcion;
y = f(x);
plot(x, y);
xlabel(’x’);
ylabel(’f(x)’)
grid on
Ver Metodo Grafico.m. La implementación de la función
function y = f1(x)
y = sin (10*x) + cos(3*x);
La corrida del método es
Metodo_Grafico(@f1, 0:0.01:5)
if(f(mitad) == 0)
r = mitad;
return;
3.3. MÉTODO DE LA FALSA POSICIÓN (REGULA FALSI) 37
end;
if f(inicio)*f(mitad) < 0
fin = mitad;
else
inicio = mitad;
end;
end;
r= mitad;
3.2.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.000000 0.500000 1.000000 -1.000000 -0.377583 0.459698
1 0.500000 0.750000 1.000000 -0.377583 0.018311 0.459698
2 0.500000 0.625000 0.750000 -0.377583 -0.185963 0.018311
3 0.625000 0.687500 0.750000 -0.185963 -0.085335 0.018311
4 0.687500 0.718750 0.750000 -0.085335 -0.033879 0.018311
5 0.718750 0.734375 0.750000 -0.033879 -0.007875 0.018311
6 0.734375 0.742188 0.750000 -0.007875 0.005196 0.018311
7 0.734375 0.738281 0.742188 -0.007875 -0.001345 0.005196
8 0.738281 0.740234 0.742188 -0.001345 0.001924 0.005196
9 0.738281 0.739258 0.740234 -0.001345 0.000289 0.001924
10 0.738281 0.738770 0.739258 -0.001345 -0.000528 0.000289
11 0.738770 0.739014 0.739258 -0.000528 -0.000120 0.000289
12 0.739014 0.739136 0.739258 -0.000120 0.000085 0.000289
13 0.739014 0.739075 0.739136 -0.000120 -0.000017 0.000085
Para ver los resultados correr
Biseccion(@f2, 0, 1)
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.6 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.
function r = Regula_Falsi(f, a, b)
if(f(c) == 0)
r = c;
return;
end;
if f(a)*f(c) < 0
b = c;
else
a = c;
end;
end;
r= c;
3.3.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.4597
2 0.7363 0.7389 1.0000 -0.0047 -0.0002 0.4597
3 0.7389 0.7391 1.0000 -0.0002 0.0000 0.4597
4 0.7391 0.7391 1.0000 0.0000 0.0000 0.4597
5 0.7391 0.7391 1.0000 0.0000 0.0000 0.4597
6 0.7391 0.7391 1.0000 0.0000 0.0000 0.4597
la ecuación f (x) = 0 de tal modo que x quede del lado izquierdo de la ecuación
f (x) ≡ x − g(x) = 0
Lo que es igual a
x = g(x)
x2 − 2x + 3 = 0
x2 + 3
x=
2
x = g(x)
x = sen(x) + x
Las iteraciones del método dado un punto inicial x(0) , se pueden llevar a cabo calculando
x(k+1)
x(k+1) = g(x(k) )
while 1
x2 = g(x1)
error = abs((x2-x1)/x2);
if(error < 0.0001) break
else x1 = x2;
end;
end;
r = x2;
3.4. MÉTODO DE ITERACIÓN DE PUNTO FIJO 41
3.4.1. Ejemplo
Utilizando la iteración de punto fijo calcular las raı́ces f (x) = e−x − x, para un valor inicial
x(0) = 0.
De acuerdo con lo revisado la función g(x) = e−x y las iteraciones hasta convergencia
son:
k x(k) error(x(k) )
1 1.0000 1.0000
2 0.3679 1.7183
3 0.6922 0.4685
4 0.5005 0.3831
5 0.6062 0.1745
6 0.5454 0.1116
7 0.5796 0.0590
8 0.5601 0.0348
9 0.5711 0.0193
.. .. ..
. . .
19 0.5672 6.7008e-05
Para hacer la ejecución correr:
Punto_Fijo(@g1, 0)
con
function y = g4(x)
y = exp(-x);
3.4.2. Ejemplo
Calcular utilizando iteración de punto fijo, un cero de la función f (x) = sen(10x) +
cos(3x)
Para llevar a cabo la iteración de punto fijo hacemos
x = sen(10x) + cos(3x) + x
La implementación en Matlab es
function y = g1(x)
y = sin(10*x)+cos(3*x) +x;
42 SE NO LINEALES
f (x(k) )
x(k+1) = x(k) −
f ′ (x(k) )
while 1
x2 = x1 - f(x1)/df(x1)
if abs((x2-x1)/x2) < 0.0001 break;
else x1 = x2;
end;
end;
r = x2;
end
3.5.1. Ejemplo
Calcular los ceros de la función f (x) = x − cos(x) utilizando el algoritmo de Newton
Raphson con x0 = 0.
3.6. APLICACIONES 43
k x(k)
0 0.0000
1 1.0000
2 0.7504
3 0.7391
4 0.7391
5 0.7391
6 0.7391
7 0.7391
En estas Figura 3.7 se muestra tres iteraciones del método. En la figura de la izquierda
mostramos la lı́nea recta que es tangente a la función f (x) (en negro) en x(0) = 0, note
que la lı́nea recta cruza el eje x en x(1) = 1. En el centro la tenemos la segunda iteración
tomando como valor inicial x(1) = 1 y solución x(2) = 0.7504. Finalmente a la izquierda
tenemos la linea recta y la solución para x(3) = 0.7391. Note como poco a poco se acerca
a la solución.
Newton_Raphson(@f2, @df2, 0)
3.6. Aplicaciones
3.6.1. Cálculo de la Raı́z Cuadrada
Hacer un algoritmo iterativo que permita hacer el cálculo de la raı́z cuadrada de A.
2
√ caso nuestra función a resolver es f (x) = x − A. La solución cuando f (x) = 0
Para este
es x = A. Para nuestro calculo hacemos
44 SE NO LINEALES
f (x) = x2 − A
f ′ (x) = 2x
2
f (x(k) ) x(k) − A
x(k+1) = x(k) − = x(k)
−
f ′ (x(k) ) 2x(k)
2
(k+1) x(k) + A
x =
2x(k)
La función en Matlab queda:
function y = Raiz_Cuadrada(A)
r = A;
while abs(r*r - A) > 0.00001
r = (r*r+A)/(2*r);
end;
end
3.6. APLICACIONES 45
La ecuación que modela al diodo esta dada por (3.3) y la curva se muestra en la figura
3.9:
log(Id (t)/Is + 1) ∗ Vt Si Id (t) ≥ 0
Vd (t) =
V sino Id(t) < 0
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
%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
3.6. APLICACIONES 47
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
−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
V − Ri − VL = 0
V i − Ri2 − PL = 0
V = 20; % Volts
R = 2; % ohms
P = 12; % Watts
function d = fp(I, p)
3.7. EVALUACIÓN DE POLINOMIOS 49
V = p(1);
R = p(2);
PL = p(3);
d = V*I - R*I*I - PL;
end
>> Simula_Flujos_Potencia
La solución es 3.000000
f3 (x) = a0 + a1 x + a2 x2 + a3 x3
En este caso los coeficientes del polinomio los podemos poner en un vector de datos dado co-
mo a = [a0 , a1 , a2 , a3 ], y realizar la evaluación polinomial utilizando ciclos. Adicionalmente
necesitamos de una función que nos permita calcular la potencias de x.
50 SE NO LINEALES
Manipulando el polinomio, podemos dar una forma eficiente de realizar la evaluación. Ası́
el mismo polinomio, lo podemos representar como:
function y = Evalua_Polinomio(x, a)
n = length(a);
p = 0.0;
for i=n:-1: 1
p = p.*x + a(i);
end
y = p;
3.7.1. Ejemplo
Consideremos el siguiente polinomio f2 (x) = 1 + 2x + 3x2 , el cual deseamos evaluar en
x = 20. Para implementarlo hacemos
i p = p ∗ x + ai p=0
2 p = 0 ∗ 20 + 3 p=3
1 p = 3 ∗ 20 + 2 p = 62
0 p = 62 ∗ 20 + 1 p = 1241
3.7.2. Ejemplo
Hacer la evaluación del polinomio f5 (x) = x5 − 3.5x4 + 2.75x3 + 2.125x2 − 3.875x + 1.25
en x = 2.5
3.7. EVALUACIÓN DE POLINOMIOS 51
i p = p ∗ x + ai p=0
5 p = 0.000000 ∗ 2.500000 + 1.000000 1.000000
4 p = 1.000000 ∗ 2.500000 − 3.500000 -1.000000
3 p = −1.000000 ∗ 2.500000 + 2.750000 0.250000
2 p = 0.250000 ∗ 2.500000 + 2.125000 2.750000
1 p = 2.750000 ∗ 2.500000 − 3.875000 3.000000
0 p = 3.000000 ∗ 2.500000 + 1.250000 8.750000
Los coeficientes a del polinomio original pueden tomarse para calcular los coeficiente del
nuevo polinomio, para esto tomamos en cuenta lo siguiente
a = [a0 , a1 , a2 , a3 , a4 , · · · , an−1 , an ]
Con esto tenemos que dado a los coeficientes de del polinomio de derivadas pueden ser cal-
culados y la evaluación de la derivada se reduce a una simple evaluación polinomial.
Ası́ por ejemplo considerando el polinomio f3 (x) = 1 + 2x + 3x2 + 4x3 , tenemos que
a = [1, 2, 3, 4] y los coeficiente de f3′ (x) son a′ = [2, 6, 12]. La evaluación del polinomio
resultante lo podemos hacer para x = 20 utilizando:
i p = p ∗ x + a′i p
2 0*20 +12 12
1 12*20 + 6 246
0 246*20 + 2 4922
52 SE NO LINEALES
3.7.4. Ejemplo
Dado el polinomio de la serie de Taylor correspondiente a la función f (x) = seno(x), hacer
la evaluación de la derivada en x = 0.5.
El polinomio de la serie de Taylor para seno es:
x3 x5 x7 x9
f (x) = x − + − + + ···
3! 5! 7! 9!
Los coeficientes de este polinomio son a = [0, 1, 0, −1/3!, 0, 1/5!, 0, −1/7!, 0, 1/9!], por lo
tanto los coeficientes para la derivada son a∗ = [1, 0, −3/3!, 0, 5/5!, 0, −7/7!, 0, 9/9!]
i p = p ∗ x + a∗i p
8 p = 0.000000 ∗ 0.500000 + 0.000025 0.000025
7 p = 0.000025 ∗ 0.500000 + 0.000000 0.000012
6 p = 0.000012 ∗ 0.500000 − 0.001389 -0.001383
5 p = −0.001383 ∗ 0.500000 + 0.000000 -0.000691
4 p = −0.000691 ∗ 0.500000 + 0.041667 0.041321
3 p = 0.041321 ∗ 0.500000 + 0.000000 0.020660
2 p = 0.020660 ∗ 0.500000 − 0.500000 -0.489670
1 p = −0.489670 ∗ 0.500000 + 0.000000 -0.244835
0 p = −0.244835 ∗ 0.500000 + 1.000000 0.877583
La implementación en Matlab es
function y = Evalua_Derivada(x, a)
N = length(a);
b = a(2:N).*[1:N-1];
y = Evalua_Polinomio(x, b);
end
Z
a1 2 a2 3 a3 4 an−1 n an n+1
f (x)dx = CT E + a0 x + x + x + x ··· + x + x
2 3 4 n n+1
Note que los coeficiente del nuevo polinomio a∗ pueden calcularse a partir de los coeficientes
a de polinomio original :
a = [a0 , a1 , a2 , a3 , a4 , · · · , an−1 , an ]
Ası́ la evaluación de la integral del polinomio se reduce a calcular los nuevos coeficientes
del nuevo polinomio y hacer la evaluación del nuevo polinomio con coeficientes a∗ .
En el caso de tener una integral definida debemos notar:
Z x1 x1
a1 2 a2 3 a3 4 an−1 n an n+1
f (x)dx = CT E + a0 x + x + x + x · · · + x + x
x0 2 3 4 n n+1 x0
3.7.6. Ejemplo
Para la serie de Taylor de la función seno hacer la evaluación de la integral:
Z π/2
sin(x)dx
0
x3 x5 x7 x9
f (x) = x − + − + + ···
3! 5! 7! 9!
54 SE NO LINEALES
Los coeficientes de este polinomio son a = [0, 1, 0, −1/3!, 0, 1/5!, 0, −1/7!, 0, 1/9!], por lo
tanto los coeficientes para la integral serán :
−1 −1
a∗ = [0, 0, 12 , 0, 4×3! 1
, 0, 6×5! , 0, 8×7! 1
, 0, 10×9! ]
i p = p ∗ x + a∗i p
10 p = 0.000000 * 0.000000 + 0.000000 0.000000
9 p = 0.000000 * 0.000000 + 0.000000 0.000000
8 p = 0.000000 * 0.000000 - 0.000025 -0.000025
7 p = -0.000025 * 0.000000 + 0.000000 0.000000
6 p = 0.000000 * 0.000000 + 0.001389 0.001389
5 p = 0.001389 * 0.000000 + 0.000000 0.000000
4 p = 0.000000 * 0.000000 - 0.041667 -0.041667
3 p = -0.041667 * 0.000000 + 0.000000 0.000000
2 p = 0.000000 * 0.000000 + 0.500000 0.500000
1 p = 0.500000 * 0.000000 + 0.000000 0.000000
0 p = 0.000000 * 0.000000 + 0.000000 0.000000
Con esto tenemos que el valor de la integral es igual a 1 − 0 = 1.
La solución en Matlab queda
>> a = [0, 1, 0, -1/factorial(3), 0, 1/factorial(5),
0, -1/factorial(7), 0, 1/factorial(9)]
3.7. EVALUACIÓN DE POLINOMIOS 55
ans =
1.0000
fn (x) = a0 + a1 x + a2 x2 + ... + an xn .
f (−s(k) )
s(k+1) = s(k) +
f ′ (−s(k) )
Para hacer la solución más eficiente evaluaremos la función utilizando evaluación polinomial
tanto para la función como para la derivada.
56 SE NO LINEALES
function sn = factor_NR(s0, a)
% factor_NR(s0, a) Calcula el factor (x+s) de un polinomio con
% coeficiente a dado un valor inicial s0 mediante el
% metodo de Newton-Raphson
Como ejemplo numérico calculamos la raı́z del polinomio f6 (x) = x6 −1 con un valor inicial
s(0) = 10
3.8. DEFLACIÓN DE POLINOMIOS Y DIVISIÓN DE POLINOMIOS 57
999999.000000
s1 = 10.000000 + −600000.000000 = 8.333335
334897.378558
s2 = 8.333335 + −241126.784337 = 6.944450
112156.191251
s3 = 6.944450 + −96903.735989 = 5.787052
37560.618299
s4 = 5.787052 + −38943.785362 = 4.822569
12578.711852
s5 = 4.822569 + −15651.050576 = 4.018871
4212.321940
s6 = 4.018871 + −6290.306217 = 3.349218
1410.434918
s7 = 3.349218 + −2528.533032 = 2.791411
472.088718
s8 = 2.791411 + −1016.880869 = 2.327159
157.838757
s9 = 2.327159 + −409.526161 = 1.941741
52.597923
s10 = 1.941741 + −165.618133 = 1.624156
17.355481
s11 = 1.624156 + −67.809321 = 1.368210
5.560199
s12 = 1.368210 + −28.768380 = 1.174936
1.630779
s13 = 1.174936 + −13.434500 = 1.053548
0.367497
s14 = 1.053548 + −7.787952 = 1.006360
0.038774
s15 = 1.006360 + −6.193251 = 1.000100
0.000598
s16 = 1.000100 + −6.002990 = 1.000000
fn (x) = a0 + a1 x + a2 x2 + ... + an xn .
La división sintética, es una manera de hacer lo mismo, pero de forma compacta. Ası́
tenemos para el ejemplo anterior que podemos hacer:
x+s a3 a2 a1 a0
−a3 s −(a2 − a3 s)s −(a1 − (a2 − sa3 )s)s
a3 a2 − a 3 s a1 − (a2 − a3 s)s a0 − (a1 − (a2 − a3 s)s)s
bn = an
bn−1 = an−1 − an s = an−1 − bn s
bn−2 = an−2 − (an−1 − an s)s = an−2 − bn−1 s
.. .
. = ..
b2 = a2 − b3 s
b1 = a1 − b2 s
b0 = a0 − b1 s
En general
bn = an
bk = ak − bk+1 s
∀k ∈ [n − 1, n − 2, · · · , 2, 1, 0]
n = length(a);
b = zeros(1,n);
b(n) = a(n);
disp(’residuo ’);
disp(b(1));
y = b(2:n);
Y para realizar la ejecución hacemos
Division_Sintetica([-120, -46, 79, -3, -7, 1], -4)
60 SE NO LINEALES
3.8.1. Ejemplo
Dado el polinomio f5 (x) = x5 − 7x4 –3x3 + 79x2 –46x–120 encontrar la división sintética con
el monomio (x − 4).
El resultado es
30 19 -15 -3 1
El polinomio resultante, en este caso, es f4 (x) = x4 − 3x3 –15x2 + 19x + 30. Note que el
residuo es cero y en este caso tenemos una deflación polinomial.
En algunos casos en deseable hacer la división sintética de la función fn (x) entre un po-
linomio de segundo orden de la forma f2 (x) = x2 + rx + s. Aplicando la técnica anterior
podemos realizar la división fn (x)/f2 (x) como
b5 = a5
b4 = a4 − a5 r = a4 − rb5
b3 = a3 − sa5 − (a4 − a5 r)r = a3 − sb5 − rb4
b2 = a2 − sb4 − rb3
b1 = a1 − sb3 − rb2
b0 = a0 − sb2 − rb1
3.8. DEFLACIÓN DE POLINOMIOS Y DIVISIÓN DE POLINOMIOS 61
En general para un polinomio de grado fn (x) tenemos que podemos hacer la división
sintética aplicando la siguiente recurrencia
bn = an
bn−1 = an−1 − rbn
bk = ak − rbk+1 − sbk+2
∀k ∈ {n − 2, n − 3, · · · 1}
R = b1 (x + r) + b0
n = length(a);
b=zeros(n,1);
b(n) = a(n);
b(n-1) = a(n-1) - r*b(n);
c = b(3:n);
r = [b(2); b(1)+r*b(2)];
end
3.8.2. Ejemplo
Hacer la división sintética del polinomio f3 (x) = x3 − 1, entre el polinomio f2 (x) = x2 +
3x + 4
62 SE NO LINEALES
b3 = a3 = 1.0
b2 = a2 − rb3 = 0 − (1)(3) = −3
b1 = a1 − rb2 − sb3 = 0 − (3)(−3) − (4)(1) = 5
b0 = a0 − rb1 − sb2 = −1 − (3)(5) − (4)(−3) = −4
R = b1 (x + r) + b0 = 5(x + 3) − 4 = 5x + 11
b =
-3
1
r =
11
5
3.8.3. Ejemplo
Hacer la división sintética del polinomio f5 (x) = x5 −3.5x4 +2.75x3 +2.125x2 −3.875x+1.25,
entre el polinomio f2 (x) = x2 + x − 2
De acuerdo con la sucesión
b5 = a5 = 1.0
b4 = a4 − rb5 = −3.5 − (1)(1) = −4.5
b3 = a3 − rb4 − sb5 = 2.75 − (1)(−4.5) − (−2)(1) = 2.75 + 6.5 = 9.25
b2 = a2 − rb3 − sb4 = 2.125 − (1)(9.25) − (−2)(−4.5) = −16.125
b1 = a1 − rb2 − sb3 = −3.875 − (1)(−16.125) − (−2)(9.25) = 30.75
b0 = a0 − rb1 − sb2 = 1.25 − (1)(30.75) − (−2)(−16.125) = −61.75
3.8. DEFLACIÓN DE POLINOMIOS Y DIVISIÓN DE POLINOMIOS 63
La siguiente tabla muestra el resumen de los coeficientes calculados. Note que el residuo
no es cero.
1 -3.5 2.75 2.125 -3.875 1.25
x2 + x − 2 -1 6.5 -18.25 34.625 -63.00
1 -4.5 9.25 -16.125 30.75 -61.75
b =
-16.1250
9.2500
-4.5000
1.0000
r =
-31.0000
30.7500
3.8.4. Ejemplo
Repetir el problema anterior pero suponiendo que f2 (x) = x2 − x − 2
b5 = a5 = 1.0
b4 = a4 − rb5 = −3.5 − (−1)(1) = −2.5
b3 = a3 − rb4 − sb5 = 2.75 − (−1)(−2.5) − (−2)(1) = 2.25
b2 = a2 − rb3 − sb4 = 2.125 − (−1)(2.25) − (−2)(−2.5) = −0.625
b1 = a1 − rb2 − sb3 = −3.875 − (−1)(−0.625) − (−2)(2.25) = 0
b0 = a0 − rb1 − sb2 = 1.25 − (−1)(0) − (−2)(−0.625) = 0
con un residuo
R = 0 × (x − 1) + 0 = 0
La corrida en Matlab es
b =
-0.6250
2.2500
-2.5000
1.0000
r =
0
0
f2 (n) = a0 + a1 x + a2 x2
p
−a1 ± a21 − 4a2 a0
x=
2a2
Para calcular la división de polinomios, hacemos uso de la división sintética. Ası́ dado
bn = an
bn−1 = an−1 − rbn
bi = ai − rbi+1 − sbi+2
Una manera de determinar los valores de r y s que hacen cero el residuo es utilizar el
Método de Newton-Raphson. Para ello necesitamos una aproximación lineal de b1 y b0
respecto a r y s la cual calculamos utilizando la serie de Taylor
donde los valores de rk y sk están dados y calculamos los nuevos valores rk+1 y sk+1 que
hacen a b1 (rk+1 , sk+1 ) y b0 (rk+1 , sk+1 ) igual a cero. El sistema de ecuaciones que tenemos
que resolver es:
cn = bn
cn−1 = bn−1 − rcn
ci = bi − rci+1 − sci+2
donde
∂b0 (rk , sk )
= c1 (rk , sk )
∂r
∂b1 (rk , sk )
= c2 (rk , sk )
∂r
∂b0 (rk , sk )
= c2 (rk , sk )
∂s
∂b1 (rk , sk )
= c3 (rk , sk )
∂s
El código en Matlab para implementar la división sintética y calculo de los coeficientes del
polinomio de derivadas es
function [b, c] = Division_Sintetica_Derivadas(a, r, s)
n = length(a);
b=zeros(n,1);
c=zeros(n,1);
b(n) = a(n);
b(n-1) = a(n-1) - r*b(n);
c(n) = b(n);
68 SE NO LINEALES
c2 (rk , sk ) c3 (rk , sk ) drk −b1 (rk , sk )
=
c1 (rk , sk ) c2 (rk , sk ) dsk −b0 (rk , sk )
−1
rk+1 rk c2 (rk , sk ) c3 (rk , sk ) −b1 (rk , sk )
= −
sk+1 sk c1 (rk , sk ) c2 (rk , sk ) −b0 (rk , sk )
N = length(a);
res =[r0, s0]’;
for k=1:100
[b, c] = Division_Sintetica_Derivadas(a, res(1), res(2));
r = res(1);
s = res(2);
c = b(3:N);
3.10.1. Ejemplo 1
Dado el polinomio f5 (x) = x5 − 3.5x4 + 2.75x3 + 2.125x2 − 3.875x + 1.25, determinar los
valores de r y s que hacen el resido igual a cero. Considere r0 = 1 y s0 = −2.
3.10. MÉTODO DE BAIRSTOW PARA LA SOLUCIÓN DE POLINOMIOS 69
Solución.
Iteración 1
La división sintética con el polinomio f2 (x) = x2 + x − 2.0 da como resultado
El polinomio de derivadas es
−1
r1 1 −43.875 16.750 −30.7500 −1.7637
= − =
s1 −2 108.125 −43.875 61.7500 −7.4034
Iteración 2
La división sintética con el polinomio f2 (x) = x2 −1.7637x−7.4034 da como resultado
−1
r2 −1.7637 27.62800 14.5427 −51.75640 −1.7164
= − =
s2 −7.4034 208.1484 27.6280 −105.68578 −3.9343
Iteración 3
70 SE NO LINEALES
−1
r3 −1.7164 13.8350 7.4418 −12.6547 −1.5997
= − =
s3 −3.9343 65.6792 13.8350 −28.1881 −2.4507
En resumen
k r s b1 b0
0 1.0000 - 2.0000 30.7500 -61.7500
1 -1.7637 -7.4034 51.7564 105.6858
2 -1.7164 -3.9343 12.6547 28.1881
3 -1.5997 -2.4507 2.8996 8.1547
4 -1.3335 -2.1867 0.7601 2.5222
5 -1.1183 -2.1130 0.2719 0.6077
6 -1.0271 -2.0232 0.0431 0.1119
7 -1.0017 -2.0015 0.0028 0.0063
8 -1.0000 -2.0000 0.0000 0.0000
9 -1.0000 -2.0000 0.0000 0.0000
10 -1.0000 -2.0000 0.0000 0.0000
La solución es:
c =
-0.6250
2.2500
-2.5000
1.0000
r =
-1.0000
s =
-2.0000
raices = [];
N = length(b);
while(N > 3)
[b r s] = Itera_Bairstow(b, r0, s0);
[x1, x2] = Formula_General(1, r, s);
raices = [raices, x1, x2];
N = N-2;
end;
if length(b) == 3
[x1, x2] = Formula_General(b(3), b(2), b(1));
raices = [raices, x1, x2];
end;
72 SE NO LINEALES
if length(b) == 2
raices = [raices, -b(1)/b(2)];
end;
3.10.2. Ejemplo 2
Dado el polinomio f5 (x) = x5 − 3.5x4 + 2.75x3 + 2.125x2 − 3.875x + 1.25, calcular la raı́ces
utilizando el Método de Bairstow con valores iniciales r0 = 1 y s0 = −2
Iteración 1
Calculamos la deflación Polinomial
x1 = 2.0000
x2 = −1.0000
Iteración 2
x3 = 1 + j0.5
x4 = 1 − j0.5
Iteración 3
Finalmente
f1 (x) = x − 0.5
ans =
3.10.3. Ejemplo 3
Dado el polinomio f5 (x) = x5 − 3.5x4 + 2.75x3 + 2.125x2 − 3.875x + 1.25, determinar las
raı́ces de este polinomio. Considere r0 = −1 y s0 = −1.
Iteración 1
x1 = 0.5
x2 = −1.0
Iteración 2
x3 = 1.0 + j0.5
x4 = −1.0 − j0.5
Iteración 3
74 SE NO LINEALES
f1 (x) = (x − 2)
x5 = 2;
Todas la raı́ces de f5 (x) son x = [0.5, −1.0, (1.0 + j0.5), (1 − j0.5), 2]. Para correr en Matlab
dar
>> Bairstow([1.25, -3.875, 2.125, 2.75, -3.5, 1], -1, -1)
ans =
En forma compacta
Ax = b
En las siguientes subsecciones se analizarán algunos de los métodos para resolver el 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
75
76 SISTEMAS DE ECUACIONES
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
(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
PN (t)
(t+1) bk – l=1,l̸=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
4.1. SOLUCIÓN DE SISTEMAS LINEALES 77
(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;
78 SISTEMAS DE ECUACIONES
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
ans =
2.0000
4.1. SOLUCIÓN DE SISTEMAS LINEALES 79
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
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;
A(1,4) = -1;
A(1,5) = 1;
80 SISTEMAS DE ECUACIONES
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)
ans =
0.1520
0.4620
0.7120
0.9620
1.6160
ans =
0.1520
0.4620
0.7120
0.9620
1.6160
4.1. SOLUCIÓN DE SISTEMAS LINEALES 81
b1 − a1,2 x2 − a1,3 x3
x1 =
a1,1
b1 − a1,2 x2 − a1,3 x3
a2,1 + a2,2 x2 + a2,3 x3 = b2
a1,1
a2,1 ∗ a1,2 a1,3 b1
a2,2 − x2 + a2,3 − a2,1 x3 = b2 − a2,1
a1,1 a1,1 a1,1
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
a1,1 a1,2 a1,3 x1 b1
0 a′2,2 a′2,3 x2 = b′2
0 a′3,2 a′3,3 x3 b′3
82 SISTEMAS DE ECUACIONES
a1,1 a1,2 a1,3 x1 b1
0 a′2,2 a′2,3 x2 = b′2
0 0 a′′3,3 x3 b′′3
a1,1 a1,2 a1,3 a1,4 x1 b1
0 a2,2 a2,3 a2,4 x2 b2
=
0 0 a3,3 a3,4 x3 b3
0 0 0 a4,4 x4 b4
Podemos resolverlo empezando con la ecuación 4 cuya solución es la más simple, para luego
solucionar 3 y ası́ sucesivamente. Este esquema de solución queda
b4
x4 =
a4,4
b3 − a3,4 x4
x3 =
a3,3
b2 − a2,3 x3 − a2,4 x4
x2 =
a2,2
b1 − a1,2 x2 − a1,3 x3 − a1,4 x4
x1 =
a1,1
PN
bk − i=k+1 ak,i xi
xk =
ak,k
function x = Eliminacion_Gaussiana(A, b)
N = length(b);
x = zeros(N,1);
for k =[Link]N-1
for n = k+[Link]N
b(n) = b(n) - A(n,k)*b(k)/A(k,k);
for m=N:-1:k
A(n,m) = A(n,m) - A(n,k)*A(k,m)/A(k,k);
end;
end;
end;
for k=N:-1:1
suma = 0;
for m=k+[Link]N
suma = suma +A(k,m)*x(m);
end;
x(k) = (b(k)-suma)/A(k,k);
end;
que corresponde a un sistema triangular superior que podemos solucionar utilizando susti-
tución hacia atrás.
Ejemplo 1
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
84 SISTEMAS DE ECUACIONES
−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
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
4.1. SOLUCIÓN DE SISTEMAS LINEALES 85
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
b4
x4 =
a4,4
2.0397
x4 = = 0.4037
5.0529
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]’)
ans =
-0.0508
0.2372
0.6707
0.4037
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
Segundo paso
5x + 2y + 1z = 3
0x − 0y − (38/11)z = −(153/11)
La solución en Matlab es
4.1. SOLUCIÓN DE SISTEMAS LINEALES 87
3.0000
-11.2000
-13.9091
ans =
-0.6579
1.1316
4.0263
Ejemplo 3
Resolver el sistema de ecuaciones
3 −1 −1 x 0
−1 1 0 y = 1
−1 0 1 z 1
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
0
1.0000
1.5000
88 SISTEMAS DE ECUACIONES
ans =
2.0000
3.0000
3.0000
1 2 6
4 8 −1
−2 3 5
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
1
-3
Inf
4.1. SOLUCIÓN DE SISTEMAS LINEALES 89
ans =
NaN
NaN
NaN
1
3
-3
ans =
0.0057
0.1371
0.1200
4.1.5. 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.
a1,1 a1,2 a1,3 x1 b1
a2,1 a2,2 a2,3 x2 = b2
a3,1 a3,2 a3,3 x3 b3
Ax = Ib
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
a1,1 a1,2 a1,3 1 0 0 b1
a2,1 a2,2 a2,3 0 1 0 b2 (4.4)
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 (4.5)
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 (4.6)
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 (4.7)
a1,1 a1,1 a1,1
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 = − 1 0
b2
a1,1
a
0
a a1,2
a3,2 − 3,1
a a1,3
a3,3 − 3,1 x3 − a3,1 0 1 b3
a1,1 a1,1 1,1
a1,2 a1,3 1 1
1 a1,1 a1,1 a1,1 0 0 a1,1 b1
a a1,2 a a1,3 a2,1 a2,1
a2,2 − 2,1 a2,3 − 2,1 − a1,1 b2 − a1,1
0 a1,1 a1,1 1 0 b1
a3,1 a1,2 a3,1 a1,3 a3,1 a3,1
0 a3,2 − a1,1 a3,3 − a1,1 − a1,1 0 1 b3 − a1,1 b1
Segundo Paso
Para un sistema equivalente
1 a′1,2 a′1,3
′ ′
x1 c1,1 0 0 b1
0 a′2,2 a′2,3 x2 = c′2,1 1 0 b′2 (4.8)
0 a′3,2 a′3,3 x3 c′3,1 0 1 b′3
Despejamos x2 de la ecuación 2
Sustituyendo en la ecuación 1
!
c′2,1 b′1 + b′2 − a′2,3 x3
x1 + a′1,2 + a′1,3 = c′1,1 b′1
a′2,2
! !
a′1,2 a′2,3 a′1,2 c′2,1 a′2,1 ′
x1 + 0x2 + a′1,3 − x3 = c′1,1 b′1 − − b
a′2,2 a′2,2 a′2,2 2
Sustituyendo en la ecuación 3
! !
a′3,2 a′2,3 a′3,2 c′2,1 a′3,2 ′
0x1 + 0x2 + a′3,3 − x3 = c′3,1 − b′1 − b 2 + b′ 3
a′2,2 a′2,2 a′2,2
92 SISTEMAS DE ECUACIONES
En forma matricial
Tercer Paso
Dado el sistema equivalente
1 0 a′′1,3
′′
c1,1 c′′1,2 0
′′
x1 b1
0 1 a′′2,3 x2 = c′′2,1 c′′2,2 0 b′′2
0 0 a′′3,3 x3 c′′3,1 c′′3,2 1 b′′3
a′′ ′′
1,3 c3,2 a′′
c′′1,1 c′′1,2 − ′′ − a1,3
′′
a3,3
b′′1
1 0 0 x1 3,3
a′′ a′′
′′
c′′2,1
0 1 0 x2 = ′′
c2,2 − a2,3′′ − 2,3
a′′ b2
3,3 3,3
0 0 1 x3 b′′3
′′
c′′
c3,1 3,2 1
a′′
3,3
′′
a3,3 a′′
3,3
Como resultado tenemos que la matrix c′′ es la inversa de nuestro sistema y b′′ 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
4.1. SOLUCIÓN DE SISTEMAS LINEALES 93
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
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
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
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
4.1. SOLUCIÓN DE SISTEMAS LINEALES 95
a3,1 a3,2 a3,M
a1,1 − a1,3 a1,2 − a1,3 0 ··· a1,M − a1,3
a3,3 a3,3 a3,3
a a3,2 a
a −
2,1 a2,3 a3,1
3,3
a2,2 − a2,3 a3,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.
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
96 SISTEMAS DE ECUACIONES
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
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
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
4.1. SOLUCIÓN DE SISTEMAS LINEALES 97
19123
El tercer renglón se normaliza entonces al dividirlo entre 191
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
El código en Matlab es
function [x, Ainv] = Gauss_Jordan(A, b)
N = length(b);
for n = 1:N;
B(n,:) = B(n, :)./B(n,n);
for m = 1: N
if n ~= m
B(m,:) = B(m,:) - B(n,:)*B(m,n);
end;
end;
end;
x = B(:,2*N+1);
98 SISTEMAS DE ECUACIONES
Como resultado tenemos en primer termino la matriz identidad resultado de las elimina-
ciones, en segundo lugar la matriz inversa y finalmente la solución del sistema de ecuacio-
nes.
>> [b, Ainv] = Gauss_Jordan([3, -0.1, -0.2; 0.1 7 -0.3; 0.3 -0.2 10],
[7.85; -19.3; 71.4])
x =
3.0000
-2.5000
7.0000
Ainv =
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
1 1 2 −1
1 3 −6 7
2 −1 2 0
Primer iteración
4.1. SOLUCIÓN DE SISTEMAS LINEALES 99
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
Ejemplo 3
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 iteración
100 SISTEMAS DE ECUACIONES
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
4.1. SOLUCIÓN DE SISTEMAS LINEALES 101
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
b =
-0.0508
0.2372
0.6707
0.4037
Ainv =
f1 (x1 , x2 , · · · , xn ) = 0
f2 (x1 , x2 , · · · , xn ) = 0
.. .
. = ..
fn (x1 , x2 , · · · , xn ) = 0
f1 (x, y) = x2 + xy − 10 = 0
f2 (x, y) = y + 3xy 2 − 57 = 0
podemos graficar la curva de nivel con la función contour de matlab. El código correspon-
diente es el siguiente y en la Figura 4.11 se muestra las curvas de nivel obtenidas.
xa = [-10:0.01:10];
ya = [-10:0.01:10];
[x,y] = meshgrid(xa, ya);
contour(x,y, f1,[0,0],’k’)
hold on
contour(x,y, f2,[0,0],’r’)
4.2. MÉTODOS PARA SISTEMAS NO LINEALES 103
f1 (x1 , x2 , · · · , xn ) = 0
f2 (x1 , x2 , · · · , xn ) = 0
.. .
. = ..
fn (x1 , x2 , · · · , xn ) = 0
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 )
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;
104 SISTEMAS DE ECUACIONES
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
√
f1 (x1 , x2 ) = x1 − 10 − x1 x2
r
57 − x2
f2 (x1 , x2 ) = x2 −
3x1
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)
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)
4.2. MÉTODOS PARA SISTEMAS NO LINEALES 105
Iteración 2
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
3x 3(1.9405)
1
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
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.
√
f1 (x1 , x2 ) = x1 − 2x1 + x2 − 0.5
r
4 − x21
f2 (x1 , x2 ) = x2 −
4
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)
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
4.2. MÉTODOS PARA SISTEMAS NO LINEALES 107
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.
(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
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:
108 SISTEMAS DE ECUACIONES
En forma compacta
h i−1
x(t+1) = x(t) − J(x(t) ) f (x(t) )
while 1
x2 = x1 - inv(J(x1))* f(x1);
if norm((x1-x2)./x2) < 0.0001 break;
else x1 = x2;
end;
end;
r = x2;
Ejemplo 1
Resolver el siguiente sistema de ecuaciones dado por y cuya gráfica se muestra en la Figura
4.11
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
10
x2
0
-5
-10
-10 -5 0 5 10
x1
" # −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
4.2. MÉTODOS PARA SISTEMAS NO LINEALES 111
Note que solamente 4 iteraciones son necesarias para llegar a la solución x = [2, 3]T
function f = f1(x)
n = length(x);
f = zeros(n,1);
El Jacobiano
function J = J1(x)
n = length(x);
J = zeros(n, n);
y la ejecución
ans =
2.0000
3.0000
Ejemplo 2
Resolver el siguiente sistema de ecuaciones dado por las ecuaciones y cuya solución gráfica
se muestra en la figura 4.12
112 SISTEMAS DE ECUACIONES
x2
0
-1
-2
-2 -1 0 1 2
x1
El Jacobiano es
2x1 − 2 −1
J(x) =
2x1 8x2
y el arreglo de funciones
" # −1
(1)
x1 0 −2 −1 −0.5 −0.25
(1) = − =
x2 1 0 8 0.0 1.00
4.2. MÉTODOS PARA SISTEMAS NO LINEALES 113
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
n = length(x);
f = zeros(n,1);
ans =
-0.2222
0.9938
4.2.4. Ejemplo
Para el circuito que se muestra en la figura 4.13, 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.
V − R1 (i1 + i2 ) − L1 /i1 = 0
V − R1 (i1 + i2 ) − R2 i2 − L2 /i2 = 0
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);
V = 100;
4.2. MÉTODOS PARA SISTEMAS NO LINEALES 115
R1 = 2;
R2 = 3;
L1 = 100;
L2 = 200;
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;
I =
116 SISTEMAS DE ECUACIONES
1.0728
2.3185
Ejemplo
f (x(k) )
x(k+1) = x(k) −
f ′ (x(k+1) )
10
(k+1) (k) x(k) −1
x =x − 9
10 x(k)
l0 = l1 + l2 (5.9)
l1 l2
= (5.10)
l0 l1
l1 l2
= (5.11)
l1 + l2 l1
l2 2
Definimos R = l1√y sustituimos en la ecuación 5.11 para obtener R + R − 1, cuya solución
positiva es R = 5−1
2 . R es definido como la razón dorada y fue un número ampliamente
utilizado por los Griegos en su arquitectura.
El algoritmo de la razón dorada es:
La implementación en Matlab es:
function y = Razon_Dorada(f, xl, xu)
R = (sqrt(5) -1)/2
t = 0;
119
120 OPTIMIZACIÓN
5.1.1. Ejemplo 1
Use la búsqueda de la sección dorada para encontrar el mı́nimo de la función f (x) =
2
−2seno(x) + x10 en el intervalo [0,4].
Los resultados de este ejemplo se muestran en la siguiente tabla
k xl fl x2 f2 x1 f1 xu fu d
0 0.0000 0.0000 1.5279 -1.7647 2.4721 -0.6300 4.0000 -3.1136 2.4721
1 0.0000 0.0000 0.9443 -1.5310 1.5279 -1.7647 2.4721 -0.6300 1.5279
2 0.9443 -1.5310 1.5279 -1.7647 1.8885 -1.5432 2.4721 -0.6300 0.9443
3 0.9443 -1.5310 1.3050 -1.7595 1.5279 -1.7647 1.8885 -1.5432 0.5836
4 1.3050 -1.7595 1.5279 -1.7647 1.6656 -1.7136 1.8885 -1.5432 0.3607
... ... ... ... ... ... ... ... ... ...
17 1.4267 -1.7757 1.4271 -1.7757 1.4274 -1.7757 1.4278 -1.7757 0.0007
ans =
1.4276
5.1.2. Ejemplo 2
2
Use la búsqueda de la sección dorada para encontrar el mı́nimo de la función f (x) = −xe−x
en el intervalo [0,4].
Los resultados de este ejemplo se muestran en la siguiente tabla
k xl fl x2 f2 x1 f1 xu fu
0 0.0000 0.0000 1.5279 -0.1480 2.4721 -0.0055 4.0000 -0.0000
1 0.0000 0.0000 0.9443 -0.3871 1.5279 -0.1480 2.4721 -0.0055
2 0.0000 0.0000 0.5836 -0.4151 0.9443 -0.3871 1.5279 -0.1480
3 0.0000 0.0000 0.3607 -0.3167 0.5836 -0.4151 0.9443 -0.3871
4 0.3607 -0.3167 0.5836 -0.4151 0.7214 -0.4287 0.9443 -0.3871
... ... ... ... ... ... ... ... ...
22 0.7071 -0.4289 0.7071 -0.4289 0.7071 -0.4289 0.7072 -0.4289
>> Razon_Dorada(@fun_2, 0, 4)
ans =
0.7071
1
q(x) = f (x(k) ) + f ′ (x(k) )(x − x(k) ) + f ′′ (x(k) )(x − x(k) )2
2
llamaremos a x(k+1) el mı́nimo de q(x). Para calcularlo, es necesario que la segunda derivada
f ′′ (x(k) ) sea positiva. La función q(x) es entonces estrictamente convexa y tiene un mı́nimo
único en xk+1 dado por
q ′ (x(k+1) ) = 0
Esto da lugar sistema lineal de ecuaciones:
f ′ (x(k) )
x(k+1) = x(k) − (5.12)
f ′′ (x(k) )
La implementación en Matlab es
5.2. OPTIMIZACIÓN NO-RESTRINGIDA. MÉTODO DE NEWTON 123
5.2.2. Ejemplo 1
x2
Calcular el máximo de la función f (x) = −2seno(x) + 10 utilizando el método de Newton
con un valor inicial x0 = 0.5.
x
f ′ (x) = − 2cos(x)
5
1
f ′′ (x) = + 2sin(x)
5
k = 0
(1) −1.6551
x = 0.5000 − = 1.9282
1.1588
k = 1
1.0854
x(2) = 1.9282 − = 1.4047
2.0735
k = 2
−0.0495
x(3) = 1.4047 − = 1.4275
2.1725
k = 3
8.20126 × 10−5
x(4) = 1.4275 − = 1.4275
2.1795
k = 4
2.02094 × 10−10
x(5) = 1.4275 − = 1.4275
2.1795
ans =
1.4276
function y = der_fun_1(x)
y = x/5 - 2*cos(x);
function y = d2f1(x)
y = 1/5 + 2*sin(x);
5.2.3. Ejemplo 2
2
Dada la función f (x) = −xe−x , calcular el mı́nimo utilizando el método de Newton para
un valor inicial x0 = −1.0.
2
f ′′ (x) = (6x − 4x3 )e−x
k = 0
(1) 0.3679
x = −1.0000 − = −0.5000
−0.7358
k = 1
(2) −0.3894
x = −0.5000 − = −0.7000
−1.9470
k = 2
−0.0123
x(3) = −0.7000 − = −0.7071
−1.7325
k = 3
(4) −0.0001
x = −0.7071 − = −0.7071
−1.7156
k = 4
−0.0000
x(5) = −0.7071 − = −0.7071
−1.7155
ans =
-0.7071
function y = der_fun_2(x)
y = (2*x^2 -1)*exp(-x^2);
function y = der2_fun_2(x)
y = (6*x-4*x^3)*exp(-x^2);
1
q(x) = f (x(k) ) + ∇f T (x(k) )(x − x(k) ) + (x − x(k) )T ∇2 f (x(k) )(x − x(k) )
2
126 OPTIMIZACIÓN
∇q (x(k+1) ) = 0
h i−1
x(k+1) = x(k) − ∇2 f (x(k) ) ∇f (x(k) ) (5.13)
while 1
d= inv(H(x1))* g(x1);
x2 = x1 - d;
if norm(d) < 0.0001 break;
else x1 = x2;
end;
end;
r = x2;
5.2.5. Ejemplo 1
2 2 (0) (0)
Minimizar es f (x) = x1 e−x1 −x2 , con un punto inicial x1 = −0.5, x2 = −0.1 aplicando el
Método de Newton.
" 2 2
#
(1 − 2x21 )e(−x1 −x2 )
∇f (x1 , x2 ) = 2 2
−2x1 x2 e(−x1 −x2 )
function g = Gradiente_f1(x)
n = length(x);
g = zeros(n, 1);
val_e = exp(-x(1)^2-x(2)^2);
g(1) = val_e*(1-2*x(1)^2);
g(2) = val_e*(-2*x(1)*x(2));
El Hessiano es:
" 2 2 2 2
#
2 (−6x1 + 4x31 )e(−x1 −x2 ) (−2x2 + 4x21 x2 )e(−x1 −x2 )
∇ f (x1 , x2 ) = 2 2 2 2
(−2x2 + 4x21 x2 )e(−x1 −x2 ) (−2x1 + 4x1 x22 )e(−x1 −x2 )
function H = Hessiano_f1(x)
n = length(x);
H = zeros(n, n);
val_e = exp(-x(1)^2-x(2)^2);
H(1,1) = val_e*(-6*x(1)+4*x(1)^3);
H(1,2) = val_e*(-2*x(2)+4*x(1)^2*x(2));
H(2,1) = H(1,2);
H(2,2) = val_e*(-2*x(1)+4*x(1)*x(2)^2);
Primer iteración
" # −1
(1)
x1 −0.5000 1.9276 0.0771 0.3855 −0.7049
(1) = − =
x2 −0.1000 0.0771 0.7556 −0.0771 0.0230
Segunda iteración
128 OPTIMIZACIÓN
" # −1
(2)
x1 −0.7049 1.7199 −0.0002 0.0038 −0.7071
(2) = − =
x2 0.0230 −0.0002 0.8564 0.0197 −0.0000
Tercer iteración
" # −1
(3)
x1 −0.7071 1.7155 0.0000 0.0180 −4 −0.7071
(3) = − × 10 =
x2 −0.0000 0.0000 0.8578 −0.2114 −0.0000
k f (x(k) )
0 -0.3855
1 -0.4287
2 -0.4289
3 -0.4289
ans =
-0.7071
0.0000
0.5 2
0.4
0.3 1
0.2
0.1 0
0
−1
−0.1
−0.2
−2
−0.3 4
−0.4 2
−3
−0.5 0
−4 −3 −2
−2 −1 0 1 −4
2 3 −4 −4 −3 −2 −1 0 1 2 3 4
4
2 2
Figura 5.14: Mı́nimo para la función f (x1 , x2 ) = x1 e(−x1 −x2 )
5.2.6. Ejemplo 2
(0)
Minimizar es f (x) = −3.5x1 − 2x2 − x22 + x41 + 2x1 x2 + x42 , con un punto inicial x1 = 1.0,
(0)
x2 = 1.0 aplicando el Método de Newton.
El vector gradiente para esta función es:
12x21
2 2
∇ f (x1 , x2 ) =
2 −2 + 12x22
function H = Hessiano_f2(x)
n = length(x);
H = zeros(n, n);
130 OPTIMIZACIÓN
H(1,1) = 12*x(1)^2;
H(1,2) = 2;
H(2,1) = 2;
H(2,2) = -2+12*x(2)^2;
Primer iteración
" # −1
(1)
x1 1 12 2 2.5000 0.8190
(1) = − =
x2 1 2 10 2.0000 0.8362
Segunda iteración
" # −1
(2)
x1 0.8190 8.0485 2.0000 0.3695 0.7820
(2) = − =
x2 0.8362 2.0000 6.3909 0.3044 0.8001
Tercer iteración
" # −1
(3)
x1 0.7820 7.3385 2.0000 0.0132 0.7807
(3) = − =
x2 0.8001 2.0000 5.6828 0.0129 0.7983
cuarta iteración
" # −1
(4)
x1 0.7807 7.3139 2.0000 0.1610 −4 0.7807
(4) = − × 10 =
x2 0.7983 2.0000 5.6483 0.3115 0.7983
ans =
5.2. OPTIMIZACIÓN NO-RESTRINGIDA. MÉTODO DE NEWTON 131
-0.7807
0.7983
En la figura 5.15 se muestras las curvas de nivel correspondientes a este ejemplo las cuales
fueron determinadas utilizando el siguiente código.
1.5
0.5
x2
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1
x(k+1) = −A−1 b
5.2.8. Ejemplo 3
Dada la función f (x) = 10x21 + x22 y un punto inicial x(0) = [1, 2]T , calcular el mı́nimo
utilizando el algoritmo de Newton.
Para esta función el gradiente y el Hessiano es:
20x1
∇f (x) =
2x2
y el Hessiano es:
20 0
A(x) =
0 2
Aplicando la ecuación 5.13 tenemos
−1
(1) 1 20 0 20 0
x = − =
2 0 2 4 0
5.2.9. Ejemplo 4
2 2 (0) (0)
Minimizar es f (x) = x1 e−x1 −x2 , con un punto inicial x1 = −0.5, x2 = −0.5 aplicando el
Método de Newton.
Primer iteración
" # −1
(1)
x1 −0.5000 1.5163 0.3033 0.3033 −1.0000
(1) = − =
x2 −0.5000 0.3033 0.3033 −0.3033 1.0000
Segunda iteración
" # −1
(2)
x1 −1.0000 0.2707 0.2707 −0.1353 −1.2500
(2) = − =
x2 1.0000 0.2707 −0.2707 0.2707 1.7500
5.2. OPTIMIZACIÓN NO-RESTRINGIDA. MÉTODO DE NEWTON 133
Tercer iteración
" # −1
(3)
x1 −1.2500 −0.0031 0.0729 −0.0208 −1.3535
(3) = − =
x2 1.7500 0.0729 −0.1256 0.0429 2.0314
Cuarta iteración
" # −1
(4)
x1 −1.3535 −0.0046 0.0280 −0.0069 −1.4416
(4) = − =
x2 2.0314 0.0280 −0.0507 0.0142 2.2629
0.5
0.4
0.3
0.2
0.1
−0.1
4
−0.2
−0.3 2
−0.4 0
−0.5 −2
−4 −3 −2 −1 0 1 −4
2 3 4
2 2
Figura 5.16: Función f (x) = x1 e−x1 −x2 .
caso la solución es un polinomio de grado 4 en cuyo caso tendrá cuatro soluciones con
x∗ = 3. Si aplicamos el método de Newton tendremos una convergencia como se muestra
a continuación
k f (x(k) )
1 2.2500
2 2.4375
3 2.5781
4 2.6836
5 2.7627
6 2.8220
7 2.8665
8 2.8999
9 2.9249
10 2.9437
.. ..
. .
44 3.0000
45 3.0000
Si uno desea aplicar el método a una función arbitraria encontraremos algunas dificultades,
esencialmente debido al hecho de que no encontraremos convergencia global. Si el punto
inicial x(0) esta muy lejos de x∗ el método no podrá converger.
Para comenzar, la aproximación de f (x) dada por q(x) es solamente valida a la vecindad de
xk , el tamaño de paso puede ser controlado a través de una formula iterativa de tipo:
5.2. OPTIMIZACIÓN NO-RESTRINGIDA. MÉTODO DE NEWTON 135
−1
x(k+1) = x(k) − λ(k) ∇2 f (x(k) ) ∇f (x(k) )
donde λ(k) es un escalar, seleccionado por ejemplo, de tal manera que ||x(k+1) −x(k) || no sea
muy largo. También podemos seleccionarlo tal que λ(k) minimice ϕ(λ(k) ) = f (x(k) +λ(k) d(k) )
en la dirección de:
−1
d(k) = − ∇2 f (x(k) ) ∇f (x(k) )
5.2.11. Ejemplo 5
5
Minimizar f (x) = (x−3)5 utilizando el método de Newton y la modificación planteada
utilizando un valor λ(k) = 2.9. Considere un punto inicial x(0) = 2.
k Newton Newton con λ = 2.9
0 2.00000 2.00000
1 2.33333 2.96667
2 2.55556 2.99889
3 2.70370 2.99996
4 2.80247 3.00000
5 2.86831 3.00000
6 2.91221 3.00000
7 2.94147 3.00000
8 2.96098 3.00000
9 2.97399 3.00000
10 2.98266 3.00000
11 2.98844 3.00000
f ′′ (x(k) )f ′ (x(k) )
x(k+1) = x(k) −
f ′′ (x(k) )f ′′ (x(k) ) − f ′ (x(k) )f ′′′ (x(k) )
Ejemplo
(x−3)5
Para la función f (x) = 5 para un valor inicial x0 = 2 la solución es:
La segunda derivada
f ′′ (x) = 4(x − 3)3
y la tercer derivada
f ′′′ (x) = 12(x − 3)2
(1) × (−4)
x(1) = 2 − =3
(4)2 − (1) ∗ (12)
Note que la solución se da en una sola iteración. Para ejecutar hacer
Maximizar
Z = c1 x1 + c2 x2 + · · · + cn xn
donde cj es el pago por cada unidad de la j-esima actividad que se lleva a cabo y xj es la
magnitud o cantidad de la j-esima actividad.
donde ai,j cantidad de i-ésima fuente que se consume por cada unidad de j-ésima actividad
y bi cantidad de la i-ésima fuente que ésta disponible.
Suponga que una planta procesadora de gasolina recibe cada semana una cantidad fija
de materia prima. Esta última se procesa en dos tipos de gasolina de calidad regular y
premium. Estas clases de gasolina son de alta demanda; es decir, se tiene garantizada
su venta y se obtiene diferentes utilidades para la compañı́a. Si embargo su producción
involucra ambas restricciones, tiempo y almacenaje en sitio. Por ejemplo, sólo una de las
clases se puede producir a la vez, y las instalaciones están abiertas solo 80 horas. Además,
existe un limite de almacenamiento para cada uno de los productos. Todos estos factores
se en listan en la siguiente tabla:
138 OPTIMIZACIÓN
7x1 + 11x2 ≤ 77
Las restricciones restantes se pueden desarrollar en una forma similar, la formulación total
resultante para la PL está dada por
Maximizar Z = 150x1 + 175x2
Sujeta a
7x1 + 11x2 ≤ 77
10x1 + 8x2 ≤ 80
x1 ≤ 9
x2 ≤ 6
x1 ≥ 0
x2 ≥ 0
5.3. METODO SIMPLEX 139
7x1 + 11x2 ≤ 77
Se puede definir una variable h1 como la cantidad de gasolina cruda que no se usa para
una nivel de producción en particular (x1 , x2 ). Si esta cantidad se agrega al lado izquierdo
de la restricción, forma la relación exacta.
7x1 + 11x2 + h1 = 77
Ahora se reconoce lo que la variable de holgura nos indicaba. Si esta es positiva, significa
que tiene algo de holgura para esta restricción. Esto es, se cuenta con algo más de recursos
que no han sido utilizado por completo. Si es negativa nos indica que nos hemos excedido en
la restricción. Finalmente, si es cero denota que se cumplió exactamente con la restricción.
Es decir, se dispuso de todo el recurso puesto que esta es la condición donde las lı́neas de
restricción se interceptan, la variable de holgura proporciona un medio para determinar
puntos extremos.
Una variable de holgura diferente se desarrolla para cada ecuación restringida resultando
en lo que se llama versión completamente aumentada.
Z = 150x1 + 175x2
Sujeta a
7x1 +11x2 +h1 = 77
10x1 +8x2 +h2 = 80
x1 +h3 =9
x2 +h4 = 6
Advierta como se han formulado de igual modo las cuatro ecuaciones, de tal manera que
las incógnitas están alineadas en las columnas. Se hizo ası́ para resaltar que ahora se trata
de un sistema de ecuaciones algebraicas lineales.
140 OPTIMIZACIÓN
1. Utilizando la forma estándar, determinar una solución básica factible inicial igualando
a las n − m variables igual a cero (el origen).
function x = Simplex(C, A, b)
N = length(b);
B = [A, diag(ones(N,1)), b];
B = [B; C, zeros(1, N+1)]
M = length(B(1,:));
while 1
[mm, m] = max(B(N+1,:)) ;
[mn, n] = imin(B(1:N, M)./B(1:N, m));
if mm == 0
break;
end;
for k = 1: N+1
if n ~= k
B(k,:) = B(k,:) - B(n,:)*B(k,m);
end;
end;
end;
x=B;
5.3. METODO SIMPLEX 141
5.3.3. Ejemplo 1
Mostrar la región de factibilidad del problema de la Gasolina planteado y resolver utilizando
el método de Simplex.
Z = 150x1 + 175x2
Sujeta a
7x1 +11x2 +h1 = 77
10x1 +8x2 +h2 = 80
x1 +h3 =9
x2 +h4 = 6
Figura 5.17: Para el ejemplo 1, en verde área de Factibilidad con la función objetivo en
rojo para Z = 100.
5.3. METODO SIMPLEX 143
de holgura tiene los valores que hacen las desigualdades ver como igualdades
7.00 11.00 1 0 0 0 77.00
10.00 8.00 0 1 0 0 80.00
1.00 0 0 0 1 0 9.00
0 1.00 0 0 0 1 6.00
150.00 175.00 0 0 0 0 0
Note que en nuestra formulación estándar la función objetivo aparece con los coeficientes
negativos, por esta razón repetiremos el procedimiento hasta que no tengamos ningún
término negativo en el renglón correspondiente a la función de costo.
Primer iteración
Dado que el costo mayor es 175 comenzamos por introducir la variable x2 por lo tanto
la columna es la 2 y sacar la variable de holgura h4 dado que esta restricción es la mas
cercana del punto básico. Por lo tanto el pivote es 4, 2
7.00 11.00 1 0 0 0 77.00 7
10.00 8.00 0 1 0 0 80.00
10
1.00 0 0 0 1 0 9.00 ∞
0 1.00 0 0 0 1 6.00 6
150.00 175.00 0 0 0 0 0
0 −11.00
7.00 0 1.00 0 11.00
10.00 0 0 1.00 0 −8.00 32.00
1.00 0 0 0 1.00 0 9.00
0 1.00 0 0 0 1.00 6.00
150.00 0 0 0 0 −175.00 −1050.00
Segunda Iteración
La siguiente variable a introducir es x1 dado que el mayor de los costos es 150 y sacar la
variable de holgura h1 dado que es la mas cercana de nuestro punto extremo. Por lo tanto
tenemos el elemento 1, 1 como pivote
144 OPTIMIZACIÓN
0 −11.00
7.00 0 1.00 0 11.00 1.57
10.00 0 0 1.00 0 −8.00 32.00
3.2
1.00 0 0 0 1.00 0 9.00
9
0 1.00 0 0 0 1.00 6.00 ∞
150.00 0 0 0 0 −175.00 −1050.00
0 −1.57
1.00 0 0.14 0 1.57
0 0 −1.43 1.00 0 7.71 16.29
0 0 −0.14 0 1.00 1.57 7.43
0 1.00 0 0 0 1.00 6.00
0 0 −21.43 0 0 60.71 −1285.71
Tercer iteración
En la función de Costo tenemos un término positivo de 60.71 corresponden a la variable de
holgura h4 , por lo tanto introducimos la variable h4 y sacamos la variable h2 . Por lo tanto
el pivote para esta iteración es 2, 6
0 −1.57
1.00 0 0.14 0 1.57 −1.00
0 0 −1.43 1.00 0 7.71 16.29
2.11
0 0 −0.14 0 1.00 1.57 7.43
4.73
0 1.00 0 0 0 1.00 6.00 6
0 0 −21.43 0 0 60.71 −1285.71
0 −0.15
1.00 0.20 0 0 4.89
0 0 −0.19 0.13 0 1.00 2.11
0 0 0.15 −0.20 1.00 0 4.11
0 1.00 0.19 −0.13 0 0 3.89
0 0 −10.19 −7.87 0 0 −1413.89
Dado que no tenemos ningún termino de costo negativo, la solución la podemos obtener
de la formulación de la siguiente manera. Para x1 la solución esta en el renglón 1 dado
que hay un 1 en la posición 1,1 y es x1 = 4.89. Para x2 el 1 esta en la posición 4 por lo
tanto la solución es x2 = 3.89. Para h1 no tenemos un 1 en la columna correspondiente y
de igual manera para h2 , por lo tanto tenemos h1 = 0 y h2 = 0. Para h3 tenemos un 1
en el renglón 3 por lo tanto h3 = 4.11. Finalmente para h4 tenemos un 1 en la posición
5.3. METODO SIMPLEX 145
B =
B =
5.3.4. Ejemplo 2
Maximizar la función
Z = 6x1 + 4x2
Sujeta a
2x1 + 2x2 ≤ 160
x1 + 2x2 ≤ 120
4x1 + 2x2 ≤ 280
2.00 2.00 1.00 0 0 160.00
1.00 2.00 0 1.00 0 120.00
4.00 2.00 0 0 1.00 280.00
6.00 4.00 0 0 0 0
Primer iteración
Tomamos como pivote el renglón 3 y columna 1
2.00 2.00 1.00 0 0 160.00 80.0
1.00 2.00 0 1.00 0 120.00 120.0
4.00 2.00 0 0 1.00 280.00 70
6.00 4.00 0 0 0 0
0 1.00 1.00 0 −0.50 20.00
0 1.50 0 1.00 −0.25 50.00
1.00 0.50 0 0 0.25 70.00
0 1.00 0 0 −1.50 −420.00
Segunda Iteración
5.3. METODO SIMPLEX 147
Figura 5.18: Ejemplo 2, Área de Factibilidad en verde y en rojo la función objetivo con
Z=350.
148 OPTIMIZACIÓN
0 1.00 1.00 0 −0.50 20.00
0 0 −1.50 1.00 0.50 20.00
1.00 0 −0.50 0 0.50 60.00
0 0 −1.00 0 −1.00 −440.00
Dado que no tenemos términos negativos en el renglón de costos terminamos con solución
x1 = 60, x2 = 20, h1 = 0, h2 = 20, h3 = 0 y costo Z = 440
B =
2 2 1 0 0 160
1 2 0 1 0 120
4 2 0 0 1 280
6 4 0 0 0 0
ans =
>>
5.3.5. Ejemplo 3
Maximizar la función
5.3. METODO SIMPLEX 149
Z = 3x1 + 4x2
Sujeta a
2x1 + 3x2 ≤ 1200
2x1 + x2 ≤ 1000
4x2 ≤ 800
2.00 3.00 1.00 0 0 1200.00
2.00 1.00 0 1.00 0 1000.00
0 4.00 0 0 1.00 800.00
3.00 4.00 0 0 0 0
Primer Iteración
Tomamos como pivote el renglón 3 y columna 2
2.00 3.00 1.00 0 0 1200.00 400.00
2.00 1.00 0 1.00 0 1000.00 1000.00
0 4.00 0 0 1.00 800.00 200.00
3.00 4.00 0 0 0 0
2.00 0 1.00 0 −0.75 600.00
2.00
0 0 1.00 −0.25 800.00
0 1.00 0 0 0.25 200.00
3.00 0 0 0 −1.00 −800.00
Segunda Iteración
Tomamos como pivote el renglón 1 y columna 1
2.00 0 1.00 0 −0.75 600.00 300.00
2.00
0 0 1.00 −0.25 800.00 400.00
0 1.00 0 0 0.25 200.00 ∞
3.00 0 0 0 −1.00 −800.00
1.00 0 0.50 0 −0.38 300.00
0 0 −1.00 1.00 0.50 200.00
0 1.00 0 0 0.25 200.00
0 0 −1.50 0 0.12 −1700.00
Tercer iteración
Tomamos como pivote el renglón 2 la columna 5
1.00 0 0.50 0 −0.38 300.00 −789.47
0 0 −1.00 1.00 0.50 200.00
400.00
0 1.00 0 0 0.25 200.00 800.00
0 0 −1.50 0 0.12 −1700.00
1.00 0 −0.25 0.75 0 450.00
0 0 −2.00 2.00 1.00 400.00
0 1.00 0.50 −0.50 0 100.00
0 0 −1.25 −0.25 0 −1750.00
B =
B =
hold off;
Figura 5.19: Ejemplo 3, área de Factibilidad en verde y en rojo función objetivo con Z =
1400.
Ajuste de curvas
y = a 0 + 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
Una estrategia, para ajustar a la mejor lı́nea, es minimizar la suma al cuadrado de los
errores para todos los datos disponibles
153
154 AJUSTE DE CURVAS
n
X n
X
E(a) = e2i = (yi − a0 − a1 xi )2
i=1 i=1
e1
e2
E(a) = [e1 , e2 , · · · , en ]
..
.
en
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
∂E(a) ∂
= [y − M a]T [y − M a] = 0
∂a ∂a
∂E(a)
= −2M T [y − M a] = 0
∂a
[M T M ]a = M T y
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
N = length(x);
M = ones(N, 1);
for k=1:grado
M = [M,x.^k];
end
a=inv(M’*M)*M’*y;
end
6.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
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
7.00 28.00 a0 24.00
=
28.00 140.00 a1 119.50
y
6
x
1 2 3 4 5 6 7
Ma = y
[M T M ]a = M T y (6.14)
a = [M T M ]−1 [M T y]
158 AJUSTE DE CURVAS
6.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
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
60
50
40
30
20
x
1 2 3 4 5
f (x1 ) − f (x0 )
f (x) = f (x0 ) + (x − x0 )
x1 − x 0
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
end
6.2.1. Ejemplo 1
Dados los puntos
160 AJUSTE DE CURVAS
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;
q(x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 )
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
por lo tanto
b0 = f (x0 )
despejando tenemos
f (x1 ) − f (x0 )
b1 =
x1 − x 0
f (x1 ) − f (x0 )
f (x2 ) = f (x0 ) + (x2 − x0 ) + b2 (x2 − x0 )(x2 − x1 )
x1 − x 0
despejando a b2
162 AJUSTE DE CURVAS
x2 −x0
f (x2 ) − f (x1 ) + (f (x1 ) − f (x0 )) 1 − x1 −x0
b2 =
(x2 − x0 )(x2 − x1 )
x2 −x1
f (x2 ) − f (x1 ) − (f (x1 ) − f (x0 )) x1 −x0
b2 =
(x2 − x0 )(x2 − x1 )
finalmente obtenemos
4 x4 f (x4 ) f (x4 )
6.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
El coeficiente b2 tenemos
1.791759−1.386294
6−4 − 0.4620981
b2 = = −0.0518731
6−1
En Forma Tabular
q(x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 )
En Matlab dar
ans =
0.5658
6.3.2. Ejemplo 2
Dados los puntos
6.3. INTERPOLACIÓN CUADRÁTICA 165
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 Forma Tabular los coeficientes para interpolar son:
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);
166 AJUSTE DE CURVAS
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;
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
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 − x j
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 − x k
Estas diferencias pueden usarse para evaluar los coeficientes en la ecuación 6.15, las cuales
entonces se sustituirán en la siguiente ecuación, para ası́ obtener el polinomio de interpo-
lación
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;
6.4. FORMULAS DE INTERPOLACIÓN DE NEWTON 169
end;
end
6.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
b0 b1 b2 b3
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
6.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;
170 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;
n
X
fn (x) = Li (x)f (xi )
i=0
donde
6.5. INTERPOLACIÓN DE POLINOMIOS DE LAGRANGE 171
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
n
Y x − xj
Li (x) =
xi − x j
j=0
j ̸= i
x − x1 x − x0
f1 (x) = f (x0 ) + f (x1 )
x0 − x 1 x1 − x 0
suma = 0;
172 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
6.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
x−4 x−1
f1 (x) = 0+ 1.386294 = −0.462098 + 0.462098x
1−4 4−1
2−4 2−1
f1 (2) = 0+ 1.386294 = 0.4620981
1−4 4−1
Como era de esperar, ambos resultados concuerdan con los que se obtuvieron utilizando
interpolación utilizando polinomios de Newton.
6.5. INTERPOLACIÓN DE POLINOMIOS DE LAGRANGE 173
6.5.2. Ejemplo 2
Dados los vectores x = [1, 2] y y = [2, 4] calcular la interpolación de grado n = 1.
x−2
L0 (x) = =2−x
1−2
x−1
L1 (x) = =x−1
2−1
6.5.3. Ejemplo 3
Dados los puntos x = [1, 2, 3] y sus correspondientes valores y = [3, 5, 10] determinar:
a) El polinómio de Lagrange de grado n = 2
(x − 2)(x − 3) x2 − 5x + 6
L0 (x) =
(1 − 2)(1 − 3) 2
(x − 1)(x − 3)
L1 (x) = = −x2 + 4x − 3
(2 − 1)(2 − 3)
(x − 1)(x − 2) x2 − 3x + 2
L2 (x) = =
(3 − 1)(3 − 2) 2
174 AJUSTE DE CURVAS
Sustituyendo tenemos
c) Calcular el polinomio de grado 2 utilizando mı́nimos cuadrados para los puntos da-
dos.
1 1 1
M = 1 2 4
1 3 9
La solución es
a = (M T M )−1 M T y
−1
3 6 14 18 4
a = 6 14 36 43 = −2.5
14 36 98 113 1.5
El polinomio es :
6.5.4. Ejemplo 4
Para los valores de x = [0.5, 1.3, 2.1, 2.9] y y = [1.25, 4.3, 3.9, 2.1] determinar:
a) La aproximación polinomial de grado n = 3
Sustituyendo tenemos
1.0000 0.5000 0.2500 0.1250
1.0000 1.3000 1.6900 2.1970
M =
1.0000
2.1000 4.4100 9.2610
1.0000 2.9000 8.4100 24.3890
176 AJUSTE DE CURVAS
−1
4.0000 6.8000 14.7600 35.9720 11.5500 −3.3191
6.8000 14.7600 35.9720 93.0948 20.4950 11.6203
a= =
14.7600 35.9720 93.0948 249.6967 42.4395 −5.2979
35.9720 93.0948 249.6967 685.4319 96.9382 0.6673
4.5
3.5
2.5
1.5
1
0.5 1 1.5 2 2.5 3
Figura 6.26: Aproximación de tercer grado para cuatro puntos utilizando interpolación de
Lagrange
6.5.5. Ejemplo 5
Dados los valores en la siguiente tabla calcular:
x f (x)
1 0.0000
2 0.6931
3 1.0986
4 1.3863
6.5. INTERPOLACIÓN DE POLINOMIOS DE LAGRANGE 177
f3 (1.4) = 0.322619
178 AJUSTE DE CURVAS
Diferenciación e Integración
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
179
180 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 (x0 ). Recordemos que el área de un rectángulo esta dado como
A = (x1 − x0 )f (x0 )
7.3.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 x0 = 0 x1 = 0.8 .
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 x1 N
X −1
f (x0 )dx ≈ f (x0 + kh) ∗ h
x0 k=0
h = abs(x1-x0)/N;
suma = 0;
for k=0:N-1
7.4. INTEGRACIÓN UTILIZANDO LA REGLA TRAPEZOIDAL 181
I = suma*h;
end
(x − x0 )
f (x) ≈ f1 (x) = f (x0 ) + (f (x1 ) − f (x0 ))
x1 − x 0
f (x1 ) + f (x0 )
I≈ (x1 − x0 )
2
7.4.1. Ejemplo
Dados los vectores x = [1, 2] y y = [2, 4] :
a) calcular la función de interpolación de grado n = 1.
x−2
L0 (x) = =2−x
1−2
x−1
L1 (x) = =x−1
2−1
Z 2 Z 2
2
I= f1 (x)dx = 2xdx = x2 1
= 22 − 12 = 3
1 1
7.4.2. 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 x0 = 0 x1 = 0.8 .
En Matlab la solución es
>> 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
7.5. INTEGRACIÓN POR EL MÉTODO DE REGLA SIMPSON 1/3 183
PN −1
f (x0 ) + 2 k=1 f (x0 + kh) + f (x1 )
I ≈ (x1 − x0 )
2N
h = abs(x1-x0)/N;
suma = f(x0);
for k=1:N-1
suma = suma + 2*f(x0+h*k);
end;
suma = suma + f(x1);
I = (x1-x0)*suma/(2*N);
end
Z x2 Z x2
f (x0 ) + 4f (x1 ) + f (x2 )
I= f (x) ≈ f2 (x)dx = (x2 − x0 )
x0 x0 6
184 DIFERENCIACIÓN E INTEGRACIÓN
7.5.1. Ejemplo
Dados los puntos x = [1, 2, 3] y sus correspondientes valores y = [3, 5, 10] determinar:
(x − 2)(x − 3) x2 − 5x + 6
L0 (x) =
(1 − 2)(1 − 3) 2
(x − 1)(x − 3)
L1 (x) = = −x2 + 4x − 3
(2 − 1)(2 − 3)
(x − 1)(x − 2) x2 − 3x + 2
L2 (x) = =
(3 − 1)(3 − 2) 2
Sustituyendo tenemos
Z 3 Z 3
3
I= f2 (x)dx = (1.5x2 − 2.5x + 4)dx = (0.5x3 − 1.25x + 4x) 1
1 1
7.5.2. 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 x0 = 0 x2 = 0.8 .
Para este caso x1 = (0.8 − 0)/2 = 0.4 y f (x1 ) = f (0.4) = 2.4560.
ans =
1.3675
La implementación para N divisiones de tamaño h en Matlab es:
function I = Simpson13(f, ini, fin, N)
h = abs(fin - ini)/N;
suma =0;
for k=0:N-1
x0 = ini + h*k;
x1 = x0 + h/2;
x2 = x0 + h;
suma = suma + f(x0) + 4*f(x1) + f(x2);
end;
I = h*suma/6;
end
(x − x1 )(x − x2 )(x − x3 )
f (x) ≈ f3 (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 )
Z x3 Z x3
f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )
I= f (x)dx ≈ f3 (x)dx = (x3 − x0 )
x0 x0 8
7.6.1. Ejemplo
Para los valores de x = [0.5, 1.3, 2.1, 2.9] y y = [1.25, 4.3, 3.9, 2.1] determinar:
a) La aproximación polinomial de grado n = 3
Sustituyendo tenemos
Z 2.9 Z 2.9
I= f3 (x)dx = (0.667321x3 − 5.297847x2 + 11.620281x − 3.319081)dx
0.5 0.5
2.9
I = (0.166683x4 − 1.765949x3 + 5.810140x2 − 3.319081x) 0.5
7.6.2. 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 x0 = 0 x3 = 0.8 .
ans =
1.5192
188 DIFERENCIACIÓN E INTEGRACIÓN
suma =0;
for k=0:N-1
x0 = ini + 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
7.7. Ejemplos
7.7.1. Ejemplo 1
Dada la función
Z 0.8 0.8
25 2 200 3 675 4 900 5 400 6
I(x) = f (x)dx = 0.2x + x − x + x − x + x = 1.6405
0 2 3 4 5 6 0
Método de Barras
Dado x0 = 0, x1 = 0.8 y f (x0 ) = 0.2 tenemos
La integral de f0 (x) es
7.7. EJEMPLOS 189
Z 0.8
I= 0.2 dx = 0.2x|0.8
0 = 0.2 × (0.8 − 0) = 0.16
0
Regla Trapezoidal
Dado x0 = 0, x1 = 0.8, f (x0 ) = 0.2 y f (x1 ) = 0.232 la aproximación del polinomio de
Lagrange es :
x − x1 x − x0
f1 (x) = f (x0 ) + f (x1 )
x0 − x 1 x1 − x 0
x − 0.8 x−0
f1 (x) = 0.2 + 0.232
0 − 0.8 0.8 − 0
cuya integral es
Z 0.8
0.8
I= (0.2 + 0.04x) dx = (0.2x + 0.02x2 ) 0
= 0.1728
0
Z 0.8
0.8
I= (0.2 + 11.24x − 14x2 ) dx = (0.2x + 5.62x2 − 4.666667x3 ) 0
= 1.367466
0
(x−0.266666)(x−0.533333)(x−0.8) (x−0)(x−0.533333)(x−0.8)
f3 (x) = (0−0.266666)(0−0.533333)(0−0.8) 0.2 + (0.266666−0)(0.266666−0.533333)(x1−0.8) 1.43272+
(x−0)(x−0.266666)(x−0.8) (x−0)(x−0.266666)(x−0.533333)
(0.533333−0)(0.533333−0.266666)(0.533333−0.8) 3.48718 + (0.8−0)(0.8−0.266666)(0.8−0.533333) 0.232
cuya integral es
Z 0.8
I= (0.2 − 4.582222x + 48.888888x2 − 53.888889x3 ) dx =
0
0.8
I = (0.2x − 2.291111x2 + 16.296296x3 − 13.472222x4 ) 0
= 1.51917
En la siguiente tabla se muestra un resumen de la la solución para cada uno de los métodos
ası́ como el error numérico respecto al valor real de la integral.
Método Solución error
Barras 0160000 |1.6405 − 0.160000| = 1.4805
R. Trapezoidal 0.172800 |1.6405 − 0.172800| = 1.4677
R. Simpson 1/3 1.367466 |1.6405 − 1.367466| = 0.2731
R. Simpson 3/8 1.519170 |1.6405 − 1.519170| = 0.1214
Note que el mejor desempeño lo tiene la Regla Simpson 3/8 en virtud de utilizar polinomios
de grado 3 para llevar a cabo la integración.
Figura 7.27: Ejemplo de las funciones aproximadas para cada uno de los métodos
192 DIFERENCIACIÓN E INTEGRACIÓN
7.7.2. Ejemplo 2
Para la función del Ejemplo 1 calcular los errores de integración cuando el intervalo de
integración 0 ≤ x ≤ 0.8 se divide en N = 5 intervalos del mismo tamaño h = 0.16. Mostrar
la gráficas de las funciones aproximadas en cada intervalo. Concluir
En la siguiente tabla se muestran los
Método Solución error
Barras 1.5373 |1.6405 − 1.5373| = 0.10318
R. Trapezoidal 1.5399 |1.6405 − 1.5399| = 0.1006
R. Simpson 1/3 1.0852 |1.6405 − 1.6401| = 0.0004
R. Simpson 3/8 1.6405 |1.6405 − 1.6403| = 0.0002
En la Figura 7.28 se muestra de izquierda a derechas las gráficas correspondientes a cada
aproximación considerando que el intervalo se divide en 5.
Figura 7.28: Ejemplo de las funciones aproximadas para cada uno de los métodos dividiendo
el intervalo de integración en 5
7.7.3. Ejemplo 3
Para la función del Ejemplo 1 calcular los errores de integración cuando el intervalo de
integración 0 ≤ x ≤ 0.8 se divide en N = 100 intervalos del mismo tamaño h = 0.008.
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
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.
ϕ = 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
193
194 ECUACIONES DIFERENCIALES ORDINARIAS
x = zeros(N+1, 1);
y = zeros(N+1, D);
x(1) = x0;
y(1,:) = y0;
for k = 1:N
y(k+1,:) = y(k,:) + f(x(k), y(k,:))*h;
x(k+1) = x(k) + h;
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.
y ′(k+1) + y ′(k)
ȳ ′ =
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
x = zeros(N+1, 1);
y = zeros(N+1, D);
x(1) = x0;
y(1,:) = y0;
for k = 1:N
%Predictor
k1 = f(x(k), y(k,:) );
k2 = f(x(k)+h, y(k,:) + k1*h);
%corrector
%disp([k1, k2, y + k1]);
x(k+1) = x(k) + h;
y(k+1,:) = y(k,:) + (k1 + k2 )*h/2.0;
end;
end
Después este valor es utilizado para calcular la pendiente en el punto medio y hacer la
solución en todo el intervalo como
La implementación en Matlab es
function [x,y] = Punto_Medio(x0, y0, h, N, f)
D = length(y0);
x = zeros(N+1, 1);
y = zeros(N+1, D);
x(1) = x0;
y(1,:) = y0;
for k = 1:N
xm = x(k) + h/2;
ym = y(k,:) + f(x(k), y(k,:))*h/2;
x(k+1) = x(k) + h;
y(k+1,:) = y(k,:) + f(xm, ym)*h;
end;
end
k1 = f (x(k) , y (k) )
3 3
k2 = f (x(k) + h, y (k) + k1 h)
4 4
8.5. RUNGE-KUTTA 3ER ORDEN 197
(k+1) (k) 1 2
y =y + k1 + k2 h
3 3
function [x,y] = RK(x0, y0, h, N, f)
D = length(y0);
x = zeros(N+1, 1);
y = zeros(N+1, D);
x(1) = x0;
y(1,:) = y0;
for k = 1:N
k1 = f(x(k), y(k,:));
k2 = f(x(k) + 3*h/4, y(k,:)+3*h*k1/4);
x(k+1) = x(k) + h;
y(k+1,:) = y(k,:) + (k1/3 + 2*k2/3)*h;
end;
end
k1 = f (x(k) , y (k) )
1 1
k2 = f (x(k) + h, y (k) + k1 h)
2 2
1
y (k+1) = y (k) +
(k1 + 4k2 + k3 ) h
6
function [x, y] = RK3(x0, y0, h, N, f)
D = length(y0);
x = zeros(N+1, 1);
y = zeros(N+1, D);
x(1) = x0;
y(1,:) = y0;
for k = 1:N
k1 = f(x(k), y(k,:));
k2 = f(x(k) + h/2, y(k,:) + k1*h/2);
k3 = f(x(k) + h, y(k,:) - k1*h + 2*k2*h);
x(k+1) = x(k) + h;
y(k+1,:) = y(k,:) + (k1 + 4*k2 + k3)*h/6;
end;
end
k1 = f (x(k) , y (k) )
1 1
k2 = f (x(k) + h, y (k) + k1 h)
2 2
1 1
k3 = f (x(k) + h, y (k) + k2 h)
2 2
k4 = f (x(k) + h, y (k) + k3 h)
8.7. EJEMPLO 199
1
y (k+1) = y (k) + (k1 + 2k2 + 2k3 + k4 ) h
6
x = zeros(N+1, 1);
y = zeros(N+1, D);
x(1) = x0;
y(1,:) = y0;
for k = 1:N
k1 = f(x(k), y(k,:));
k2 = f(x(k) + h/2, y(k,:) + k1*h/2);
k3 = f(x(k) + h/2, y(k,:) + k2*h/2);
k4 = f(x(k) + h, y(k,:) + k3*h);
x(k+1) = x(k) + h;
y(k+1,:) = y(k,:) + (k1 + 2*k2 + 2*k3 + k4)*h/6;
end;
end
8.7. Ejemplo
8.7.1. Una ecuación diferencial sencilla
Consideremos la ecuación diferencial
Z Z
dy
dx = f (x, y)dx = −0.5x4 + 4x3 − 10x2 + 8.5x
dx
200 ECUACIONES DIFERENCIALES ORDINARIAS
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.
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]’;
[x, abs(solucion_real(x)-y1),...
abs(solucion_real(x)-y2),...
abs(solucion_real(x)-y3),...
abs(solucion_real(x)-y4),...
abs(solucion_real(x)-y5),...
abs(solucion_real(x)-y6)]
function dy = ecuacion(x, y)
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
202 ECUACIONES DIFERENCIALES ORDINARIAS
8.7.2. Circuito RL
Para el circuito RL que se muestra en la figura 8.29 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
Solución
La ecuación diferencial para este problema es
di V − Ri
=
dt L
R = 2.0; % Ohms
L = 0.3; % Henries
V = 10.0; % Volts
di = (V - R *i)/L;
end
La solución real de la ecuación diferencial, en función del tiempo, es:
V
i(t) = 1 − e(−Rt/L)
R
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 8.30 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.)
% Numero de puntos
N = 100;
% Valores Iniciales
tini = 0;
tfin = 1;
204 ECUACIONES DIFERENCIALES ORDINARIAS
i_ini = 0;
h = (tfin-tini)/N;
xlabel(’tiempo (seg.)’);
ylabel(’Corriente (Amp.)’);
title(’Solucion Circuito RL’);
R = 2.0; % Ohms
L = 0.3; % Henries
V = 10.0; % Volts
t = [t0:h:N*h]’;
Tau = L/R;
ic = 5 *(1 - exp(-t/Tau));
end
function di = Circuito_RL(t, i)
%Parametros
8.7. EJEMPLO 205
R = 2.0; % Ohms
L = 0.3; % Henries
V = 10.0; % Volts
di = (V - R *i)/L;
end
end
d2 x
dt2
=0
d2 y
dt2
=a
Las ecuaciones diferenciales en el caso del tiro parabólico pueden representarse por un
sistema de ecuaciones con cuatro variable el desplazamiento horizontal x, el desplazamiento
vertical y, la velocidad horizontal vx y la velocidad vertical vy tal como se muestra a
continuación:
dx
dt = vx
d2 x
dt2
=0≡
dvx
=0
dt
dy
dt = vy
d2 x
dt2
= −g ≡
dvy
= −g
dt
dv(1) = v(2);
dv(2) = 0;
dv(3) = v(4);
206 ECUACIONES DIFERENCIALES ORDINARIAS
dv(4) = -g;
end
Para el caso de el tiro parabólico encontrar la solución del sistema de ecuaciones diferen-
ciales, para un tiempo tf = 2.0 seg, asumiendo los siguientes inicialmente x(0) = 0 mts,
y0 = 5 mts, vx (0) = 10 m/s, vy (0) = 10 m/s y un incremento h = 0.001.
La solución en matlab es:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Datos
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = floor((tf - t0)/h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Graficas de la solución
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(d(:,1), d(:,3));
grid on
title(’Tiro Parabólico’);
xlabel(’x (mts.)’);
ylabel(’y (mts.)’);
La solución es
8.7. EJEMPLO 207
>> tiro_parabolico
x(2) = 0 + 10 ∗ 2 = 20
En la figura 8.31 se muestra la gráfica del desplazamiento x(t) vs y(t) para la solución del
tiro parabólico.
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:
208 ECUACIONES DIFERENCIALES ORDINARIAS
Tiro Parabólico
11
10
9
y (mts.)
5
0 5 10 15 20 25
x (mts.)
dx
=v
dt
dv F − Cv − Kx
=
dt M
% Numero de puntos
N = 1000;
% Valores Iniciales
tini = 0;
tfin = 40;
d_ini = 2;
v_ini = 0;
h = (tfin-tini)/N;
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.)
subplot(2,2,1)
plot(t, x(:,1));
xlabel(’tiempo (seg.)’);
ylabel(’Desplazamiento (m.)’);
subplot(2,2,2)
plot(t, x(:,2));
xlabel(’tiempo (seg.)’);
ylabel(’Velocidad (m/s.)’);
subplot(2,2,3)
plot(x(:,1), x(:,2));
xlabel(’Desplazamient (m.)’);
ylabel(’Velocidad (m/s.)’);
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);
Considere el circuito serie alimentado por una fuente VS , con una resistencia R y un capa-
citor C tal como se muestra en la Figura ??. Para este circuito la ecuación que lo rige esta
dado por:
212 ECUACIONES DIFERENCIALES ORDINARIAS
Z
1
Vs − Ri − idt = 0
C
dq 1
Vs − R − q=0
dt C
dq q(h) − q(h − 1)
i= ≈
dt h
8.7. EJEMPLO 213
ic = (q(2:N)-q(1:N-1))/h;
tc = [t0:h:tf];
tc = tc(1:N-1);
function Circuito_RC()
% Numero de puntos
N = 10000;
% Valores Iniciales
t0 = 0;
tf = 0.13;
q0 = 0;
h = (tf-t0)/N;
subplot(2,1,1);
plot(t, q(:,1));
xlabel(’tiempo (seg.)’);
ylabel(’Carga (Coulombs.)’);
title(’Solucion Circuito RC’);
% Calculo de la derivada
ic = (q(2:N)-q(1:N-1))/h;
tc = [t0:h:tf];
tc = tc(1:N-1);
subplot(2,1,2);
plot(tc, ic);
xlabel(’tiempo (seg.)’);
ylabel(’Corriente (Amp.)’);
title(’Solucion Circuito RC’);
214 ECUACIONES DIFERENCIALES ORDINARIAS
function dq = CtoRC(t, q)
%Parametros
R = 2.0; % Ohms
C = 0.01; % Farads
V = 10.0; % Volts
dq = (V*C - q)/(R*C);
end
end
Figura 8.35: Solución numérica del circuito RC. En la parte de arriba la carga con respecto
al tiemo y abajo la corriente en función del tiempo.
Ecuaciones diferenciales parciales
dy
C1 + C2 y = C 3
dx
yk − yk−1
C1 + C2 yk = C3
h
215
216 ECUACIONES DIFERENCIALES PARCIALES
(C1 + C2 h)y1 − C1 y0 = C3 h
(C1 + C2 h)y2 − C1 y1 = C3 h
(C1 + C2 h)y3 − C1 y2 = C3 h
(C1 + C2 h)y4 − C1 y3 = C3 h
..
.
(C1 + C2 h)yN − C1 yN −1 = C3 h
(C1 + C2 h) 0 0 ··· 0 y1 C3 h + C1 y0
−C1 (C1 + C2 h) 0 ··· 0
y2
C3 h
0 −C1 (C1 + C2 h) · · · 0
y3 =
C3 h
.. .. .. .. .. .. ..
. . . . . . .
0 0 0 −C1 (C1 + C2 h) yN C3 h
A = zeros(N-1, N-1);
b = zeros(N-1, 1);
for k=1:N-1
A(k,k) = C(1) + C(2)*h;
b(k) = C(3)*h;
if k>1
A(k, k-1) = -C(1);
else
b(k) = b(k)+ C(1)*y0;
end;
end;
y = [y0; A\b];
end
9.1. FÓRMULA DE INTEGRACIÓN POR EL MÉTODO EXPLÍCITO DE DIFERENCIAS DIVIDIDAS FINITA
9.1.1. Ejemplo
Consideremos el caso del circuito RL planteado en el capı́tulo anterior, el cual tenı́a una
resistencia R = 2 ohms, una inductancia L = 0.3 H y estaba alimentado por una fuente de
voltaje de V = 10 volts dado por la ecuación diferencial
di
L + Ri = V
dt
0.3020 0 0 0 i1 0.0100
−0.3000 0.3020 0 0 i2
0.0100
=
0 −0.3000 0.3020 0 i3 0.0100
0 0 −0.3000 0.3020 i4 0.0100
0
i1
i2 0.0331
=
i3 0.0660
0.0987
i4
0.1311
4.5
3.5
3
i(t)
2.5
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
11.1. Tarea 1
Escribir código en Matlab que permita imprimir, utilizando un solo asterisco, condicionales
y ciclos: a) Un Triángulo
*
**
***
****
*****
******
*******
b Un cuadro
*******
* *
* *
* *
* *
* *
*******
11.2. Tarea 2
Dada la función f (x) = sen10x + cos3x hacer:
2.- Utilizando el método gráfico determinar los mı́nimos en el intervalo [10, 15].
219
220 TAREAS
11.3. Tarea 3
Dada la función f (x) = sen(10x) + cos(3x) la cual tiene un cruce por cero en el intervalo
x∗ ∈ [−3.6, −3.2], calcular el cruce por cero en esta región utilizando los métodos de Bisec-
ción, Regla Falsa, Punto Fijo y Newton Raphson. Reportar una tabla para N iteraciones
como la siguiente para cada uno de los métodos.
Biseccion R. Falsi P Fijo N Raphson
k xBk f (xBk ) xRk f (xRk ) xP k f (xP k ) xN k f (xN k )
0 xB0 f (xB0 ) xR0 f (xR0 ) xP 0 f (xP 0 ) xN 0 f (xN 0 )
1 xB1 f (xB1 ) xR1 f (xR1 ) xP 1 f (xP 1 ) xN 1 f (xN 1 )
2 xB2 f (xB2 ) xR2 f (xR2 ) xP 2 f (xP 2 ) xN 2 f (xN 2 )
.. .. .. .. .. .. .. .. ..
. . . . . . . . .
N xBN f (xBN ) xRN f (xRN ) xP N f (xP N ) xN N f (xN N )
donde N es el número máximo de iteraciones en el cual todos los métodos converge. Hacer
una gráfica con los valores de la Tabla de k vs f (xB k), f (xR k), f (xP k) y f (xN k) con la
solución de todos los métodos.
11.4. Tarea 4
Escribir el código en Matlab para realizar la división sintética de un polinomio fn (x) =
an xn + an−1 xn−1 · · · + a1 x + a0 con un monomio x + s.
11.5. Tarea 5
a) Calcular las raı́ces del polinomio x4 − 1 utilizando el método de Bairstow.
2
b) Dada la función 0.1 + (x − 3)e−(x−3) determinar, la expanción en Serie de Taylor con
N = 5 y a = 0. Para el polinomio resultante calcular las raı́ces utilizando el Método de
Bairstow.
11.6. Tarea 6
1.- Para el sistema de ecuaciones
10 x + 2 y - z = 5
2 x + 15 y = 2
-x + 20 z = 3
11.7. TAREA 7 221
Resolver a mano todas las operaciones para calcular el sistema Triangular Superior con el
método de Eliminación Gaussiana
11.7. Tarea 7
a) Calcular el mı́nimo de la función f (x) = (x − 2)4 utilizando el algoritmo de Razón
Dorada en el intervalo [1, 5] b) Calcular el mı́nimo de la función f (x1 , x2 ) = (x1 − 2)4 +
(x1 − 2)2 x22 + (x2 + 1)2 , utilizando el método de Newton.
En ambos casos hacer al menos tres iteraciones a mano y mostrar el resultado que da la
implementación en Matlab.
11.8. Tarea 8
Buscar una aplicación en Ingenierı́a de los temas revisados en los capı́tulos 1 al 5.
11.9. Tarea 9
Dado el sistema de ecuaciones
11.10. Tarea 10
Dadas la funciones a optimizar Z1 y Z2 , dadas por la siguiente expresiones,
222 TAREAS
Z1 = 6x1 + 4x2
Sujeta a
2x1 + 2x2 ≤ 160
x1 + 2x2 ≤ 120
4x1 + 2x2 ≤ 280
Z2 = 3x1 + 4x2
Sujeta a
2x1 + 3x2 ≤ 1200
2x1 + x2 ≤ 1000
4x2 ≤ 880
calcular los puntos del polı́gono de la región de factibilidad (para cada problema) y evaluar
la función de costo correspondiente en cada uno de ellos.
11.11. Tarea 11
a) Para la implementación en Matlab del algoritmo Simplex hacer una modificación que
permita hacer que esta función regrese:
[z, x, h] = Simplex[c, A, b]
donde z es el valor de la función de costo, x el valor de la variables en la solución y h el
valor de las variables de holgura.
b) Optimizar utilizando el método simplex, el siguiente problema. Todos los cálculos de-
berán ser realizados a mano y comprobar que la solución es la correcta utilizando el algo-
ritmo implementado en Matlab.
11.12. Tarea 12
Dados los puntos
x ln(x)
1 0.0000
2 0.6931
3 1.0986
4 1.3863
5 1.6094
6 1.7918
a) Encontrar los polinomios de mı́nimos cuadrados fk (x) de grado k = 1, 2, 3, 4, 5,
b) Realizar la interpolación utilizando interpolación de Newton de grado k = 1, 2, 3, 4, 5
y
c) Llenar la siguiente tabla, donde se hace la comparación del valor real de la función y los
valores de los polinomios de mı́nimos cuadrados y la interpolación de Newton
Grado Polinomio Polinomio Interpolación Interpolación
k fk (1.5) Error fˆk (1.5) Error
1 f1 (1.5) (0.4055 − f1 (1.5))2 fˆ1 (1.5) (0.4055 − fˆ1 (1.5))2
2 f2 (1.5) (0.4055 − f2 (1.5))2 fˆ2 (1.5) (0.4055 − fˆ2 (1.5))2
3 f3 (1.5) (0.4055 − f3 (1.5))2 fˆ3 (1.5) (0.4055 − fˆ3 (1.5))2
4 f4 (1.5) (0.4055 − f4 (1.5))2 fˆ4 (1.5) (0.4055 − fˆ4 (1.5))2
5 f5 (1.5) (0.4055 − f5 (1.5))2 fˆ5 (1.5) (0.4055 − fˆ5 (1.5))2