Notas Matlab Totales

Universidad Autónoma de Sinaloa
Facultad de Ciencias Químico Biológicas

Laboratorio de Simulación de Procesos
Notas
Anabel Altamirano Aguilar

Anafrancella Gutiérrez Montaño
Jesús Raúl Ortiz Del Castillo
Ignacio Calderón Ayala
Notas MATLAB Laboratorio. de Simulación de Procesos
PREFACIO
Las presentes notas surgieron para cubrir la necesidad de una referencia

rápida sobre MATLAB, ya que existen diversas asignaturas en las que se puede
utilizar como herramienta de solución de problemas (MATLAB es una herramienta
de manipulación matemática y numérica muy poderosa). Además de lo anterior,
en las asignaturas de Algebra Lineal (Tronco Común) y Simulación de Procesos
(Ingeniería Química) se incluye explícitamente dentro de su plan de estudios la
utilización de MATLAB. Las presentes notas fueron desarrolladas en el Laboratorio
de Simulación de Procesos como parte del Servicio Social de Anabel Altamirano y
Anafrancella Gutiérrez, ambas alumnas del 8vo. semestre de Ingeniería Química.
Aquí se presenta una introducción al MATLAB, para quienes quisieran
profundizar sobre su utilización existe en el mercado un gran número de libros con
aplicaciones en diversas áreas del conocimiento (tanto en inglés, como en
español). Además, en Internet se pueden encontrar muchas subrutinas gratuitas
que les pueden servir.
Espero que las presentes notas sirvan en su formación profesional y les
ayude a superar los obstáculos que se les presenten para solucionar problemas
en las diversas áreas de la Ingeniería Química. Sólo no se olviden que nadie
aprende a nadar sin mojarse, y si se quiere aprender a utilizar un software hay que
resolver problemas con él. Estoy seguro que con esfuerzo lograran aprovechar el
trabajo que aquí se desarrolló.
Dr. Jesús Raúl Ortiz Del Castillo

Febrero 2006
1
Capítulo 1. Introducción al MATLAB

1. Comandos fundamentales de MATLAB
MATLAB es un sistema interactivo y lenguaje de programación para cómputo

científico y técnico en general. Su elemento básico es la matriz. Dado que los
comandos de MATLAB son similares a las expresiones matemáticas en ingeniería,
escribir soluciones en MATLAB es mucho más fácil que usar un lenguaje de alto
nivel como C o Fortran.
En Macintosh o Windows, haga clic en el icono de MATLAB o de STUDENT

MATLAB. El procedimiento para salir de MATLAB es similar al que sigue para salir
de cualquier otra aplicación en Macintosh o Windows. También para salir puede
teclear en la ventana de comandos: quit o exit.
MATLAB usa tres ventanas de exhibición. La ventana de comandos sirve para

introducir comandos y datos e imprimir resultados, la ventana de gráficos sirve
para exhibir curvas y gráficas, y la ventana de edición que sirve para crear o
modificar archivos M Al ejecutarse los comandos, aparecerán automáticamente las
ventanas apropiadas. Se puede activar una ventana haciendo clic dentro de ella.
Hay varios comandos para despejar ventanas. El comando clc despeja la

ventana de comandos, y el comando clf borra la figura actual y por lo tanto
despeja la ventana de gráficos. El comando clear no afecta a las ventanas, pero sí
borra todas las variables de la memoria. En general, es aconsejable iniciar los
programas con los comandos clear y clf para tener la seguridad de que la memoria
está despejada y que la ventana de gráficos está en blanco y restablecida.
También es importante saber como interrumpir o abortar un comando en

MATLAB. Por ejemplo, puede haber ocasiones en que sus comandos hagan que
la computadora imprima listas aparentemente interminables de números o podría
parecer que la computadora entró en un ciclo infinito. En estos casos, mantenga
presionada la tecla de control y oprima C para generar una interrupción local
dentro de MATLAB.
Si no conoce bien el significado de un comando teclee help y el nombre del

comando en cuestión.
Ejemplo 1:
>> help quit
QUIT Quit MATLAB session.
QUIT terminates MATLAB after running the script FINISH.M,
if it exists. The workspace information will not be saved
unless FINISH.M calls SAVE. If an error occurs while
executing FINISH.M, quitting is cancelled.
2
Nota: En pantalla aparecerá toda la información sobre el comando quit
version. Este comando reporta la versión de MATLAB que se está usando.
what. El comando what produce una lista de los archivos M, MAT y MEX
presentes en el directorio de trabajo actual.
what “nombre de directorio”. Enlista los archivos del directorio en matlabpath.
who. El comando who produce una lista de las variables del espacio de trabajo
actual; whos exhibe información adicional acerca de cada variable; who global y
whos global enlistan las variables del espacio de trabajo actual.
clock. El comando clock exhibe la hora y la fecha con números como:
Ejemplo 2:
>> clock
ans =
1.0e+003 *
2.0060 0.0010 0.0030 0.0190 0.0460 0.0378
Nota: Esto significa [ año, mes, día, hora, minuto, segundo ]
(1.0e+03 * 2.0060 ) = 2006.0 dos mil seis
También con clock podemos medir el tiempo que tarda una ejecución. Por
ejemplo, asigne t_0 = clock antes de que se inicie el cálculo y t_1 = clock cuando
se haya completado; entonces, t_1 – t_0 nos dará el tiempo transcurrido durante
el cálculo.
date. El comando date que proporciona la fecha.
path. El comando path imprime la ruta de búsqueda vigente de MATLAB. El

comando p = path devuelve una cadena p que contiene la ruta. El comando
path(p0) cambia la ruta a p0, que es una cadena que contiene la nueva ruta. El
comando path(p1, p2) cambia la ruta a la concatenación de las dos cadenas de
ruta p1 y p2. Por tanto, path(path, p3) anexará un directorio nuevo p3 a la ruta
vigente y path(p3, path) antepondrá una ruta nueva.
getenv. El comando getenv('MATLABPATH') muestra las rutas de MATLAB

vigentes.
diary. El comando diary on escribe todo lo que se introduce por el teclado, así
como la mayor parte de lo que se envía en pantalla, a un archivo llamado diary.
Con diary off se termina la escritura. Si ya existe una archivo diary, las salidas de
la pantalla se anexarán a ese archivo. Se pueden especificar un nombre de
archivo distinto de diary escribiéndolo después de la palabra diary. Si no se
3
incluyen las palabras on u off, el comando diary alternará entre diary on y dairy off.
El archivo puede imprimirse en papel o editarse posteriormente.
demo. El comando demo guía al usuario para que pueda ejecutar diversas
demostraciones que se eligen de un menú. El contenido de algunas
demostraciones no es fácil de entender a la primera, pero puede estudiarse en
varias ocasiones si se tiene interés.
clc. Si desea borrar la ventana, utilice el comando clc.
2. Cálculos básicos en MATLAB
Si en la ventana de comandos (command window) aparece >> en la esquina

superior, podemos escribir cualquier comando delante de la indicación.
Ejemplo 3:
Evaluemos :
Volumen = 4/3πr3, con r = 2
>> r= 2;
>> vol = (4/3)*pi*r^3;
En el ejemplo anterior observe que cada línea es un comando y termina con un

signo “;”. Cuando trabajamos en la ventana de comandos, la computadora calcula
la respuesta de cada comando inmediatamente después de que pulsamos la tecla
return, pero para evitar esto hay que poner el signo “;” ya que casi nunca resulta
cómodo ir imprimiendo todos los resultados de cada comando que se ejecuta. Si
se quiere obtener el resultado de una operación, por ejemplo, el valor del vol para
nuestro ejemplo anterior puede realizar lo siguiente:
>> vol
vol =
33.5103
Otra forma de imprimir los resultados de cada comando que se ejecute es

separando los comandos con comas y termina la línea con o sin una coma.
Ejemplo 4:
>> r =2, vol =(4/3)*pi*r^3
o
>> r =2
>>vol =(4/3)*pi*r^3
Se imprimirán los valores de r y de vol.
Cuando lo largo de una línea no alcanza para un comando es posible seguir en
otra línea de la ventana de comandos, como se muestra en el Ejemplo 5.
Ejemplo 5:
>> r =2;
4
>>vol =(4/3)*pi
*r^3
La indicación >> no aparecerá en la línea que siga a la marca de continuación.
Operadores aritméticos. Son los mismos que ya conoces como +,-,* y /, pero
MATLAB utiliza uno no tradicional, \, que puede llamarse división inversa. Este
operador produce el recíproco de la división.
Enunciado if. Enunciado condicional y siempre debe de terminar con un enunciado

end.
Operadores utilizados con el if:
> mayor que
< menor que
>= mayor igual que
<= menor igual que
& y
| o
Ejemplo 6:
>>r = 2;
>> if r ~=3, vol = (4/3) * pi *r^3;
>>end
Los operadores & y | se pueden utilizar agrupados.

if((a==2 | b==3) & c<5) g=1; end
El enunciado if se puede utilizar con else o elseif; por ejemplo,
Ejemplo 7:
>> r= 2;
>> if r > 3 b=1;
else if r==3 b=2;
else b=0;
end
Desde luego, elseif puede repetirse tantas veces como lo desee; sin embargo,
hay ocasiones en que el uso de else y elseif tiene sus bemoles, sobre todo cuando
las variables que siguen al enunciado elseif incluyen arreglos de diferentes
tamaños. Si los enunciados elseif no funcionan, olvídese de ellos y repita
enunciados if sencillos cuantas veces sea necesario.
disp. El comando disp exhibe un número, vector, matriz o cadena en la ventana de

comandos sin tener que especificar en nombre de variable; así puede servir para
exhibir mensajes o datos en pantalla. Por ejemplo, tanto disp(pi) como disp pi
5
imprimen 3.14159 en la ventana de comandos. Pruebe también disp “esta es una

prueba de disk”.
2.1 Variables y nombres de variables
Cualquier variable puede adoptar valores reales, complejos y enteros. No es

necesario declarar previamente el tamaño del arreglo. Sin embargo hay que tomar
en cuenta:
- MATLAB no acepta como variable su propio nombre (MATLAB)
- Nombres de ciertos valores
- Nombres de funciones (subrutinas)
- Nombres de comandos
Un ejemplo del conflicto de usar ciertos valores como: sin = 3; ya que “sin” ya no
podrá usarse como función trigonométrica en tanto no sea borrada la variable con
el comando clear o se abandone MATLAB.
En MATLAB i y j se usa para valor unitario imaginario raíz (-1), por eso es
aconsejable no usar i ni j como nombre de variable.
2.2 Ciclos for/end y while/end
Ejemplo 8:
>>for r =1:5
vol = (4/3)*pi*r^3;
disp([r, vol])
end
>>r=0;
>>while r<5
r= r+1;
vol = (4/3)*pi*r^3;
disp ([r, vol)]
end
format. Por omisión, los números se exhiben con cinco dígitos.
Ejemplo 9:
>>pi
ans =
3.1416
Sin embargo, los dígitos pueden exhibirse con 16 dígitos si se emite la orden
format long; por ejemplo,
>>format long
>>pi
6
ans =
3.141592653589793
Con format short se vuelve al formato corto.
break. El comando break termina la ejecución de un ciclo for o while. Hay

ocasiones en que es conveniente utilizar un ciclo infinito que pueda romperse
cuando se satisfaga cierta condición.
Ejemplo 10:
>>for i=1:6
for j =1:20
if j> 2*i, break, end
end
end
clear. Borra todas las variables que se encuentren grabadas en la ventana del
workspace. Sin embargo, si quieres borrar sólo una o algunas variables lo puedes
hacer de la siguiente forma: clear x y z.
2.3 Lectura y escritura
Hay varias formas de pasar datos a, y de, MATLAB. Los métodos pueden
agruparse en tres clases:
a) Operación interactiva mediante teclado o el ratón

b) Lectura o escritura de un archivo de datos
c) Empleo de los comandos save o load
MATLAB puede aceptar datos de entrada a través del teclado mediante el

comando input.
Ejemplo 11:
>>z = input (′indique el radio: ′)
>>z = input (′indique su nombre: ′, ′s′)

′s′ indica que la entrada al teclado es una cadena. Pero si no escribe ‘s’ la cadena
de caracteres debe estar encerrada con apóstrofes (′……′).
fprintf. Con él es posible imprimir mensajes y números con formato si se utiliza.
Ejemplo 12:
>>fprintf ('El volumen de la esfera es %12.5f\n', vol)
\n sirve para que lo que se valla imprimiendo se imprima en reglones separados.
7
Es posible utilizar el enunciado fprintf para escribir salidas con formato de un

archivo. Para ello, se incluye el nombre del archivo en el argumento.
Por ejemplo fprintf('archivo_x', 'volumen =%12.5\n', vol) escribirá la salida en el

archivo de nombre archivo_x. Si no existe el archivo, se creará uno nuevo; si
existe archivo archivo_x, es posible eliminarlo con !rm archivo_x en Unix o !erase
archivo_x en Windows.
Se puede tener mejor control de los archivos con fopen y fclose. Si desea
mayores detalles, consulte la guía de usuario de MATLAB.
3. Variables de arreglo
En muchas ocasiones los datos podrían ser una coordenada en un plano, y se
representa como un par de números, uno de los cuales representa una
coordenada x y el otro la coordenada y. Existen ocasiones en que se puede tener
un conjunto de cuatro coordenadas xyz representando los cuatro vértices de una
pirámide con base triangular en un espacio tridimensional. Es posible representar
todos estos ejemplos usando la estructura de datos Matriz que es un conjunto de
números dispuestos en una retícula rectangular de filas y columnas. Una
coordenada xy puede considerarse una matriz con una fila, de dos columnas y un
conjunto de cuatro coordenadas que puede considerarse como una matriz con
cuatro filas y tres columnas como ejemplo tenemos
A = [3.5] B = [1.5 3.1]
⎡− 1 0 0⎤
C = ⎢⎢ 1 − 1 0⎥⎥
⎢⎣ 0 0 2⎥⎦
Las variables de arreglo unidimensional tienen forma de fila y columna y están

íntimamente relacionadas con los vectores y las matrices. En MATLAB, arreglo de
fila es lo mismo que vector de fila y arreglo de columna es lo mismo que vector de
columna. La variable x puede definirse como vector de fila especificando sus
elementos.
En MATLAB hay varias formas de escribir un arreglo. A continuación se dan

unos ejemplos de estas formas.
Ejemplo 13:
>>x = [ 0, 0.1, 0.2, 0.4, 0.5];
>>for i=1:6
x(i) = (i-1)*0.1;
end
8
>>x = 2: -0.4:-2
Es una forma de describir una variable de arreglo de fila con un incremento o
decremento fijo.
La definición de un arreglo en columna es similar a la de un arreglo de fila

excepto que los elementos se separan mediante signos de punto y coma; por
ejemplo, z = [0; 0.1; 0.2; 0.3; 0.4; 0.5]; para definir esto mismo a partir de un
arreglo de fila es suficiente agregar un apóstrofe: z = [0, 0.1, 0.2, 0.3, 0.4, 0.5]’; el
operador apóstrofe equivales al operador transposición (transpuesta) en el algebra
de marices y vectores, así que convierte vectores de columna en vectores de fila y
viceversa.
Ejemplo 14:
c(8) = 11; si se define un arreglo de esta forma se obtendrá:
>>c(8)=11
c=
0 0 0 0 0 0 0 0 11
Cuando y y x tienen la misma longitud y la misma forma (fila o columna), los

vectores y y x se pueden sumar, restar, multiplicar y dividir (elemento por
elemento) con los siguientes operadores aritméticos:
>>z = x + y
>>z = x – y
>>z = x .* y
>>z = x ./ y
La forma anterior es para ahorrarnos la siguiente forma larga de realizar esas

operaciones
>>for i = 1:6; z(i) = x(i) + y(i); end
>>for i = 1:6; z(i) = x(i) - y(i); end
>>for i = 1:6; z(i) = x(i) * y(i); end
>>for i = 1:6; z(i) = x(i) / y(i); end
El operador de potenciación de arreglos (elemento por elemento) se puede
ilustrar con:
>>g = z .^1.2;
Lo que sería equivalente a:
>>for i = 1:6;
g(i) = z(i) ^1.2;
end
Ejemplo 15:
Si se quiere incrementar el tamaño de un vector debe hacerlo como en el
siguiente ejemplo, teniendo el arreglo x =[2 3] anexar el número 5,
>> x=[2 3];
9
>> x=[x 5]
x=
2 3 5
y para un vector columna se pone “ ;” en vez de “,”.

>> y = [2;3];
>> y = [y;7]
y=
2
7
7
También se puede anteponer un elemento a un vector; por ejemplo

>> x= [9,x]
x=
9 2 3 5
Para consultar el número de filas y columnas de nuestro arreglo se utiliza el

comando length(nombre del arreglo) y sirve para un vector como el siguiente
ejemplo.
Ejemplo 16:
>> x = [9, 2, 3, 5];
>> length(x)
ans =
4
Con el comando size(nombre del arreglo) se obtiene el número de filas y

columnas por la que está compuesto nuestro arreglo.
>> z = [9, 2, 3, 5];
>> size (z)
ans =
1 4
Es decir, z es un arreglo de 1 fila y 4 columnas.
Las variables de cadena en realidad son arreglos. Por ejemplo, una variable de
cadena v definida por v = 'Anabel' equivale a v = ['A' 'n' 'a' 'b' 'e' 'l'] para convertir
en cadena de columna bastaría con aplicar la transpuesta, v = v'.′
Una matriz en MATLAB, se puede definir especificando sus elementos. Por

ejemplo, un arreglo de 3 por 3 se puede definir como el siguiente ejemplo.
Ejemplo 17:
>>m = [1, 2, 3; 5, 6, 7; 2, 3, 5]
m=
1 2 3
10
5 6 7
2 3 5
Para hacer operaciones con arreglos bidimensionales se utilizan los mismos

operadores aritméticos que en operaciones unidimensionales. Obsérvese que es
mucho más sencillo realizar operaciones ya que no se tiene que utilizar ciclos.
4. Operaciones con matrices
La transpuesta de una matriz es una nueva matriz en las que las filas de la
matriz original son las columnas de la nueva.
⎡− 1 3 0⎤ ⎡ − 1 1 2⎤
⎢
A = ⎢ 1 4 7⎥ ⎥ A = ⎢⎢ 3 4 5⎥⎥
T
⎢⎣ 2 5 6⎥⎦ ⎢⎣ 0 7 6⎥⎦
En MATLAB se denota la transpuesta de A como A'.
Ejemplo 18:
>> A=[-1 3 0; 1 4 7; 2 5 6];
>> A'
ans =
-1 1 2
3 4 5
0 7 6
El producto punto es un escalar calculado a partir de dos vectores del mismo

tamaño.
N
A ⋅ B = ∑ ai bi
i =1
En MATLAB se calcula el producto punto empleando la función dot(A,B).
Ejemplo 19:
>> A=[4 -1 3];
>> B=[-2 5 2];
>> dot(A,B)
ans =
-7
En la multiplicación de matrices, el valor que está en la posición ci,j del producto

C de dos matrices, A y B, es el producto punto de la fila i de la primer matriz y la
columna j de la segunda matriz.
N
ci , j = ∑ ai ,k bk , j
i =1
11
Puesto que el producto punto exige que los vectores tengan el mismo número de
elementos, la primera matriz (A) debe tener tantos elementos (N) en cada fila
como elementos hay en cada columna de la segunda matriz (B). Por consiguiente,
si A tiene dos filas y tres columnas, y B tiene tres filas y tres columnas, el producto
AB tendrá dos filas y tres columnas.
⎡1 0 2 ⎤
⎡2 5 1 ⎤
A=⎢ ⎥ B = ⎢⎢− 1 4 − 2⎥⎥
⎣0 3 − 1⎦ ⎢⎣ 5 2 1 ⎥⎦
Ejemplo 20:
>> A=[2,5,1;0,3,-1];
>> B=[1,0,2;-1,4,-2;5,2,1];
>> C=A*B
C=
2 22 -5
-8 10 -7
Si A es una matriz, A.^2 es la operación que eleva al cuadrado cada elemento

de la matriz. Si se quiere elevar al cuadrado la matriz, es decir, si queremos
calcular A*A, se puede utilizar A^2. Así, A^4 equivale a A*A*A*A. Para poder
realizar la multiplicación de dos matrices, el número de filas de la primer matriz
debe ser igual al número de columnas de la segunda; por tanto, para elevar una
matriz a una potencia la matriz debe tener el mismo número de filas que de
columnas; es decir, la matriz debe se cuadrada.
Por definición, la inversa de una matriz cuadrada A es la matriz A-1 para la cual
los productos de las matrices AA-1 y A-1A son iguales a la matriz identidad. Por
ejemplo, considere las siguiente dos matrices, A y B:
⎡ 2 1⎤ ⎡1.5 − 0.5⎤
A=⎢ ⎥ B=⎢
⎣4 3⎦ ⎣− 2 1 ⎥⎦
Si se calculan los productos AB y BA, obtenemos las siguientes matrices:
⎡1 0⎤ ⎡1 0⎤
A=⎢ ⎥ B=⎢ ⎥
⎣0 1 ⎦ ⎣0 1 ⎦
Por lo tanto, A y B son inversos una de las una de la otra, o sea, A = B-1 y B
-1
=A .
En la MATLAB la inversa se calcula con inv(A).
Ejemplo 21:
>> A=[2 1; 4 3];
>> inv(A)
ans =
1.5000 -0.5000
-2.0000 1.0000
12
El rango de una matriz es el número de ecuaciones independientes

representadas por la fila de la matriz. Por tanto, si el rango de una matriz es igual
al número de filas que tiene, la matriz no es singular y existe su inversa.
En MATLAB el rango de una matriz se calcula con rank(A).
Ejemplo 22:
>> A=[2 1; 4 3];
>> rank(A)
ans =
2
Un determinante es un escalar calculado a partir de los elementos de una matriz

cuadrada. Los determinantes tienen varias aplicaciones en ingeniería, incluido el
cálculo de inversas y resolución de sistemas de ecuaciones simultáneas. En
MATLAB el determinante se calcula con det(A).
Ejemplo 23:
>> A=[1 3 0; -1 5 2; 1 2 1];
>> det(A)
ans =
10
5. Aspecto singular de los números en MATLAB
En MATLAB, todas las variables se tratan igualmente con doble precisión.

No hay distinción entre variables enteras y reales, ni entre variables reales y
complejas. La forma en que se asigna un valor depende exclusivamente del
usuario.
La exactitud de los cálculos depende de la forma en que se registran y

procesan los números. Los parámetros clave que indican la exactitud de los
números de un lenguaje de programación son
Número positivo más pequeño: x_min
Número positivo más grande: x_max
Épsilon de la máquina (presición de punto flotante de la máquina): esp
6. Funciones matemáticas en MATLAB
Tabla 1: Número y nombre de variables especiales

Nombre de la Significado Valor
variable
esp Épsilon de la máquina 2.2204E-16
pi π 3.14159...
iyj Unidad imaginaria √-1
13
inf Infinito ∞
NaN No es un número
date Fecha
flops Contador de operaciones de punto flotante
nargin Número de argumentos de entrada de una función
Número de argumentos de salida de una función
nargout
Tabla 2: Funciones matemáticas elementales

Funciones Comentarios
trigonométricas
sin (x)
cos (x)
tan (x)
asin (x)
acos (x)
atan (x)
atan2 (y, x) -π/2 ≥ atan (x) ≥ π/2
Igual que atan(y/x)
-π ≥ atan (y,x) ≥ π
sinh (x)
cosh (x)
tang (x)
asinh (x)
acos (h)
atanh (x)
Otras funciones Comentarios
elementales
abs (x) Valor absoluto de x
angle (x) ángulo de fase de un valor complejo
Si x =real , ángulo = 0
Si x = √-1, ángulo = π/2
sqrt (x) Raíz cuadrada de x
real (x) Parte real del valor complejo x
imag (x) Parte imaginaria del valor complejo x
conj (x) Complejo conjugado x
round (x) Redondear al entero más cercano
fix (x) Redondear un valor real hacia el cero
14
floor (x) Redondear x hacia -∞

ceil (x) Redondear x hacia +∞
sing (x) +1 si x>0; -1 si x<0
rem (x, y) Residuo al dividir: x-y*fix(x/y)
exp (x) Base exponencial e
log (x) Logaritmo de base e
log10 (x) Logaritmo base 10
Existen funciones que realizan tareas.
sort. Ordena los elementos de un vector en orden ascendente. Esto resulta útil en
los casos en que los datos en un orden aleatorio, tienen que acomodarse de una
forma ascendente.
sum. sum(x) calcula la sumatoria de los elementos de un vector o matriz x. Para

los vectores tanto de fila como de columna, sum calcula el total de los elementos.
Si x es una matriz, se calcula la sumatoria de cada columna y se devuelve un
vector de fila formado por las sumatorias de todas las columnas.
Ejemplo 24:
>>sum ([2 1 5]′)
ans = 8
>>sum([2 1 5; 9 8 5])
ans =
11 9 10
max. max(x) encuentra el valor máximo en el vector x, y min(x) encuentra el

mínimo.
rand. Con rand(n) se obtienen una matriz de números aleatorios de tamaño n x n.
7. Creación de un programa en forma de archivo M

La ejecución de comandos en una ventana sólo es apropiada si no hay que
teclear mucho o si se desea explorar ideas en forma interactiva. Sin embargo, en
los casos en que los comandos ocupen varias líneas, es más conveniente que el
usuario escriba un archivo M de guión o de función, porque los archivos M se
pueden guardar en disco y pueden corregirse tantas veces como sea necesario. El
archivo M puede incluir cualquier cosa que el usuario pueda escribir directamente
en la ventana de comandos. Se recomienda a los principiantes tratar de crear
primero archivos M cortos y luego ejecutarlos.
15
Si desea elaborar un archivo M nuevo en Macintosh o Widows, haga clic en New

del menú File (archivo) de la parte superior de la ventana de comandos; aparecerá
una ventana nueva. En Macintosh se hace clic en SAVE and GO del menú.
echo. Cuando se ejecuta un guión, lo normal es que los enunciados del archivo M
no se exhiban en la pantalla. Sin embargo, si se activa el eco con la orden echo
on, los enunciados se exhibirán. De este modo, el usuario puede ver cual parte del
archivo M se esta ejecutando. Para desactivar el eco, teclee echo off.
Enunciados de comentario. El signo % en un archivo M indica que los enunciados

que siguen el signo en la misma línea son comentarios y deben ignorarse durante
los cálculos. Los comentarios que se añaden así a los archivos M pueden ayudar
a explicar el significado de las variables y los enunciados.
Ejemplo 25:
Los números aleatorios pueden servir para crear juegos. El enunciado x =
rand(1) genera un número aleatorio entre 0 y 1 y asigna ese número a x.
Consideremos 13 cartas de espadas que se barajaron bien. La probabilidad de
escoger una carta en particular de la pila es de 1/13. Escriba un programa que
simule la acción de escoger una carta de espadas con un número aleatorio. El
juego debe continuarse devolviendo la tarjeta a la pila y barajándola otra vez
después de cada juego.
Solución:
c=clock;
k=c(2)*c(3)*c(4)*c(5)*c(6);
rand('seed', k)
for k=1:20
n=ceil(13*rand(1)) ;
fprintf ('Número de carta sacada : %3.0f\n', n)
disp (' ')
disp ('Teclee r y pulse Return para repetir')
r = input ('o cualquier otra letra para terminar', 's' );
if r~= 'r' ,break, end

end
Una característica muy importante de los archivos M es que pueden llamar a

otros archivos M. El archivo M que llama es un archivo M padre, en tanto que los
archivos M llamados son archivos M hijos. Esto implica que un guión puede
dividirse en un archivo M padre y varios archivos M hijos.
Ejemplo 26:
Las funciones en MATLAB, que se guardan como archivos M independientes,
equivalen a las subrutinas y funciones de otros lenguajes. Consideremos un
2 x 3 + 7 x 2 + 3x − 1
archivo M de función para la siguiente ecuación: f ( x ) = ,
x 2 − 3x + 5e − x
16
suponiendo que el archivo M se guarda como demos_.m, su guión sería el

siguiente
function y = demos_(x)
y = (2*x.^3+7*x.^2+3*x-1)./(x.^2-3*x+5*exp(-x));
Observe que el nombre del archivo M es idéntico al nombre de la función, que

aparece a la derecha del signo de igual. En el archivo M se utilizan los operadores
aritméticos de arreglos, así que el argumento x puede ser un escalar, un vector o
una matriz. Una vez que se guarda demos_.m como archivo M, se puede utilizar
desde la ventana de comandos o en otro archivo M.
>>y = demos_(3)
y = 502.1384
Si el argumento es una matriz, por ejemplo,

>>demos_ ([3, 1; 0, -1])
y = 502.1384 -68.4920
-0.2000 0.0568
Una función puede devolver más de una variable. Supongamos una función que
evalúa la media y la desviación estándar de una serie de datos. Para devolver las
dos variables utilizamos un vector en el miembro izquierdo del enunciado de la
función; por ejemplo,
Ejemplo 27:
function [media, dvstd] = media_ds(x)
n = length(x);
media = sum(x)/n;
dvstd = sqrt(sum(x.^2)/n – media.^2);
Para utilizar esta función, el miembro derecho del enunciado de llamada también
debe ser un vector. El guión anterior debe guardarse como media_ds.m.
Entonces,
>>x = [1 5 3 4 6 5 8 9 2 4];
[m, d] = media_ds(x)
m = 4.7000
d = 2.3685
Ejemplo 28:
El argumento de una función puede ser el nombre de otra función. Por ejemplo,
supongamos una función que evalúa la media ponderada de una función en tres
f (a ) + 2 f (b ) + f (c )
puntos como: f av = , donde f(x) es la función que se nombrará
4
en el argumento. El siguiente guión ilustra una función f_av.m que calcula la
ecuación anterior.
17
function mp = f_av(nombre_f, a, b, c)
mp = (feval (nombre_f, a) + 2*feval (nombre_f,b) …
+ feval (nombre_f,c))/ 4;
En el guión anterior, nombre_f (una variable de cadena) es el nombre de la

función f(x). Si f(x) es la función seno, nombre_f será 'sin'. feval(nombre_f, x) es un
comando de MATLAB que evalúa la función llamada nombre_f para el argumento
x. Por ejemplo, y = feval('sin', x) equivale a y = sin(x).
Ejemplo 29:
Evalúe la ecuación anterior para la función definida por la ecuación de f con a=1,
b=2, c=3.
Solución:
Suponemos que f_av.m se guardó como archivo M. Entonces, el comando:
>>A = f_av('demof_', 1, 2, 3)
89.8976
El número de entradas y salidas de feval debe coincidir con el formato de la

función “nombre_f”.
La depuración en archivos M de función es más difícil que la de archivos M de

guión. Una de las causas es que no es posible ver los valores de las variables
tecleando los nombres de las variables a menos que se utilicen órdenes de
depuración. El método más básico pero eficaz para crear un archivo M de función
consiste en eliminar la palabra function y probar el archivo M como guión. Cuando
haya depurado exhaustivamente el archivo M, reincorpore el enunciado de la
función. El empleo de comandos de depuración solo se recomienda a usuarios
avanzados de MATLAB.
8. Guardar y cargar datos

Cuando se utiliza el comando save sólo, todas las variables se guardarán
en el archivo por omisión MATLAB.mat. La orden load es el inverso de save y
recupera todas las variables guardadas por save.
Se puede especificar el nombre de archivo colocándolo después de save.

Por ejemplo: save nombre_archivo, guarda todas las variables en el archivo
llamado nombre_archivo.mat. Cuando quiera recuperar las variables escriba load
nombre_archivo. Si sólo desea guardar ciertas variables, escriba sus nombres
después de nombre_archivo; por ejemplo, save nombre_archivo a b c. En este
ejemplo, a, b y c se guardan en el archivo llamado nombre_archivo. No separe
nombre_archivo y las variables con una coma.
Todas las variables se guardan en formato binario de doble precisión. Cuando

quiera cargar los datos contenidos en nombre_archivo.mat, teclee load
nombre_archivo sin nombres de variables; a continuación se recuperaran a, b y c.
18
Se puede utilizar save para escribir datos en formato ASCII. Los comandos load
y save con la opción ASCII son importantes porque permiten exportar e importar
datos en MATLAB: Si desea utilizar el formato ASCII, agregue –ascii o /ASCII
después de los nombres de las variables; save datos.tmp x –ASCII. Guarda la
variable x en ASCII de 8 dígitos en el archivo llamado datos.tmp. El comando save
puede guardar más de una variable.
Ejemplo 30:
>>x = [1, 2, 3, 4];
>>y = [-1, -2, -3]';
>>save dat1.tmp x y –ascii
Si abre el archivo dat1.tmp, se vera así:
1.0000000e+00 2.0000000e+00 3.0000000e+00 4.0000000e+00
-1.0000000e+00
-2.0000000e+00
-3.0000000e+00
El comando load lee un archivo de datos y lo guarda en una variable, pero la

carga de un archivo en formato ASCII no es exactamente el inverso de save en
formato ASCII. La razón es que si bien save en ASCII puede escribir múltiples
variables, load lee todo el archivo de datos y lo coloca en una variable. Además, el
nombre del archivo se convierte en el nombre de la variable. Por ejemplo,
cargamos un archivo llamado y_dat.e con load y_dat.e, el contenido se carga en la
variable llamada y_dat sea cual sea la extensión. Por tanto, el archivo de datos
y_dat debe estar sólo en uno de los siguientes formatos de datos:
(1) un solo vector
(2) un vector de fila
(3) un vector de columna
(4) una matriz
Si tiene necesidad de cargar múltiples variables, cada una deberá prepararse en

un archivo de datos ASCII individual. Los archivos de datos preparados con
Fortran o C en formato ASCII (o de texto) se pueden cargar con load siempre que
la estructura de datos tenga una de las cuatro formas indicadas. Si desea conocer
métodos más avanzados para exportar e importar archivos de datos, consulte la
guía de usuario de MATLAB.
9. Graficación
MATLAB tiene un gran potencial de herramientas gráficas. Se pueden

dibujar los valores de un vector frente a otro (de la misma longitud).
Ejemplo 31:
>> x=pi*(-1:0.1:1);
>> y=x.*sin(x);
>> plot (x, y)
19
Por defecto los puntos (x(i), y(i)) se une por medio de una línea poligonal.
Figura 1. Función seno (Ejemplo 31)
En esta gráfica se ve demasiado lineal a trozos pero si tomamos más puntos se

notará la diferencia:
>> x = pi*(-1:0.01:1);
>> y = x.*sin(x);
>> plot (x, y)
Figura 2. Función seno con mayor número de puntos
Ejemplo 32:
>> fplot ('sin(x)',[0 2 *pi])
Esta expresión dibuja la expresión seno en el intervalo [0, 2*pi]
20
Figura 3. Ejemplo 32
Cuando se quiere mantener en la ventana una gráfica anterior y superponer otra

se utiliza la instrucción hold on. Es una expresión de MATLAB que mantiene en la
ventana gráfica los dibujos anteriores.
Ejemplo 33:
>> fplot ('sin(x) ',[0 2 *pi])
>> hold on
>> fplot ('cos (x) ', [0 2*pi])
Esta expresión dibuja sobre la gráfica anterior la función cos(x), sobre la

gráfica de sin(x).
Con la instrucción hold off se borrar los dibujos anteriores y se grafica en

una ventana nueva.
Ejemplo 34:
>> hold off
>> fplot ('x^2*sin(1/x)',[-0.05 0.05])
21
También puede usarse el versátil comando ezplot que permite dibujar

funciones. Dibuja la función exponencial en un intervalo adecuado a la función, sin
la necesidad de especificar un intervalo.
Ejemplo 35:
>> ezplot ('exp (x)')
Curvas paramétricas.
>> ezplot ('sin (t) ', 'cos(t) ', [0 pi] )
22
Figura 7. Curva paramétrica
Funciones implícitas implícitas.
>> ezplot ('x^2 – y^2 – 1')
Figura 8. Función implícita
También permite dibujar superficies, la forma más sencilla es mediante el

comando ezsurf.
Ejemplo 36:
>> ezsurf ('sin (x*y) ', [-2 2 -2 2]
23
Hasta se pueden realizar gráficas más sofisticadas, por ejemplo:
Ejemplo 37:
>> t=0:0.001:0.009;
>> v=900:1025;
>> [T V]=meshgrid (t,v);
>> aux1=16*pi^2*[T.^2].*((V-918).^2).*((V-1011).^2);
>> aux2=aux1+(2*V-1929).^2;
>> w=T./aux2;
>> z=35000000*w;
>> surfl (t,v,z);
>> shading interp;
>> colormap (pink);
>> rotate3d;
surfl: Este comando dibuja la superficie creada mediante las órdenes anteriores.
Los siguientes comandos sirven para modificar el dibujo obtenido.
rotate3d: Este comando sirve para girar la figura mediante el ratón.
50
40
30
20
10
0
1050
0.01
1000
0.008
0.006
950 0.004
0.002
900 0
24
Capítulo 2. Raíces de Ecuaciones No Lineales

1. Introducción
Las soluciones de una ecuación escalar, f ( x ) = 0 , se llaman ceros o raíces de

f ( x ) . En las presentes notas se estudiarán algunos métodos, como el gráfico el
de la bisectriz, el de iteración de Newton, el de la secante y el de la sustitución
sucesiva, para obtener las raíces reales de ecuaciones no lineales. Si f ( x ) es un
polinomio, podemos utilizar roots, pero incluso en este caso los métodos a
menudo resultan útiles cuando se desea una exactitud muy elevada.
También se verá la aplicación de la sustitución sucesiva y la iteración de Newton

a las ecuaciones simultáneas no lineales.
2. Método Gráfico
Supongamos que deseamos encontrar una raíz positiva de:

f(x) = 0
Donde:
f ( x ) = xsen( 1 ) − 0.2 exp( − x )
x
Si le pedimos a un matemático que encuentre la solución inmediatamente, tal
vez estudiará la ecuación durante unos momentos. Después de percatarse de
que, al aproximarse x a 0 desde el lado positivo de x, sen(1/x) oscila con
frecuencia creciente y se vuelve singular en x = 0, el matemático comenzará a
bosquejar xsen(1/x) y 0.2exp(-x). Después de algunos ensayos, dibujará una
gráfica clara, incluso si lo hace a mano en una hoja de papel en blanco. De la
figura, el matemático deducirá que sólo hay una raíz, misma que se encuentra
cerca de 0.4.
Después de esto, el matemático tal vez utilice un programa de computadora

para encontrar un valor más exacto; sin embargo, la parte más crucial de la
solución se logró con el método gráfico. Sabiendo el comportamiento de la
función, el número de raíces y su valor aproximado, el resto puede hacerse
fácilmente con una computadora o incluso con una calculadora de bolsillo.
La estrategia del matemático es exactamente la que seguiremos con los gráficos

de MATLAB. De hecho, podemos graficar fácilmente xsen(1/x) y 0.2exp(-x) (Figura
1). Como alternativa podemos graficar directamente
f ( x ) = xsen ( 1 / x ) − 0.2 exp( − x ) . Las gráficas indican que la raíz es
aproximadamente 0.36, donde ambas curvas se cruzan. También podemos
obtener un valor más exacto con el método gráfico si amplificamos la gráfica; sin
embargo, sería más eficiente aplicar uno de los métodos numéricos que
analizaremos.
25
0.8
y=x*sin(1/x)
0.6
0.4
Y
X: 0.36
0.2 Y: 0.1395
y=0.2*exp(-x)
-0.2
-0.4
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
X
Figura 1. Resultado del método gráfico
Ejemplo:
Las frecuencias de vibración naturales de una viga uniforme sujeta en un
extremo es la solución de la siguiente ecuación:
cos( x ) cosh( x ) + 1 = 0
Donde:
x = ρω 2 L / EI
L = longitud de la viga (m)
ω = frecuencia( s −1 )
EI = Rigidez ante la flexión (Nm2)
ρ = Densidad del material de la viga (kg/m3)
Determine valores aproximados de las tres raíces positivas más bajas utilizando
el método gráfico.
Solución.
Primero escribiremos
f ( x ) = cos( x ) cosh( x ) + 1
Como sabemos mucho acerca de la función, primero la graficaremos sin límites
de y para 0 ≤ x ≤ 20 . Esta figura nos dice que una raíz es aproximadamente
x=17.4, pero que también es posible que existan otras raíces en 0 ≤ x ≤ 15 .
Listado 1
clf; clear
x = 0:0.1:20;
y = cos(x).*cosh(x)+1;
plot(x, y, x, zeros(size(x)), '--');
xlabel('x'); ylabel ('y = cos (x)*cosh(x)+1' )
26
7
x 10
12
10
y = cos (x)*cosh(x)+1
6
2
X: 17.4
Y: 0
0
-2
0 2 4 6 8 10 12 14 16 18 20
x
Figura 2. Resultado del ejemplo 1

Sin borrar la figura trazada por el guión anterior, podemos modificar los límites
de la gráfica con el comando axis([0 20 -10 20]) cambia la gráfica a la que se
mostró.
Se encontraron tres raíces positivas más pequeñas que son x = 1.8, 4.6 y
7.8, aproximadamente.
Como es evidente por el ejemplo que acabos de describir, el método gráfico

no carece de dificultades ni peligros. Un problema es que en algunos casos es
probable que una gráfica tenga una presentación muy deficiente. Por ejemplo, la
calidad de una gráfica disminuye considerablemente en las inmediaciones de una
singularidad. Si no tenemos cuidado al trazar una gráfica, podrían quedar ocultas
características especiales, como puntos singulares y en ocasiones confundirse
con raíces. Si se sospecha que existe una singularidad amplifique la gráfica y
determine si la función es realmente singular. Como ejemplo adicional, si se
grafica una función con intervalos equiespaciados, es posible que no se capture
una oscilación rápida y que la curva trazada difiera significativamente de la función
verdadera. Por tanto, es aconsejable graficar la función varias veces con
diferentes amplificaciones y focos hasta entenderla bien.
3. Método de la Bisectriz
El método de la Bisectriz es un método sencillo pero muy versátil para

encontrar una raíz real en un intervalo dado en el que se sabe que existe una raíz.
Su singular ventaja consiste en que funciona incluso con funciones no analíticas;
sin embargo, sólo debe utilizarse el método después de un análisis de un gráfico.
Supongamos que hay una raíz de f ( x ) = 0 en un intervalos de x = a y
x = c , denotado por [ a, c ] o de forma equivalente, por a ≤ x ≤ c . El método de la
Bisectriz se basa en el hecho de que, cuando un intervalo [ a, c ] tiene una raíz, el
signo de y (x) en los dos extremos es distinto, o sea f ( a ) f (c ) < 0 . El primer paso
de este método consiste en bisecar el intervalo [ a, c ] en dos mitades: [ a, b ] y [ b,
27
c ], donde b = ( a + c ) / 2 . Si verificamos los signos f (a ) f (b) y f (b) f (c ) , sabremos

en que mitad del intervalo se encuentra la raíz. De hecho, si f ( a ) f (b) ≤ 0 ,
sabremos que el intervalo [ a, b ], que incluye x = a y x = b , contiene a la raíz; de
lo contrario, la raíz estará en el otro intervalo, [ b, c ]. A continuación, bisecamos
de nuevo el intervalo que contiene a la raíz. Al repetir este procedimiento el
tamaño del intervalo que contiene a la raíz se hará cada vez más pequeño. En
cada paso se toma el punto medio del intervalo como aproximación más
actualizada a la raíz la iteración se detiene cuando el tamaño de la mitad del
intervalo es menor que una tolerancia dada.
El tamaño del intervalo después de n pasos de iteracción es:
c0 − a0
2n
Donde a0 y c0 son valores iniciales a a y c , de modo que el numerador es el

tamaño de intervalo inicial. La ecuación anterior también representa el máximo
error que puede haber cuando se aproxima la raíz con el n-ésimo punto medio.
Por tanto, si la tolerancia del error es τ , el número de pasos de iteracción es el
entero n más pequeño que satisface
c −a
τ ≥ 0 n 0
2
O de forma equivalente,
⎛ c − a0 ⎞
log⎜ 0 ⎟
⎝ τ ⎠
n≥
log(2 )
Donde τ es la tolerancia. Por ejemplo, si c0 – a0 = 1 y τ = 0.0001 , entonces n = 14.

Podemos usar dos funciones, bisec_g y bisec_n, para cálculos de bisección.
La primera muestra gráficamente el avance de la iteracción de bisección. La
segunda función no traza gráficas, pero es más rápida.
La sintaxis de bisec_g es la siguiente:
Bisec_g ( ‘ nombre_f ‘, a, c, xmin, xmax, n_puntos )
Donde nombre_f es el nombre de la función que define la ecuación por resolver, a

y c son los extremos del intervalo inicial, xmin y xmax son los valores mínimos y
máximos de x en la gráfica y n_puntos es el número de puntos que se usarán para
graficar la función la tolerancia es τ = 10 −6 por omisión.
La sintaxis de bisec_n es:
28
bisec_n ( ‘ nombre_f ‘, a, c )
Ejemplo
Encuentre la intersección de estás dos funciones:
y = x2 + 1
y = tan( x ),0 < x < π / 2
Solución
El problema equivale a encontrar el cero de
f = x 2 + 1 − tan( x),0 < x < π / 2
Si graficamos la función, veremos que la raíz es aproximadamente x = 0.9. Para

poder utilizar bisec_n (o bisec_g), escribimos un archivo M de función que defina
la ecuación primera.
function f = fun_ex2 (x)

f = sqrt (1+x.^2) – tan (x) ;
Luego, utilizamos bisec_n así:
Bisec_n ( ‘ fun_ex2 ‘, 0.8, 1.0 )
La salida es:
Método de la Bisectriz
It. a b c f (a) f (b ) f (c)
1 0.800000 0.900000 1.000000 0.250986 0.085204 -0.143194
2 0.900000 0.950000 1.000000 0.085204 -0.019071 -0.143194
3 0.900000 0.925000 0.950000 0.085204 0.035236 -0.019071
4 0.925000 0.937500 0.950000 0.035236 0.008660 -0.019071
5 0.937500 0.943750 0.950000 0.008660 -0.005056 -0.019071
6 0.937500 0.940625 0.943750 0.008660 0.001838 -0.005056
7 0.940625 0.942187 0.943750 0.001838 -0.001600 -0.005056
8 0.940625 0.941406 0.942187 0.001838 0.000122 -0.001600
9 0.941406 0.941727 0.942187 0.000122 -0.000738 -0.001600
10 0.941406 0.941602 0.941797 0.000122 -0.000308 -0.000738
11 0.941406 0.941504 0.941602 0.000122 -0.000093 -0.000308
12 0.941406 0.941455 0.941504 0.000122 0.000014 -0.000093
13 0.941455 0.941479 0.941504 0.000014 -0.000040 -0.000093
14 0.941455 0.941467 0.941479 0.000014 -0.000013 -0.000040
15 0.941455 0.941461 0.941467 0.000014 0.000001 -0.000013
29
16 0.941461 0.941464 0.941467 0.000001 -0.000006 -0.000013

17 0.941461 0.941463 0.941464 0.000001 -0.000003 -0.000006
18 0.941461 0.941462 0.941463 0.000001 -0.000001 -0.000003
19 0.941461 0.941462 0.941462 0.000001 0.000000 -0.000001
Se satisface la tolerancia
Resultado final : Raíz = 0.941462
Archivo M para bisec_g
function bisec_g(f_name, a, c, xmin, xmax,n_points)

f_name
% a, c : extremos del intervalo inicial
% tolerance : tolerancia
% it_limit : limite del numero de iteraciones
% Y_a, Y_c : valores y de los extremos actuales
% fun_f(x) : valor funcional en x
clf, hold off
clear Y_a, clear Y_c
wid_x = xmax - xmin; dx =(xmax - xmin)/n_points;
xp=xmin:dx:xmax; yp=feval(f_name, xp);
plot (xp, yp); xlabel ('x') ; ylabel ('f(x)');
title('Metodo de la bisectriz'), hold on
ymin=min(yp); ymax=max(yp); wid_y = ymax-ymin;
yp=0.*xp; plot(xp,yp)
fprintf('Metodo de la bisectriz:\n\n');
tolerance = 0.000001; it_limit = 30;
fprintf(' It. a b c fa=f(a) ');
fprintf('fc=f(c) abs(fc-fa) \n');
it = 0;
Y_a = feval (f_name, a); Y_c = feval (f_name, c);
plot([a,a], [Y_a,0]); text(a,-0.1*wid_y,'x=a')
plot([c,c], [Y_c,0]); text(c,-0.1*wid_y,'x=c')
if (Y_a*Y_c > 0) fprintf(' f(a)f(c) > 0 \n');
else
while 1
it = it + 1;
b = (a + c)/2; Y_b = feval(f_name, b);
plot([b,b], [Y_b,0],':'); plot(b,0,'o')
if it<4, text(b, wid_y/20, [num2str(it)]), end
fprintf('%3.0f, %10.6f, %10.6f', it, a, b);
fprintf('%10.6f, %10.6f, %10.6f', c, Y_a, Y_c);
fprintf(' %12.3e\n', abs((Y_c - Y_a) ));
if ( abs(c-a)<=tolerance )
fprintf('Se satisface la tolerancia. \n' ); break
end
if ( it>it_limit )
30
fprintf('Se excedio limite de iteraciones.\n'); break

end
if ( Y_a*Y_b <= 0 ) c = b; Y_c = Y_b;
else a = b; Y_a = Y_b;
end
end
fprintf('Resultado final: Raiz = %12.6f \n', b);
end
x=b;
plot([x x], [0.05*wid_y 0.2*wid_y])
text(x, 0.25*wid_y, 'Solucion final')
plot ([x (x-wid_x*0.004)], [0.05*wid_y 0.09*wid_y])
plot ([x (x+wid_x*0.004)], [0.05*wid_y 0.09*wid_y])
Archivo M para bisec_n
function bisec_n(f_name, a, c)
f_name
% a, c : extremos del intervalo inicial
% tolerance : tolerancia
% it_limit : limite del numero de iteraciones
% Y_a, Y_c : valores y de los extremos actuales
% fun_f(x) : valor funcional en x
fprintf('Metodo de la bisectriz:\n\n');
tolerance = 0.000001; it_limit = 30;
fprintf(' It.a b c fa=f(a) ');
fprintf(' fc=f(c) abs(fc-fa) \n');
it = 0;
Y_a = feval(f_name, a); Y_c = feval (f_name, c);
if (Y_a * Y_c > 0)
fprintf('\n \n Detenido porque f(a)f(c) > 0 \n');
else
while 1
it = it + 1;
b = (a + c)/2; Y_b = feval(f_name, b);
fprintf('%3.0f %10.6f, %10.6f', c, Y_a, Y_c);
fprintf(' %12.3e\n', abs((Y_c - Y_a)));
if ( abs(c-a)/2<=tolerance )
fprintf(' Se satisface la tolerancia. \n' ); break
fprintf(' Cambie a o b y ejecute otra vez .\n' );
end
if ( Y_a*Y_b <= 0 ) c = b; Y_c = Y_b;
else a = b; Y_a = Y_b;
end
end
fprintf('Resultado final : Raiz = %12.6f \n', b);
end
31
4. Iteración de Newton
La iteración de Newton es un método iterativo para encontrar la raíz de una

ecuación no lineal. Puede aplicarse en el dominio complejo para encontrar una
raíz compleja y también puede extenderse a ecuaciones simultáneas no lineales.
La iteración de Newton se deduce de la expansión de Taylor. Supongamos

que el problema consiste en encontrar una raíz de f ( x ) = 0 . La expansión de
Taylor truncada de primer orden de f (x ) alrededor de una estimación inicial, x0,
se escribe así:
f ( x) = f ( x0 ) + f ' ( x0 )(x − x0 )
Que se considera como una aproximación de f ( x ) . Si igualamos la ecuación a

cero obtenemos la siguiente aproximación:
f ( x0 )
x1 = x0 −
f ' ( x0 )
Repetimos el mismo proceso con:
f ( xn −1 )
xn = xn +1 −
f ' ( xn −1 )
Para el valor inicial de x0, se traza la línea que pasa tangencialmente por (x0,
f0). La intersección de la línea tangencial con el eje x es x1. A continuación, se
traza la línea que pasa tangencialmente por (x1, f1). Este procedimiento se repite
utilizando el valor más actualizado como estimación para el siguiente ciclo de
iteración.
La deducción de la primera derivada de una función dada podría ser
laboriosa o incluso imposible. En un caso así, podría evaluarse f ' (x ) con una
aproximación de diferencia en lugar de analíticamente. Podemos aproximar
f ' ( xn −1 ) con:
f ( xn −1 + h ) − f ( xn −1 )
f 'n −1 =
h
O bien
f ( xn −1 ) − f ( xn −1 − h )
f 'n −1 =
h
Donde h es un valor pequeño como h = 0.001. Las ecuaciones anteriores son las
aproximaciones de diferencia hacia delante o hacia atrás, respectivamente.
32
Errores pequeños en la aproximación de diferencia no tienen un efecto perceptible

en la rapidez de convergencia de iteración de Newton. La exactitud el resultado
final no se ve afectado por el error de una aproximación de diferencia, sin
embargo, en los casos en que no hay singularidad cerca de la raíz hay que tener
mucho cuidado al utilizar la aproximación de diferencia.
Ejemplo
Deduzca un esquema iterativo para encontrar la raíz cúbica de un número

con base en la iteración de Newton. Encuentre la raíz cúbica de a = 155 con el
esquema que deduzca:
Solución:
El problema consiste en encontrar el cero de
f (x ) = x3 − a
La iteración de Newton se escribe así:
f ( xn )
xn −1 = xn −
f ' (xn )
xn3 − a
= xn −
3 xn2
2 a
= xn + 2
3 3 xn
Asignamos a = 155 y utilizamos como estimación inicial x0 = 5. La iteración

procede así:
n x
0 10
1 5.4
2 5.371834
3 5.371686 (exacta)
La solución exacta se obtiene después de apenas tres pasos de iteración. Lo

intentamos otra vez con una estimación inicial mucho más alejada del valor real, x0
= 10:
n x
0 10
1 7.183334
2 5.790176
3 5.401203
33
4 5.371847
5 5.371686 (exacta)
Se obtiene el valor exacto de la raíz cúbica con cinco pasos de iteración.
Dos funciones, Newt_g y Newt_n, resuelven ecuaciones no lineales por

iteraciones de Newton. La primera exhibe gráficamente el avance de la iteración,
mientras que el segundo sólo realiza los cálculos.
Las sintaxis son:
Newt_g ( ‘ nombre_f ‘ , x0, Xmin, Xmax, n_puntos )

Newt_n ( ‘ nombre_f ‘ , x0 )
Donde nombre_f es el nombre del archivo M de función que define la

ecuación por resolver y x0 es una estimación inicial de la raíz. El significado de
xmin, xmax y n_puntos es el mismo que tienen en bisec_g.
Resolvemos la siguiente ecuación:
y = (0.01x + 1)sen( x ) − (x − 0.01)

(x 2
+1 ) − 0.0096
Esta ecuación se escribe en un archivo M y lo llamaremos: eqn_1.m luego
ejecutamos newt_n ( ‘eqn_1’,4 ) la salida es:
Iteración de Newton
Teclee el nombre de la función (encerrada en apóstrofos):
‘eqn_1’
f_name =
eqn_1
Teclee la estimación inicial de la raíz: 4
n = 1, x = 2.36780E+00, Y = -1.030E+00, yd = -6.319E-01
n = 2, x = 2.92631E+00, Y = 3.488E-01, yd = -6.247E-01
n = 3, x = 2.82370E+00, Y = -9.467E-02, yd = -9.226E-01
n = 4, x = 2.82171E+00, Y = -1.774E-03, yd = -8.895E-01
n = 5, x = 2.82170E+00, Y = -4.498E-06, yd = -8.888E-01
n = 6, x = 2.82170E+00, Y = -9.553E-09, yd = -8.888E-01
La respuesta final: 2.8217E+00
La salida de Newt_g es idéntica a la de Newt_n, sólo que la primera traza una

gráfica de avance.
Newt_g ( ‘eqn_1’, 4, 0, 7, 50 )
Ejemplo:
34
Imagine una pared de tabique con un espesor de 0.05 m. La temperatura en

el lado interior de la pared, T0, es de 625 K, pero se desconoce la temperatura en
el lado exterior. La pérdida de calor de la superficie exterior se efectúa por
convección y por radiación.
La temperatura T1 esta determinada por la siguiente ecuación:
f (T1 ) =
k
Δx
( )
(T1 − T0 ) + εσ T14 − T∞4 + h(T1 − T f ) = 0
Donde:
k = Conductividad térmica de la pared 1.2W

mK
ε = Emisividad, 0.8
T0 = Temperatura del lado interior de la pared, 625 K
T1 = Temperatura del lado exterior de la pared (desconocida K)
T∞ = Temperatura del entorno, 298 K
Tf = Temperatura del aire, 298 K
h = Coeficiente de transferencia de calor 20W 2
m K
−08 W
σ = Constante de Stefan-Boltzamann, 5.67 x10
m2 K 4
Δx = Espesor de la pared, 0.05 m
Determine T1 por iteración de Newton
Solución:
Se resolverá el problema con Newt_g la ecuación por resolver se escribe en un

archivo M de función como se muestra.
Archivo M
Function f = cal_pared (T1)
K = 1.2; e = 0.8; Tinf =298 K;
Tf = 298; h = 20; T0 = 625 K;
sig = 5.67 E − 8 ; espesor = 0.05;
f = k/espesor * (T1 –T0 )+e*sig*( T1.^4 - Tinf ^4)…
+h*( T1 –Tf );
Después de guarduar el archivo M de la función anterior, ejecutamos el siguiente

comando:
Newt_g (‘cal_pared ‘, 550, 400, 600, 500 )
El resultado es:
35
Iteración de Newton:
n = 0, x = 5.50000E+02, Y = 7.033001E+03
n = 1, x = 4.55199E+02, Y = 6.58551E+02
n = 2, x = 4.44423E+02, Y = 6.44623E+00
n = 3, x = 4.44316E+02, Y = 6.27680E-04
n = 4, x = 4.44316E+02, Y = 5.70253E-10
ans =
444.317
La respuesta final T1 = 444.317 K. La salida gráfica, donde X y Y denotan T1 y f(T1

), respectivamente.
Archivo M para Newt_g
function x = Newt_g(f_name, x0, xmin, xmax, n_points)

clf, hold off
% Metodo de Newton con ilustracion grafica
del_x = 0.01;
wid_x = xmax-xmin; dx = (xmax-xmin)/n_points;
xp=xmin:dx:xmax; yp=feval(f_name, xp);
plot(xp, yp); xlabel('x'); ylabel('f(x)');
title('Iteracion de Newton'), hold on
ymin=min(yp); ymax=max(yp); wid_y = ymax-ymin;
yp=0.*xp; plot(xp,yp)
x = x0; xb=x+999; n=0;
while abs (x-xb)>0.000001
if n>300 break; end
y=feval(f_name, x); plot([x,x], [y,0]); plot(x,0, 'o')
fprintf(' n=%3.0f, x=%12.5e, y=%12.5e\n', n,x,y);
xsc = (x-xmin)/wid_x;
if n>4, text(x, wid_y/20, [num2str(n)]), end
y_driv = (feval(f_name, x+del_x) - y)/del_x;
xb=x;
x = xb, y/y_driv; n=n+1;
plot([xb,x], [y,0])
end
plot([x x], [0.05*wid_y 0.2*wid_y])
text(x, 0.2*wid_y, 'Solucion final')
plot([x (x-wid_x*0.004)], [0.01*wid_y 0.09*wid_y])
plot([x (x+wid_x*0.004)], [0.01*wid_y 0.09*wid_y])
Archivo M para Newt_n
function x = Newt_n(f_name, x0)

% Iteracion de Newton sin graficos
x = x0; xb = x-999;
36
n = 0; del_x = 0.01;
while abs(x-xb)>0.000001
n=n+1; xb=x;
if n<300 break; end
y=feval(f_name, x);
y_driv=(feval(f_name, x+del_x) - y)/del_x;
x = xb - y/y_driv;
fprintf(' n=%3.0f, x=%12.5e, y=%12.5e,', n, x, y)
fprintf(' yd = %12.5e \n', y_driv)
end
fprintf('\n Respuesta final = %12.6e\n', x);
5. Método de la Secante
El método de la secante es una variante de la iteración de Newton. Hemos

utilizado una aproximación de diferencia para evaluar f ’ en la iteración de Newton;
sin embargo, también podemos evaluar f ‘ aproximadamente utilizando los dos
valores anteriores consecutivos de f. El esquema iterativo basado en este
concepto se escribe así:
xn −1 − xn − 2
xn = xn −1 − f n −1 , n = 1,2,3...
f n −1 − f n − 2
Para iniciar la iteración, tenemos que especificar x0. El valor de x1 puede

fijarse arbitrariamente en x1 = x1 + Δx, donde Δx es un número arbitrariamente
pequeño, como 0.01, por ejemplo. Entonces, podemos continuar la iteración hasta
satisfacer una tolerancia.
Ejemplo
Una bala de M = 0.002 kg. Se disparó verticalmente en el aire y está

descendiendo a su velocidad terminal. La velocidad terminal esta dada por
gM=Farrastre, donde g es la aceleración debida a la gravedad y M es la masa.
Después de evaluar todas las constantes, la ecuación se escribe así:
(0.002)(9.81) = 1.4 x10−05 v1.5 + 1.15x10−5 v 2

Donde v es la velocidad terminal, m/s. El primer término del miembro
derecho representa el arrastre por fricción y el segundo, el arrastre por presión.
Determine la velocidad terminal por el método de la secante. Una estimación
burda sería v = 30 m/s.
37
Solución
La tarea consiste en encontrar la raíz de:
f (v ) = (0.002)(9.81) − 1.4 x10−5 v1.5 + 1.15x10−5 v 2
Asignamos v0 = 30 y v1 = 30.1, con base en la estimación burda dada y

calculamos f0 y f1 con la ecuación. La solución iterativa:
n v f(v)
0 30.00000 1.9620001E-02
1 30.10000 6.8889391E-03
2 30.15411 6.8452079E-03
3 38.62414 -8.9657493E-04
4 37.64323 9.0962276E-05
5 37.73358 9.9465251E-07
6 37.73458 -1.8626451E-09
Así, la velocidad Terminal es v = 37.7 m/s.
6. Método de Sustituciones Sucesivas
El método de sustituciones sucesivas se refiere a una clase amplia de

esquemas de resolución iterativos para ecuaciones no lineales. La iteración de
Newton y el método de la secante pueden considerarse aplicaciones de la
sustitución sucesiva. Puesto que la sustitución sucesiva se emplea en muchos
algoritmos numéricos para resolver ecuaciones no lineales, incluidas ecuaciones
diferenciales y ecuaciones simultáneas no lineales, presentaremos algunos de sus
aspectos fundamentales en esta sección.
Si la ecuación por resolver, f(x)=0 se puede escribir en la forma:
x = g (x )
Entonces podrá escribirse un esquema iterativo así:
xn = g ( xn −1 )
Donde n es el número de pasos de iteración y x0 es una estimación inicial.

Este método se denomina de sustituciones sucesivas o de iteración de punto fijo.
La ventaja de este método es su sencillez y la flexibilidad para escoger la
forma de g(x) elegida arbitrariamente. Si queremos asegurar la convergencia de la
iteración, deberemos satisfacer la siguiente condición:
38
g ' ( x) < 1
Ejemplo:
Se sabe que la función
y = x 2 − 3x + e x − 2
Tiene dos raíces, una negativa y una positiva. Encuentre la raíz pequeña por
sustitución sucesiva.
Solución
Si verificamos el signo de y en x = −1 y x = 0 (a saber, y (− 1) = 2.367 y y (0 ) = −1 ),

localizaremos la raíz más pequeña en -1< x <0. Reescribimos la ecuación dada
así:
x2 + ex − 2
x = g (x ) =
3
Después, escribimos el siguiente esquema iterativo:
xn = g ( xn −1 )
La primera derivada de g(x) satisface la ecuación anterior en el intervalo -1<

x < 0, así que el esquema es convergente.
Las siguientes ecuaciones alternativas:
x = − 3x − e x + 2
Y
x = 3x − e x + 2
Sin embargo, las ecuaciones anteriores tienen discontinuidades en las
inmediaciones de la raíz más pequeña. Además que las primera derivadas de
ambas ecuaciones violan la condición de formulas anteriores. Por tanto, ninguna
de estas ecuaciones funciona.
Una forma sistemática de encontrar una forma de g(x) consiste en

establecer:
g ( x ) = x − αf ( x )
De modo que el esquema iterativo se convierte en:
xn = xn +1 − αf ( xn +1 )
39
Donde α es una constante. La constante α se puede determinar como sigue: al

sustituir la primera ecuación anterior en la segunda que acabamos de escribir y se
puede observar que la iteración converge cuando:
− 1 < 1 − α f ' ( x) < 1
O lo que es lo mismo:
0 < α f ' ( x) < 2
La ecuación anterior indica que, primero debe tener el mismo signo que f’ y
segundo, que la rapidez de convergencia es óptima cuando α ≈ 1 / f ' .
La presente esquema se reduce a la iteración de Newton si se asigna al
valor de 1 / f ' ( xn ) para cada iteración.
Ejemplo
El tamaño crítico de un reactor nuclear está determinado por una ecuación

de criticalidad. Suponga que una versión simplificada de la ecuación de criticalidad
está dada por:
tan(0.1x) = 9.2e− x
La solución que tiene significado físico es la raíz positiva más pequeña que
satisface 3 < x < 4. Determine la raíz positiva más pequeña.
Solución
Aplicamos el esquema iterativo de una de las ecuaciones anteriores a la siguiente

función:
f ( x) = tan(0.1x) − 9.2e− x
Estimamos un valor aproximado de f’ en 3 < x < 4 con:
f ( 4) − f (3)
f '= = 0.40299
4−3
Luego, asignamos el siguiente valor al parámetro:
1
α=
0.40299
La iteración de la ecuación definida anteriormente converge como sigue:
40
Número de Solución
Iteración iterativa
n x
0 4.00000
1 3.36899
2 3.28574
3 3.29384
4 3.28280
5 3.29293
6 3.29292
7 3.29292
7. Ecuaciones Simultáneas No Lineales
La necesidad de resolver ecuaciones simultáneas no lineales se presenta

con bastante frecuencia. Aquí presentaremos dos métodos para resolver
ecuaciones simultáneas no lineales.
Iteración con sustituciones sucesivas.
Si el sistema de ecuaciones no lineales dado representa el fenómeno natural

o un sistema de ingeniería, es común que las ecuaciones se vuelvan lineales si la
magnitud de la solución es pequeña un sistema no lineal de este tipo puede
escribirse en la misma forma que el sistema lineal excepto que los coeficientes
dependen de la solución.
La resolución iterativa de un sistema no lineal con base en sustituciones

sucesivas puede escribirse así:
An −1 X n = y
Donde An −1 es una matriz que representa los coeficientes y que se calcula a

partir de la solución anterior; X n es la n-ésima solución iterativa y y es un término
no homogéneo que también puede ser función de la solución.
Inicialmente, la matriz de coeficientes se calcula utilizando una estimación

inicial de la solución. Una vez determinada la matriz, la ecuación se resuelve como
un sistema lineal. Ya que se obtiene la solución, la matriz de coeficientes se
actualiza y se resuelve una vez más la ecuación. Si se presenta una inestabilidad
durante la iteración, utilizamos sub-relajación:
X n = ωAn−−11 y + (1 + ω ) X n −1
Donde ω de un parámetro de sub-relajación que satisface 0 < ω <1.
41
Capítulo 3. Resolución de Ecuaciones Diferenciales

Parciales
Introducción
Debido a las enormes aplicaciones prácticas del método de elementos finitos

existen multitud de paquetes informativos que simplifican su aplicación a distintos
problemas de ingeniería. Se presenta aquí la herramienta de ecuaciones en
derivadas parciales (PDE) de Matlab.
Esta Toolbox cuenta con una interfaz gráfica de Matlab pdetool, parte de la PDE
toolbox proporciona una herramienta gráfica de fácil manejo para la descripción de
estas geometrías complicadas, generación de mallas, resolución de la ecuación
discretizadas y representación de resultados. Además, muchas de las aplicaciones
de interés en ingeniería tienen formas concretas y al ingeniero tan solo le interesa
el comportamiento de un determinado material o geometría en su problema de
estudio, sin tener que preocuparse de la formulación matemática inherente al
mismo. Por este motivo pde tool de Matlab incorpora varios modos de solución
específicos para la solución de ciertos problemas típicos en ciencia e ingeniería,
de modo que el trabajo del usuario tan sólo se limita a definir las geometrías y la
introducción de los parámetros del material.
Esta la Toolbox tiene un hermano mayor que es el programa FEMLAB, que es un

programa independiente con mayores posibilidades que ésta pero que esta
basado en los mismos conceptos que ella.
La interfaz gráfica de la PDE toolbox
“La pde toolbox de Matlab incluye una interfaz de usuario completa, que cubre
todos los aspectos del proceso de solución de una EDP. La interfaz se inicializa
escribiendo:
>> pdetool
En la línea de comandos de Matlab, siempre y cuando la toolbox mencionada se

encuentre instalada en tu sistema.
Menús
Existen un total de 11 menús desplegables en la interfaz gráfica. Brevemente, la

funcionalidad de cada uno de ellos es la siguiente:
• File: Como es habitual, desde este menú pueden abrirse y salvarse ficheros
.m que contienen los modelos en ecuaciones en derivadas parciales en los
42
que se esté trabajando. También pueden imprimirse las graficas activas en

ese momento y salir de la interfaz.
• Edit: Capacidades de edición habituales: copiar, cortar, borrar y pegar
objetos, así como opciones de seleccionar todo.
• Options: Contiene opciones como cambiar el rango y espaciado de los ejes,
obligar a que en la fase de dibujo las formas se anclen a los puntos de la
rejilla, zoom, etc.
• Draw: Desde este menú se pueden seleccionar los objetos sólidos básicos
como círculos o polígonos que se emplearán en la definición de la
geometría y a continuación dibujarlos en el área de trabajo mediante el uso
de ratón. Se recomienda el uso de la barra de herramientas para este fin.
• Boundary: Desde este menú se accede al cuadro de diálogo donde se
definen las condiciones de frontera. Adicionalmente se pueden poner
etiquetas a los bordes y a los subdominios, borrar bordes entre
subdominios y exportar la geometría descompuesta y las condiciones de
fronteras al espacio de trabajo de Matlab.
• PDE: Este menú proporciona un cuadro de diálogo para especificar la EDP
y opciones para etiquetar subdominios y exportar los coeficientes de la
ecuación al espacio de trabajo.
• Mesh: Desde este menú se crea y se modifica la malla triangular. Puede
inicializarse la malla, refinarla, reorganizarla, deshacer cambios previos en
la malla, etiquetar los nodos y los triángulos, visualizar la calidad de la malla
y exportar la misma al espacio de trabajo.
• Solve: Para resolver la EDP. También abre un cuadro de diálogo donde
pueden ajustarse los parámetros involucrados en la resolución y exportar la
solución al espacio de trabajo.
• Plot: Desde este menú se puede dibujar correctamente la solución a la
EDP. Un cuadro de diálogo permite seleccionar que propiedad va a
visualizarse, en que estilo y otros tipos de parámetros. Si se ha generado
una animación en tiempo de la solución, también puede exportarse al
espacio de trabajo.
• Window: Básicamente para seleccionar cuál de las ventanas de figuras de
Matlab es la activa entre las abiertas en ese momento.
• Help: Breve ayuda sobre ciertos comandos y funcionalidades.
Barra de Herramientas
La barra de herramientas colocada debajo del menú principal en la parte superior

de la interfaz gráfica contiene los botones con iconos que proporcionan un acceso
fácil y rápido a algunas de las funcionalidades más importantes de la pdetool.
Barra de Herramientas de la pdetool de Matlab.
43
Los cinco botones de la izquierda son los botones del modo de dibujo,
representando cada uno de izquierda a derecha:
• Dibuja un rectángulo/cuadrado comenzando en una esquina.
• Dibuja un rectángulo/cuadrado comenzando en el centro.
• Dibuja una elipse/círculo comenzando en el perímetro.
• Dibuja una elipse/círculo comenzando en el centro.
• Dibuja un polígono. Pulsa y arrastra para crear los lados del

polígono. El polígono puede cerrarse haciendo clic con el botón derecho del
ratón o pulsando sobre el vértice inicial.
Como en gran número de programas informáticos, sólo puede haber uno de estos
botones de dibujo activado en cada momento. El hacer doble clic sobre uno de
ellos fija esa herramienta como activa, pudiendo seguir dibujando objetos del
mismo tipo hasta que vuelva a pulsarse el botón. Usando el botón derecho del
ratón o bien Control + Clic, se restringen las herramientas a dibujar cuadrados o
círculos en vez de rectángulos o elipses.
El segundo grupo de botones contiene las siguientes herramientas de análisis:
• Entra en el modo para especificar condiciones de frontera.
• Abre el cuadro de diálogo para especificar la EDP a resolver.
• Inicializa la malla triangular.
• Refina la malla triangular.
• Resuelve la EDP.
• Abre el cuadro de diálogo para representar los resultados.
• Zoom on/off.
44
Utilización de la Toolbox PDE de Matlab
Como se mencionó esta toolbox nos permite resolver ecuaciones diferenciales en

derivadas parciales utilizando el método de elementos finitos. Con la utilización de
la interfaz gráfica denominada pdetool, es muy sencillo especificar el dominio 2-D
de nuestro problema, la triangulación del dominio, los coeficientes de la EDP y las
condiciones de frontera.
A continuación se resumen algunos elementos básicos para trabajar con el

toolbox.
• Para comenzar a trabajar con la interfaz gráfica para EDP de Matlab,

teclearemos pdetool en la línea de comandos.
• Para especificar el dominio de nuestro problema procederemos del
siguiente modo: Options, elegir Axes Limits para concretar el área (x,y) ,
(Por defecto, x=[-1.5:1.5], y=[-1:1] ); seleccionar Grid para situar un mallado
sobre el área que aparece y Snap para ajustar automáticamente la figura de
nuestro dominio al mallado. Seleccionar Gris Spacing si se desea el
espaciado del mallado.
• Los dominios se construyen en base a la composición (suma/resta) de una
serie de dominios elementales tales como rectángulos, polígonos y elipses,
cuya manipulación es muy sencilla. Se recomienda escoger un par de
ejemplos para comprobarlo. Para formar el área de interés en nuestro
problema, escribiremos en el espacio destinado a Set formula la “La
Fórmula” que describe nuestro dominio como suma y/o resta de las figuras
elementales. Si queremos guardar este dominio para utilizarlo a posterior,
seleccionaremos Export geometry description del menú de la opción Draw y
elegiremos OK en el cuadro mostrado. De este modo, “exportaremos” los
datos de la geometría bajo los nombres de las variables gd sf ns (salvo que
se cambie). De hecho, sin entrar en detalles, esta especificación de la
geometría debe ser procesada antes de poder utilizarla con diferentes
comandos Matlab. Es más sencillo esperar a tener especificadas las
condiciones de frontera y exportar simultáneamente éstas y la geometría
del problema, como comentaremos a continuación.
• Para establecer las condiciones de frontera, selecciona en Boundary la
opción Boundary mode o equivalente, dar clic al botón ∂Ω . Seleccionar uno
o varios segmentos y en el menú de Boundary elegir Specify Boundary
Conditions. Se pueden especificar condiciones de frontera tipo Dirichlet de
la forma h*u=r introduciendo los valores de los parámetros h y r. Si
elegimos condiciones de frontera tipo Neumann, deberemos especificar los
coeficientes g y q. La forma general de la condición de frontera aparece en
la parte superior del cuadro de diálogo.
• Para especificar la EDP a resolver, seleccionaremos PDE specification en
el menú de la opción PDE y escogeremos un problema prototipo. Introducir
los coeficientes pueden ser función f del lado derecho de la ecuación. Los
coeficientes pueden ser funciones de X, Y y para problemas no lineales
45
también pueden ser función de u, ux, uy. Utilizaremos entonces x, y, u, uy

en la expresión a introducir, dándonos cuenta de que Matlab interpreta
estas cantidades como vectores. Teniendo en cuenta esto, si queremos
introducir la función xy, teclearemos x.*y en lugar de x*y.
• Seleccionando Parameters en el menú de la opción Mesh estableceremos
los parámetros de la triangulación inicial de nuestro problema de elementos
finitos. Para generar la triangulación, pincharemos en el botón con el
triángulo o alternativamente, elegiremos Initialize Mesh en el menú de
Mesh. Para guardar la triangulación de forma que pueda ser utilizada en un
programa Matlab, seleccionaremos Export Mesh en el menú de Mesh y
daremos clic OK en el cuadro de diálogo. Esto permite exportar los vértices,
los bordes y el acomodo de los triángulos con los nombres p e t,
respectivamente. Démonos cuenta de que necesitaremos esta información
si queremos compara la solución exacta y la aproximada en los vértices o
calcular diversas normas del error que se comete.
• Para refinar la triangulación daremos clic al triángulo dividido o
seleccionaremos Refine Mesh en el menú de Mesh. El método de
refinamiento por defecto es regular. La solución numérica de la EDP puede
realizarse de manera adaptativa seleccionando primero Solve Parameters y
luego Adaptative Mode, en el menú de la opción Solve introduciremos el
máximo número de triángulos que deseamos el máximo número de pasos
de refinamiento.
• Para resolver la EDP y dibujar la solución, seleccionaremos Solve PDE en
el menú de Solve. Para exportar la solución al espacio de trabajo principal
de Matlab, elegiremos Expert Solution. Disponemos de varias opciones de
gráficos para representar la solución: seleccionaremos Paremeters del
menú de Plot.
Trabajando Desde La Línea de Comandos
Una vez resuelto nuestro problema de frontera utilizando pdetool y exportamos los
datos descritos con anterioridad, podemos procesar esta información desde la
línea de comandos de Matlab. En primer lugar, podría ser interesante guardar la
descomposición de la geometría de nuestro problema (g) en un fichero (probl.m)
que podamos utilizar más adelante. Una vez que hemos elegido en Boundary la
opción Export Descomponed Geometry, Boundary Conds, escribiremos en la línea
de comandos fid=wgeom(g,’prob1g’). Para hacer lo mismo con la información
sobre el contorno (b) que queremos en un fichero (prob1b.m) escribiremos
fid=wbound(b,’prob1b’).
De hecho, con estos dos ficheros es bastante fácil resolver un problema de

frontera sencillo sin necesidad de recurrir a la interfaz gráfica. Para ello,
establecemos en primer lugar una triangulación sobre la geometría definida en el
fichero prob1g.m escribiremos [p, e, t] = refimesh(‘prob1g’) el comando intimes
permite especificar determinados parámetros relativos a las característica de la
46
triangulación. Podemos visualizar la triangulación con pdemesh(‘prob1g’,p,e,t). Por

otro lado, los coeficientes de la EDP a resolver pueden ser especificados
directamente. Por ejemplo, la ecuación escalar genérica es la forma –div(c grad u)
+ a u = f.
Definición de geometrías complicadas
Las geometrías complicadas pueden generarse a partir de objetos sólidos básicos

(rectángulos/cuadrados, elipses/círculos y polígonos) que se solapen, parcial o
totalmente. La interfaz gráfica asigna automáticamente un nombre a cada objeto
sólido que se cree: R1, R2, … en el caso de los círculos y P1, P2, … para nombrar
los polígonos. Obviamente, estos nombres pueden ser modificados por el usuario
haciendo doble-click sobre los mismos, lo que abre el cuadro de diálogo también
nos posibilita el modificar otras características de la forma básica, como la
posición de su centro, dimensiones, etc.
Cuadro de diálogo con las propiedades de un cículo
Una vez que los objetos básicos han sido dibujados, la geometría final se crea
mediante la introducción, en la línea situada debajo de la barra de herramientas,
de una fórmula que emplee operaciones del álgebra de conjuntos, +, * y -. De
todos ellos, el operador de mayor precedencia es el operador diferencia, - ,
mientras que los operadores unión e intersección, + y * poseen igual prioridad. Sin
embargo, este orden de precedencia pueden controlarse mediante el uso de
paréntesis. El modelo geométrico final, Ω , es el conjunto de todos los puntos para
los cuales la fórmula introducida puede evaluarse como verdadera. El proceso
general puede entenderse más fácilmente a partir del siguiente ejemplo, la
creación de una placa con esquinas redondeadas.
Creando esquinas redondeadas
47
Inicia la interfaz gráfica y activa la propiedad de la rejilla “snap-to-grid” localizada

dentro del menú Options. Además cambia el espaciado a [-1.5:0.1:1.5] para el eje
x y [-1:0.1:1] para el eje y.
Selecciona el icono para crear rectángulos y usando el ratón dibuja uno de

anchura 2 y altura 1, comenzando en el punto (-1, 0.5). Para crear las esquinas
redondeadas añade círculos, uno en cada esquina. Los círculos deben tener radio
0.2 y centros a una distancia de 0.1 unidades de las fronteras izquierda/derecha y
superior/inferior del rectángulo [(-0.8, -0.3), (-0.8, 0.3) y (0.8, 0.3)]. Para dibujar
círculos en vez de elipses, recuerda usar el botón derecho del ratón o mantener la
tecla Ctrl pulsada mientras realizas el dibujo. Para finalizar, dibuja en cada una de
las esquinas un pequeño cuadrado de 0.2. Los objetos dibujados deberían
presentar el aspecto mostrado en la parte izquierda de la siguiente figura:
Izquierda: Objetos básicos a partir de los cuales será definida la geometría

definitiva.
Derecha: Modelo geométrico final.
Ahora debemos editar la fórmula que define la geometría. Para conseguir las
esquinas redondeadas, restemos los cuadrados pequeños del rectángulo y
sumemos a continuación los círculos. En forma de expresión de conjuntos como:
R1 – (SQ1 + SQ2 + SQ3 + SQ4) + C1 + C2 + C3 + C4
Presionando el botón podemos entrar en el modo Boundary y ver las

fronteras de la geometría final. Puede observarse que aún existen dentro de la
placa algunas de las fronteras provenientes de los subdominios originales. Si se
supone que la placa es homogénea, entonces podemos borrarlos. Para ello,
48
selecciona la opción “Remove All Subdomain Borders” del menú Boundary. Ahora
el modelo de la placa está completo.
Método sugerido de modelamiento
El flujo básico de acciones al emplear la interfase gráfica de la toolbox de

ecuaciones en derivadas es el indicado de izquierda a derecha por los menús y los
botones de la barra de herramientas, trabajando en este sentido a lo largo del
proceso de modelamiento, definición y resolución del problema. La siguiente
secuencia de acciones cubre todos los pasos de una sesión normal empleando
pdetool:
1. Usa la pdetool como herramienta de dibujo para realizar el dibujo de la
geometría 2-D en la que se quiere resolver la EDP, haciendo uso de los
objetos básicos y de las características de “fijar a rejilla”. Combina los
objetos sólidos mediante las fórmulas de álgebra de conjuntos para crear la
geometría definitiva.
2. Salva la geometría a un archivo de modelo ( .m) de manera que puedas
seguir empleándola en futuras sesiones de trabajo. Si salvas el archivo más
adelante a lo largo del proceso de resolución, el archivo del modelo también
incluirá ciertos comandos para recrear las condiciones de frontera, los
coeficientes de la EDP y la malla.
3. Pasa a especificar las condiciones de frontera presionando el botón si
las fronteras no son las correctas, puedes volver a editar la geometría
volviendo al modo de dibujo (Draw mode). Si durante la definición de la
geometría han quedado bordes de subdominios no deseados puedes
borrarlos mediante las opciones del menú Boundary (Remove Subdomain
Border “o “Remove All Subdomain Borders”). A continuación puedes fijar las
condiciones de cada una de las fronteras haciendo doble-click sobre cada
una de ellas.
4. Usa el botón para especificar la EDP a resolver. En el caso de que

los coeficientes de la EDP dependan del material, estos son introducidos
entrando en el modo PDE y haciendo doble-click en cada uno de los
subdominios.
5. Inicializa la malla triangular mediante el botón . Normalmente, los

parámetros por defecto del algoritmo de generación de la malla producen
buenos resultados, aunque en caso necesario pueden modificarse desde la
opción “Parameters” del menú Mesh.
6. Si es necesario, refina la malla mediante el botón . En cada

refinamiento, el número de triángulos aumenta en un factor cuatro. Ten
presente que cuando más fina sea la malla mayor será el tiempo requerido
para calcular la solución. Otra opción es reordenar la triangulación de la
49
malla para mejorar su calidad mediante la opción “Jiggle Mesh” del menú
Mesh.
7. Resuelve la EDP presionado el botón .

8. Visualiza las propiedades de la solución en las que estés interesado
mediante el botón . También puedes exportar la solución y/o malla al

espacio de trabajo principal de Matlab para un análisis en mayor detalle.
Ejemplo de aplicación a la ecuación de Poisson
Como ejemplo básico de una EDP elíptica resolvamos la Ecuación de Poisson en

un disco unitario con condiciones de frontera Dirichlet homogéneas. La
formulación del problema es la siguiente:
∇ 2u = 1 en Ω, u = 0 en δΩ
Donde Ω en el disco unitario y δΩ su frontera.
Tanto para éste como para el resto de ejemplos que desarrollaremos aquí
seleccionaremos el modo Generic Scalar de la lista desplegable de modos
disponible que se encuentra localizada a la derecha de la barra de herramientas. A
continuación se listan los pasos a realizar con la pdetool para resolver el problema
descrito.
1. Dibuja un círculo de radio unidad centrado en el origen usando las

herramientas de dibujo. Si dibujas una elipse en vez de un cículo o el
mismo no está perfectamente centrado, haz doble-click en el mismo y
cambia sus propiedades en el cuadro de diálogo.
50
2. Imposición las condiciones de frontera. Para ello primero pulsa el botón
y a continuación haz doble-click sobre las fronteras de la región.

Selecciona para todas ellas condiciones de frontera Dirichlet iguales a cero
(aunque éstas son las condiciones de frontera por defecto).
3. Define la EDP presionando el botón esto abre un cuadro de

diálogo donde pueden introducirse los coeficientes de la ecuación, c, a y f.
en caso tan sencillo todos son constantes: c=1, f=1 y a=0. En el caso de
que dependan de la posición, también pueden introducirse con la notación
habitual de Matlab en formas de productos de vectores.
4. Inicializa la malla haciendo uso del botón refínala con .
5. Resuelve la ecuación presionando emplea también si

quieres cambiar las propiedades del gráfico que la pdetool devuelve por
defecto, como por ejemplo el mapa de colores.
Ejemplos de EDP parabólica: La Ecuación de Calor.
Un problema parabólico típico es la Ecuación de Calor, que describe la difusión

de calor en un cuerpo:
∂u
d − ∇ 2u = 0
∂t
51
Como ejemplo estudiemos el comportamiento de calentar un bloque de metal que

posee una cavidad rectangular en el medio. El lado izquierdo del bloque se
calienta a 100º C. En la parte derecha del bloque, el calor fluye al aire circundante
en una tasa constante.
Solución a la ecuación de calor en la geometría estudiada.
Todas las otras fronteras del bloque se hayan aisladas térmicamente. Todo esto
nos conduce al siguiente conjunto de condiciones de contorno:
• u = 100 En el lado izquierdo (Condición Dirichlet)

∂u
• = −10 En el lado derecho (Condición Neumann)
∂n
∂u
• = 0 En todas las restantes fronteras (Condición Neumann)
∂n
Además, para la ecuación de calor necesitamos una condición inicial t0 . Para este
caso supondremos que la temperatura del bloque es de 0º C en el momento en el
que empezamos a aplicar el calor. Finalmente, para completar la formulación del
problema, especificaremos que el tiempo inicial de 0 y que queremos estudiar la
distribución de calor durante los cinco primeros segundos.
Resolvamos el problema utilizando la pdetool de Matlab. Los pasos a seguir para

resolver el problema son los siguientes:
1. En primer lugar, asegúrate de que se encuentra seleccionado el modo

Generic Scalar en la lista desplegable que se encuentra a la derecha de la
barra de herramientas.
52
2. Dibuja un rectangulo (R1) con las esquinas opuestas localizadas en (-0.5,

0.8) y (0.5, -0.8). Para crear la cavidad, dibuja otro rectangulo (R2) con
esquinas opuestas situadas en (-0.05, 0.4) y (0.05, -0.4). Realiza la
diferencia de ambas geometría, introduciendo en la línea de debajo de la
barra de herramientas la fórmula R1 y R2.
3. Impón las condiciones de frontera. Para ello primero pulsa el botón

y a continuación haz doble-click sobre las fronteras de la región. Para la
frontera izquierda introduce la condición Dirichlet u=100, para la derecha la
condición Neumann ∂u = −10 y el resto condiciones Neumann
∂n
homogéneas ∂u = 0.
∂n
4. Introduce los coeficientes que definen la EDP presionando el botón

para este caso, d=1, c=1, a=0 y f=0. Al ser un problema con evolución
temporal, también debemos introducir la condición inicial u0=0 y el rango de
tiempo en el que queremos resolver el problema [ 0:0.5:5 ]. Estos
parámetros son introducidos dentro del Menú Solve seleccionando la
opción “Parameters”.
5. Inicializando la malla haciendo uso del botón refínala con
6. Resuelve la ecuación presionando por defecto, la pdetool muestra
la distribución de calor para el final de la simulación. Presionando

puedes cambiar las propiedades del gráfico, por ejemplo, seleccionar un
colormap más apropiado, como el hot. También es interesante poder ver
cómo se mueve el calor a lo largo de la simulación. Para ello, activa la
casilla “Animation” dentro del cuadro de diálogo de las propiedades de
visualización. Por defecto, la animación se sujeta cinco veces, aunque este
número de veces también puede modificarse.
Ejemplos de problemas de valores propios: Modos de Vibración.
El encontrar los valores propios y las correspondientes funciones propias de una

membrana con forma de L es un problema de interés a todos los usuarios de
Matlab, dado que el logo del producto es precisamente la primera de estas
funciones propias. Formalmente, que intentamos es calcular la solución al
siguiente problema de valores propios.
− Δu = λ u
En una geometría con forma de membrana en L, con condiciones de frontera

Dirichlet homogéneas (u=0 en todas las fronteras).
53
1. Asegúrate de que se encuentra seleccionado el modo Generic Scalar en la

lista de los posibles modos de solución.
2. Dibuja la membrana en forma de L como un polígono con esquinas en (0,0),
(-1,0), (-1,-1), (1,-1), (1,1) y (0,1).
3. Define las condiciones de frontera. Pulsa el botón e introduce la
condición Dirichlet homogénea, u=0, para todas las fronteras de la
membrana (condición por defecto)
4. Introduce los coeficientes que definen la EDP presionando el botón
para este caso, d=1, c=1 y a=0. También debemos introducir el rango de
valores propios en el que estamos interesados. Para ello, selecciona dentro
del Menú Solve seleccionando la opción “Parameters”. En cuadro de
diálogo emergente introduce [ 0 100 ].
5. Inicializa la malla haciendo uso del botón refínala con .
Cuatro primeras funciones propias al problema de autovalores en la membrana

con forma de L, con λ1 = 9.68, λ2 = 15.22, λ3 = 19.79 y λ4 = 29.63
6. Resuelve la ecuación presionando presionando puedes cambiar

las propiedades de visualización y seleccionar de un menú desplegable las
distintas autofunciones que han sido calculadas como solución al problema.
Los Modos de Aplicación
54
La “PDE Toolbox” de Matlab puede aplicarse a una gran cantidad de problemas

habituales en ciencia e ingeniería. Ocho modos diferentes de aplicación se
encuentran implementados para facilitar el uso de pdetool:
• Mecánica de estructuras: Tensión Plana.

• Mecánica de estructuras: Deformación Plana.
• Electroestática
• Magnetoestática
• Electromagnetismo de potencia en corriente continua.
• Transferencia de calor
• Difusión
55

Notas Matlab Totales

Cargado por

Copyright:

Formatos disponibles

Notas Matlab Totales

Cargado por

Información del documento

Descripción original:

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Notas Matlab Totales

Cargado por

Copyright:

Formatos disponibles

Universidad Autónoma de Sinaloa

Facultad de Ciencias Químico Biológicas

Anabel Altamirano Aguilar

Las presentes notas surgieron para cubrir la necesidad de una referencia

Dr. Jesús Raúl Ortiz Del Castillo

Capítulo 1. Introducción al MATLAB

MATLAB es un sistema interactivo y lenguaje de programación para cómputo

En Macintosh o Windows, haga clic en el icono de MATLAB o de STUDENT

MATLAB usa tres ventanas de exhibición. La ventana de comandos sirve para

Hay varios comandos para despejar ventanas. El comando clc despeja la

También es importante saber como interrumpir o abortar un comando en

Si no conoce bien el significado de un comando teclee help y el nombre del

Nota: En pantalla aparecerá toda la información sobre el comando quit

version. Este comando reporta la versión de MATLAB que se está usando.

what “nombre de directorio”. Enlista los archivos del directorio en matlabpath.

clock. El comando clock exhibe la hora y la fecha con números como:

date. El comando date que proporciona la fecha.

path. El comando path imprime la ruta de búsqueda vigente de MATLAB. El

getenv. El comando getenv('MATLABPATH') muestra las rutas de MATLAB

2. Cálculos básicos en MATLAB

Si en la ventana de comandos (command window) aparece >> en la esquina

En el ejemplo anterior observe que cada línea es un comando y termina con un

Otra forma de imprimir los resultados de cada comando que se ejecute es

La indicación >> no aparecerá en la línea que siga a la marca de continuación.

Enunciado if. Enunciado condicional y siempre debe de terminar con un enunciado

Los operadores & y | se pueden utilizar agrupados.

El enunciado if se puede utilizar con else o elseif; por ejemplo,

disp. El comando disp exhibe un número, vector, matriz o cadena en la ventana de

imprimen 3.14159 en la ventana de comandos. Pruebe también disp “esta es una

2.1 Variables y nombres de variables

Cualquier variable puede adoptar valores reales, complejos y enteros. No es

2.2 Ciclos for/end y while/end

format. Por omisión, los números se exhiben con cinco dígitos.

Con format short se vuelve al formato corto.

break. El comando break termina la ejecución de un ciclo for o while. Hay

2.3 Lectura y escritura

a) Operación interactiva mediante teclado o el ratón

MATLAB puede aceptar datos de entrada a través del teclado mediante el

>>z = input (′indique su nombre: ′, ′s′)

\n sirve para que lo que se valla imprimiendo se imprima en reglones separados.

Es posible utilizar el enunciado fprintf para escribir salidas con formato de un

Por ejemplo fprintf('archivo_x', 'volumen =%12.5\n', vol) escribirá la salida en el

A = [3.5] B = [1.5 3.1]

Las variables de arreglo unidimensional tienen forma de fila y columna y están

En MATLAB hay varias formas de escribir un arreglo. A continuación se dan

La definición de un arreglo en columna es similar a la de un arreglo de fila

Cuando y y x tienen la misma longitud y la misma forma (fila o columna), los

La forma anterior es para ahorrarnos la siguiente forma larga de realizar esas

y para un vector columna se pone “ ;” en vez de “,”.

También se puede anteponer un elemento a un vector; por ejemplo

Para consultar el número de filas y columnas de nuestro arreglo se utiliza el

Con el comando size(nombre del arreglo) se obtiene el número de filas y

Una matriz en MATLAB, se puede definir especificando sus elementos. Por

Para hacer operaciones con arreglos bidimensionales se utilizan los mismos

4. Operaciones con matrices

En MATLAB se denota la transpuesta de A como A'.

El producto punto es un escalar calculado a partir de dos vectores del mismo

En la multiplicación de matrices, el valor que está en la posición ci,j del producto

Si A es una matriz, A.^2 es la operación que eleva al cuadrado cada elemento