Clase 02

Computación Científica
Rubén Agapito
Verano 2017
b
Índice general
Prefacio iii
1 Justificación del Curso 1

1.1 Soluciones en Computadora a Problemas Matemáticos . . . . . . . . . . . . . . . . . . . . 1
1.1.1 ¿Por qué debemos aprender a usar un software matemático? . . . . . . . . . . . . 1
1.1.2 Soluciones analíticas versus soluciones numéricas . . . . . . . . . . . . . . . . . . 6
2 Fundamentos de Programación en MATLAB y Visualización Científica 9

2.1 Elementos Esenciales de MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.1 Variables y Constantes en MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.2 Estructuras de Datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Fundamentos de Programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Entrada y Salida de Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.2 Problemas Propuestos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2.3 Bucles For y While . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.4 Problemas Propuestos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.2.5 Ploteo en Dos Dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Referencias 43
i
ii
Prefacio
Los contenidos de estos apuntes de clase están bajo constante supervisión. Si existe alguna sugeren-
cia o error, por favor mandar un correo electrónico al autor: [email protected]
Rubén Agapito
San Miguel, febrero de 2017
Lima, Perú.
iii
iv Prefacio
Justificación del Curso
1
Cuando resolvemos un problema en Ciencias e Ingeniería un modelo matemático es establecido pri-
mero, luego, una solución a este modelo puede ser obtenido. Métodos numéricos y analíticos serán
presentados en este libro para resolver esta clase de problemas.
Sección 1.1
Soluciones en Computadora a Problemas Matemáticos
1.1.1
¿Por qué debemos aprender a usar un software matemático?
Sabemos que una derivación manual de soluciones a problemas matemáticos es una habilidad útil
cuando los problemas no son tan complicados. Sin embargo, para una gran variedad de problemas,
una solución manual es a veces tan laboriosa o incluso imposible. Por lo tanto, computadoras deben
ser empleadas para resolver estos problemas. Existen básicamente dos maneras de resolver estos
problemas en una computadora. La primera es implementar algoritmos numéricos existentes usando
lenguajes de programación como C o Fortran. La otra es usar lenguajes matemáticos específicos con
una buena reputación, como Matlab, Mathematica y Maple.
Algoritmos numéricos únicamente pueden ser usados para manipular problemas que involucren
data numérica, mientras que problemas simbólicos como resolver x3 + ax + c = d, donde a, c, d no
son números sino variables simbólicas, no pueden resolverse con algoritmos numéricos. Un lenguaje
matemático con capacidades de cálculo simbólico debe ser usado en su lugar.
Cuando enfrentemos un problema, soluciones analíticas son exploradas primero, y si estas son
difíciles de conseguir, soluciones numéricas son obtenidas.
Presentamos ahora algunos ejemplos que ilustren la necesidad de usar un lenguaje matemático
como Matlab.
Ejemplo 1
sen x d4 f (x)
Si una función f (x) está dada por f (x) = , ¿cómo obtenemos manualmente ?
x2 + 4x + 3 dx4
Solución. Podemos derivar sucesivamente cuatro veces a la función f como se ha aprendido en
los cursos de cálculo. Sin embargo, el procedimiento es más conveniente realizarlo en una compu-
tadora. La implementación en un script de Matlab es la siguiente:
1
2 Justificación del Curso
syms x
f = sin ( x ) /( x ^2+4* x +3) ;
y = diff (f ,x ,4) ;
y = simplify ( y ) ;
pretty ( y )
Una vez compilado obtenemos
2 3 4 5 6
(1581 sin(x) − 2448 cos(x) − 4080 x cos(x) − 296 x cos(x) + 1040 x cos(x) + 552 x cos(x) + 112 x
7 2 3 4 5 6
cos(x) + 8 x cos(x) − 3120 x sin(x) − 3120 x sin(x) − 1094 x sin(x) − 32 x sin(x) + 72 x
7 8 2 5
sin(x) + 16 x sin(x) + x sin(x) − 5640 x cos(x) + 192 x sin(x))/(x + 4 x + 3)
en el Command Window. Es fácil ver que esta expresión significa que
d4 f ( x) ¡
= 1581 sin( x) − 2448 cos( x) − 4080 x2 cos( x) − 296 x3 cos( x) + 1040 x4 cos( x) + 552 x5 cos( x)
d x4
+ 112 x6 cos( x) + 8 x7 cos( x) − 3120 x2 sin( x) − 3120 x3 sin( x) − 1094 x4 sin( x) − 32 x5 sin( x)
¢5
+ 72 x6 sin( x) + 16 x7 sin( x) + x8 sin( x) − 5640 x cos( x) + 192 x sin( x) / x2 + 4 x + 3 .
¢¡
Es obvio que una derivación manual podría ser un trabajo tedioso y laborioso. Es más, un leve
error en la manipulación podría dar resultados desastrosos en la respuesta final. Si usamos un
lenguaje matemático el trabajo tedioso y posiblemente no confiable puede ser evitado. Por ejemplo,
en Matlab podemos calcular exactamente d50 f (x) / dx50 , y puede ser obtenido en aproximadamen-
te ¡un segundo!. Cabe señalar que el cálculo de esta derivada no es visualizado completamente
en Matlab, debido a que el output excede la cantidad máxima de 25000 caracteres. En caso sea
necesario obtener la respuesta completa añadimos las siguientes líneas a nuestro script:
diary ( ' mivar . txt ')

y
diary off
Esto creará el archivo mivar.txt en donde se almacenará completamente el contenido de la varia-

ble y.
Ejemplo 2
En algunas aplicaciones el cálculo de las raíces de un polinomio es requerido. El teorema de Abel-

Ruffini afirma que no hay una fórmula general para obtener las raíces de polinomios de grado
cinco o mayor que use solo radicales y/o operaciones aritméticas. Se sabe que el problema puede
ser resuelto numéricamente usando algoritmos conocidos.a Resuelva la ecuación polinomial
134 4 135 3 1215 2 729 729

x6 + 9x5 + x + x + x + x+ = 0.
4 2 16 16 64
Solución. En la primera línea del siguiente script se resuelve numéricamente el problema, sin
1.1. SOLUCIONES EN COMPUTADORA A PROBLEMAS MATEMÁTICOS 3
embargo, en la segunda línea intentamos resolver el problema simbólicamente:
p = [1 9 135/4 135/2 1215/16 729/16 729/64]; roots ( p ) % sol . numerica

p1 = poly2sym ( p ) ; solve ( p1 ) % solucion analitica
La respuesta numérica se muestra a continuación:
ans =
-1.5052 + 0.0030i
-1.5052 - 0.0030i
-1.5000 + 0.0060i
-1.5000 - 0.0060i
-1.4948 + 0.0030i
-1.4948 - 0.0030i
la cual no es confiable, ya que se tuvo éxito en la respuesta simbólica, la cual indica que la solución
es real con multiplicidad algebraica cinco:
ans =
-3/2
-3/2
-3/2
-3/2
-3/2
-3/2
a Ver https://en.wikipedia.org/wiki/Root-finding_algorithm#Finding_all_roots_at_once
Ejemplo 3
¿Recuerda Ud. cómo encontrar el determinante de una matriz n × n?

Solución. En un curso de Algebra Lineal, se enseña que una manera de calcular el determinante
de una matriz n × n es por el método de expansión de Laplacea (o expansión en cofactores), el
cual implica el cálculo de los determinantes de n matrices de tamaño (n − 1) × (n − 1). A su vez,
el determinante de cada matriz (n − 1) × (n − 1) se obtiene al calcular los determinantes de n − 1
matrices de tamaño (n − 2) × (n − 2). Eventualmente se termina con determinantes de matrices de
tamaño 1 × 1, los cuales son los mismos escalares. Así, puede concluirse que la solución analítica
del determinante de cualquier matriz existe.
Se puede probar que el número de operaciones aritméticas que se requiere en el método de ex-
pansión de Laplace es (n − 1)(n − 1)! + n. Por ejemplo, cuando n = 25, el número de operaciones punto
flotante (flops, en inglés) es 9.679 × 1027 , lo cual equivale a 5580 años de cálculo en mainframes
capaces de realizar 55 000 billones (5.5 × 1016 ) de flops por segundo (esta velocidad corresponde al
mainframe más rápido del mundo del año 2014). Por lo tanto, el método de expansión en cofactores,
aunque elegante e instructivo, no es prácticamente útil. En aplicaciones reales, el determinante de
matrices de tamaño aún más grande es necesitado (n À 25), pero es obvio que no es posible aplicar
el método de cofactores.
En cursos de análisis numérico varios algoritmos han sido diseñados para calcular el determi-
nante de una matriz. Sin embargo, debido a errores de redondeo cometidos por la computadora,
estos algoritmos pueden resultar no confiables cuando la matriz es muy cercana a ser singular (o
no invertible). Por ejemplo, consideremos la matriz de Hilbert dada por

 
1 1/2 1/3 ··· 1/n
 1/2 1/3 1/4 · · · 1/(n + 1) 
 
H =  .. .. .. .. .. .
 . . . . .


1/n 1/(n + 1) 1/(n + 2) ··· 1/(2n − 1)
Para n = 25, un valor erróneo, det(H) = 0, es obtenido, inclusive si precisión doble es usado. Esto
es lo que conseguimos con el siguiente código:
>> n = 80; H = hilb ( n ) ; det ( H )
Por otro lado, si cálculo simbólico es usado,
>> H = sym ( hilb ( n ) ) ; d = det ( H )
obtenemos (en aproximadamente 4 segundos)
d =
1/9903010146699347787886767841019251068...00000000
Aquí se tienen 3790 dígitos en el denominador, y no se han mostrado todos por motivos de espacio.
Como la variable d es simbólica, calculamos su valor numérico con el comando vpa, usando la
sintaxis vpa(d), obtenemos
ans =
1.0097939769690106521663398695661e-3790
lo cual es más confiable, y nos indica que el determinante no es exactamente cero.

a Ver https://en.wikipedia.org/wiki/Laplace_expansion
Ejemplo 4
Considere la ecuación no lineal de Van der Pol
y00 + µ(y2 − 1)y0 + y = 0, y(0) = −1, y0 (0) = 1.
Cuando µ = 1000, muestre cómo resolver numéricamente la ecuación en el intervalo de tiempo

[0, 3000].
Solución. Algoritmos numéricos estándares para resolver ecuaciones diferenciales ordinarias, ta-
les como el de Runge-Kutta clásico de orden 4, pueden causar problemas numéricos. Algoritmos
especializados para ecuaciones diferenciales rígidas (stiff ordinary differential equations) deben
ser usados. Las siguientes sentencias en Matlab son suficientes para resolver esta ecuación:
>> mu = 1000;
>> f = @ (t , x ) [ x (2) ; - mu *( x (1) ^2 -1) * x (2) -x (1) ]; % describir EDO como sistema
>> [t , x ] = ode15s (f ,[0 ,3000] ,[ -1;1]) ; % obtencion de la solucion numerica
>> plot (t , x (: ,1) ) , grid on % resolucion grafica
Obtenemos la siguiente gráfica:

2.5
1.5
0.5
-0.5
-1
-1.5
-2
0 500 1000 1500 2000 2500 3000
Ejemplo 5
Resuelva el siguiente problema de programación lineal
mı́n − 2x1 − x2 − 4x3 − 3x4 − x5 ,


 2x2 + x3 + 4x4 + 2x5 ≤ 54,

sujeto a 3x1 + 4x2 + 5x3 − x4 − x5 ≤ 62,

 x , x ≥ 0, x ≥ 3.32, x ≥ 0.678, x ≥ 2.57
1 2 3 4 5
Solución. Como la minimización es un problema de optimización lineal restringida, el método

analítico sin restricciones, esto es, poner las derivadas parciales de la función objetivo con res-
pecto a cada variable de decisión x i iguales a cero, no puede ser usado. Con las herramientas de
programación lineal de Matlab, las siguientes sentencias pueden ser usadas:
>> P . f = [ -2 -1 -4 -3 -1]; P. Aineq = [0 2 1 4 2; 3 4 5 -1 -1];

>> P . Bineq = [54; 62]; P. lb =[0; 0; 3.32; 0.678; 2.57];
>> P . options = optimoptions ( ' linprog ' , ' Algorithm ' ,' dual - simplex ') ;
>> x = linprog (P )
y obtenemos
Optimal solution found.

x =
19.7850
0
3.3200
11.3850
2.5700
De aquí la solución numérica encontrada se traduce a x1 = 19.7850, x2 = 0, x3 = 3.3200, x4 =

11.3850 y x5 = 2.5700.
Si las variables de decisión están restringidas a ser números enteros, debemos usar programa-
ción entera. En M ATLAB, la solución a este nuevo problema es fácilmente encontrada al compilar
las siguientes sentencias:
>> P . solver = ' intlinprog '; P . options = optimoptions ( ' intlinprog ') ;
>> P . intcon = 1:5; x = intlinprog ( P )
lo cual nos da
x =
19.0000
0
4.0000
10.0000
5.0000
1.1.2 Soluciones analíticas versus soluciones numéricas

El desarrollo de la ciencia moderna e ingeniería depende grandemente de la matemática. Los in-
tereses de investigación de matemáticos puros son diferentes de otros científicos e ingenieros. Los
matemáticos usualmente no están muy interesados en encontrar soluciones analíticas o fórmulas
cerradas a problemas matemáticos, más bien se interesan en probar la existencia y unicidad de
las soluciones, y no en hallar las soluciones explícitas en sí. Por el contrario, los ingenieros están
más interesados en encontrar las soluciones exactas o aproximadas a problemas dados, y usualmen-
te no prestan demasiada atención a cómo fueron obtenidos los resultados, siempre que estos sean
confiables y significativos. La manera más usual de encontrar soluciones aproximadas es el cálculo
numérico.
Es bastante común encontrar que soluciones analíticas no existen en varias ramas de la mate-
mática. Por ejemplo, es bien conocido que la integral definida
Z a
2 2
p e− x dx
π 0
no tiene valor exacto. Para resolver este problema, los matemáticos introducen una función espe-
cial, la función error, denotada por erf(a), y no les importa el valor numérico para un a dado. Para
encontrar un valor aproximado, los ingenieros usan el cálculo numérico.
Otro ejemplo es el número irracional π, el cual no tiene una fórmula cerrada, pero para propósitos
prácticos lo aproximamos con 4 o 6 cifras.
Estas discusiones sugieren que soluciones numéricas aproximadas abundan en la práctica. El
conocimiento matemático teórico de uno podría no corresponder a la habilidad de poder resolver pro-
blemas matemáticos prácticos. En ingeniería y ciencias aplicadas de hoy, uno necesita usualmente
resolver problemas matemáticos que se presenten, en corto tiempo, sin entender completamente las
técnicas numéricas involucradas, inclusive durante el proceso de solución. Por lo tanto, en la actua-
lidad es una tendencia que uno se enfoque más en cómo formular el problema de manera que tenga
una forma conveniente para ser resuelto por computadora, y en la interpretación de los resultados
generados por computadora.
Técnicas numéricas ya han sido usadas en varias áreas de la Ciencia e Ingeniería. Por ejem-
plo, en mecánica, métodos de elementos finitos (FEM, de finite element methods) se han utilizado
para resolver ecuación diferenciales parciales. En ingeniería de control y aeroespacial, álgebra li-
neal numérica y soluciones numéricas de ecuaciones diferenciales ordinarias han sido utilizados
exitosamente por décadas. Para simulación de experimentos, soluciones numéricas de ecuaciones en

diferencia y ecuaciones diferenciales constituyen el núcleo del problema. En desarrollos de alta tec-
nología, procesamiento de señales digitales basado en la transformada rápida de Fourier (FFT, de
fast Fourier transform), es una rutina diaria. No hay duda que si uno domina una o más herramien-
tas computacionales prácticas, un mejoramiento significativo en la habilidad de resolver problemas
matemáticos es esperado.
Fundamentos de Programación en MATLAB y Visualización
2
Científica
En este libro nos concentraremos en presentar M ATLAB con sus aplicaciones en la resolución de
problemas de matemática aplicada.
Considerado como lenguaje de programación, M ATLAB tiene las siguientes ventajas:
1. Claridad y alta eficiencia. M ATLAB es un lenguaje altamente integrado. Unas pocas líneas de
código en M ATLAB equivalen a cientos de líneas de código fuente de otros lenguajes. Por lo tanto,
un programa en M ATLAB es más confiable y fácil de mantener.
2. Computación científica. El elemento básico en M ATLAB es una matriz compleja de precisión do-
ble. Manipulaciones matriciales pueden realizarse directamente. Funciones predefinidas en M AT-
LAB que implementan herramientas de cálculo numérico pueden ser usadas directamente. Tam-
bién, computación simbólica provista por la caja de herramientas Symbolic Math Toolbox puede
utilizarse para obtención de fórmulas, con el fin de resolver analíticamente un problema.
3. Funciones de graficación. M ATLAB puede ser utilizado para visualizar data experimental de
una manera muy fácil. Funciones implícitas también pueden ser graficadas con facilidad. Más
aún, creación de interface gráfica de usuario (GUI, de graphical user interface) y programación
orientada a objetos, también son posibles de llevar a cabo en M ATLAB.
4. Diversos toolboxes y blocksets. Existe una enorme cantidad de toolboxes (librerías creadas
para un fin específico) y Simulink blocksets (bloques gráficos que permiten entrada y/o salida
de data, que permiten realizar simulaciones) desarrollados por programadores experimentados e
investigadores.
5. Poderoso sistema de simulación. La poderosa técnica de modelación basada en diagrama de

bloques por Simulink, puede ser utilizada para analizar sistemas con casi cualquier grado de
complejidad. En particular, bajo Simulink, bloques de control, de electrónica y de mecánica pueden
ser modelados de manera conjunta bajo el mismo marco de trabajo, lo cual actualmente no es
posible en otros lenguajes matemáticos computacionales.
9
10 Fundamentos de Programación en MATLAB y Visualización Científica
Sección 2.1
Elementos Esenciales de MATLAB
2.1.1 Variables y Constantes en MATLAB

El nombre de una variable en M ATLAB debe comenzar con una letra (mayúscula o minúscula) y
puede ser seguida con letras, números o el carácter de subrayado, con un máximo de 63 caracteres.
Por ejemplo, MiVar01, Mi_Var01 y mi_var_01 son nombres de variable válidos, mientras que 1raVar y
mi.var son inválidos, dado que el primero empieza con un número y el segundo contiene un punto, el
cual es un carácter inválido para una variable simple.1 M ATLAB es sensible a la capitalización, esto
es, variables como abc, ABC y aBc son todas diferentes.
Es importante aprender de memoria algunas palabras reservadas (keywords, en inglés) de M AT-
LAB , ya que usarlas como posibles nombres de variables (definidas por el usuario), nos traería pro-
blemas. Tenemos constantes con valores predefinidos.
(1) eps. Es el error de tolerancia para operaciones punto flotante. El valor aproximado es 2.2204 ×
10−16 . También puede ser considerado, algunas veces, como la distancia entre dos números reales
punto flotante contigüos. Obsérvese que si el valor absoluto de un número punto flotante es menor
que eps, entonces en máquina será tratado como cero.
p
(2) i y j. Ambos representan −1. Si alguno de ellos fuese sobrescrito, su valor puede ser restaurado
con el comando i=sqrt(-1), o el comando clear i
(3) inf. Es la representación en M ATLAB de +∞. También puede ser escrito como Inf. Similarmente,
−∞ puede ser tipeado como -inf. Cuando 0 es usado en el denominador, el valor inf puede ser
generado.
(4) NaN. Significa not a number, lo cual es equivalente a una expresión indeterminada, como 0/0, ∞/∞,
∞ − ∞, 1∞ , 0 · ∞. Debemos tener cuidado con las otras dos expresiones indeterminadas, ya que en
M ATLAB se tiene 00 = 1 e ∞0 = 1.
(5) pi. Es la representación en máquina, en precisión doble, del número irracional π. Si deseamos
ver con cuántas cifras está almacenado en memoria este número, digitamos format long y luego
pi.
(6) realmax. Es el número real positivo punto flotante más grande que puede ser almacenado en
máquina. Su valor aproximado es 1.7977 × 10308 .
(7) realmin. Es el número real positivo punto flotante más pequeño que puede ser almacenado en
máquina. Su valor aproximado es 2.2251 × 10−308 .
y tenemos algunas palabras reservadas usuales muy útiles en programación:
for end while if else elseif fprintf case disp

save load max min sum print doc help format
mod sin cos tan exp log abs sqrt plot
grid doc help input size title ones zeros length
1 Más tarde veremos que un punto (como también el uso de paréntesis) nos sirve para crear variables de tipo estructu-
rado, ver https://www.mathworks.com/help/matlab/matlab_prog/create-a-structure-array.html

2.1. ELEMENTOS ESENCIALES DE MATLAB 11
2.1.2 Estructuras de Datos

I. Tipo de dato Numérico
Para asegurar alta precisión en los cálculos, M ATLAB usa data de precisión doble punto flotante, la
cual es de 8 bytes (64 bits). De acuerdo al estándar IEEE, está formado por 11 bits para la parte
exponencial, 52 bits para la parte fraccionaria y 1 bit para el signo. Esto llega a representar data con
un rango de ±1.7 × 10308 , con alrededor de 15 a 17 cifras significativas. La función en M ATLAB para
definir este tipo de data es double().
En otras aplicaciones especiales, como en procesamiento de imágenes, un entero sin signo de 8
bits puede ser usado, cuya función es uint8(), el cual permite un rango numérico de números enteros
desde 0 a 255. Así, un significativo espacio en memoria puede ser ahorrado. Otros tipos de data son
int8(), int16(), int32(), uint16(), y uint32().
II. Tipo de dato Simbólico

Variables simbólicas también pueden ser definidas en M ATLAB. Estas pueden ser usadas para de-
rivación de fórmulas y para encontrar soluciones analíticas de problemas matemáticos. Antes de
encontrar una solución analítica, las variables que van a ser usadas deben ser declaradas como
simbólicas, por medio del comando syms. Su sintaxis de uso es syms listaVar listaProp, donde
listaVar es la lista de variables a ser declarada, separadas con espacios en blanco. Opcionalmente,
el tipo de las propiedades de las variables puede ser asignado por listaProp, el cual puede ser real
o positivo. Por ejemplo, si uno quiere que a y b sean variables simbólicas, la sentencia syms a←-
b puede ser usada. Además, la sentencia syms a real puede ser usada para indicar que a es una
variable real.
Suposiciones adicionales pueden ser hechas para variables simbólicas específicas, para lo cual se
usa la función assume(). Por ejemplo
>> syms t x ; assume ( t <= 4) ; assume ( x ~= 2) ;
declara que t ≤ 4 y x 6= 2. La función assumeAlso() también puede ser usada para poner más especifi-
caciones sobre variables simbólicas. Por ejemplo, si x es real y −1 ≤ x < 5, la declaración en M ATLAB
sería
>> syms x real ; assume (x >= -1) ; assumeAlso (x <5)
Alternativamente, podríamos haber usado una forma más compacta
>> syms x real ; assume (x >= -1 & x <5)
En caso hayamos definido una variable simbólica, y luego queremos saber el tipo de restricciones
impuestas sobre ella, basta usar el comando assumptions().
La función de precisión aritmética variable (variable precision arithmetic) vpa() es usada para
exhibir el valor numérico de una variable simbólica con cualquier precisión. La sintaxis de la función
es vpa(A,n) o vpa(A), donde A es la variable simbólica y n el número de dígitos requerido, con el valor
por defecto de 32 dígitos decimales.
Ejemplo 1
Exhiba los primeros 300 dígitos decimales de π.

Solución. Basta digitar lo siguiente
>> vpa ( pi ,300)
lo cual nos da
ans =
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482←-
5342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559←-
6446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104←-
5432664821339360726024914127
III. Otros tipos de datos

Aparte de los tipos de data más usados en M ATLAB, también pueden emplearse los siguientes tipos
cuando sea necesario:
1) Strings. Variables tipo string son usadas para almacenar mensajes. El string debe estar entre
diéresis. Por ejemplo
>> s = ' Hola mundo '

>> size ( s )
nos da una variable tipo string almacenada en s, de tamaño 1 × 10.
2) Arreglos multidimensionales. Arreglos tridimensionales y de dimensiones mayores son la ex-

tensión del concepto de matriz. Arreglos unidimensionales son llamados vectores y arreglos bidi-
mensionales son llamados matrices. Corchetes son usados para digitar arreglos en el Command
Window. Por ejemplo, un vector fila puede ser asignado como sigue:
>> a = [10 20 30 25 13]
a =
10 20 30 25 13
Un vector columna puede ser ingresado de varias maneras
>> b = [2; 4; 6; 8; 10]
o también
>> b = [2
4
6
8
10]
o tomando la transpuesta a un vector fila, con el operador ',
>> b = [2 4 6 8 10] '
Una matriz 3 × 3 puede ser asignada como sigue
>> A = [1 2 3; 4 5 6; 7 8 9]
Además, podríamos construir la misma matriz por concatenación (esto es, uniendo) los vectores
que representan cada columna
>> A = [[1 4 7] ' [2 5 8] ' [3 6 9] ']
En este punto de la sesión, una lista de todas las variables actualmente en memoria se consigue
con el comando who:
>> who
Your variables are:
A a ans b
o si deseamos más detalle, usamos el comando whos:
>> whos
Name Size Bytes Class Attributes
A 3x3 72 double
a 1x5 40 double
ans 1x1 8 double
b 5x1 40 double
Para accesar a un elemento individual de un arreglo, usamos paréntesis. Por ejemplo, el cuarto
elemento del vector columna b puede ser exhibido como
>> b(4)
ans =
8
Para una matriz, A(m,n) selecciona el elemento de la fila m y columna n. Por ejemplo,
>> A(3,1)
ans =
7
Hay varias funciones predefinidas que pueden ser usadas para crear matrices. Por ejemplo, las
funciones ones y zeros crean arreglos rellenados con unos y ceros, respectivamente. Ambos tienen
dos argumentos, el primero es para el número de filas y el segundo para el número de columnas.
Por ejemplo, creemos una matriz 2 × 3 de ceros:
>> C = zeros(2,3)
C =
0 0 0
0 0 0
Similarmente, la función ones puede ser usada para crear un vector fila de unos:
>> u = ones(1,4)
u =
1 1 1 1
3) Arreglos Celulares. A diferencia de una matriz, los elementos (células o celdas) de un arreglo
celular no tienen porqué ser solo del tipo numérico. Por ejemplo, en M ATLAB no es posible formar
la siguiente matriz:
>> x = [ ' hola ' 55; ' mundo ' 2]
dado que elementos del arreglo x tienen diferentes dimensiones. La entrada x(1,1) es de tipo
string y puede ser considerado como un arreglo de tamaño 1 × 4, mientras que la entrada x(1,2)
es de tipo numérico y es de tamaño 1 × 1. De allí se suscita el problema. Pero si ahora construimos
el siguiente ejemplo un poco más elaborado (nótese el uso de llaves en vez de paréntesis):
>> A = [1 3; 7.1 9]; x = { ' hola ' 55; ' mundo ' A }
obtenemos
x = 2x2 cell array
'hola' [ 55]
'mundo' [2x2 double]
Describimos el arreglo celular entrada por entrada:
• Celda (1,1) es un arreglo string de 4 caracteres.

• Celda (1,2) es un arreglo doble 1 × 1.
• Celda (2,1) es un arreglo string de 5 caracteres.
• Celda (2,2) es un arreglo doble 2 × 2.
¿Cómo invocamos las celdas de este arreglo celular? Veamos algunos ejemplos:
• x{1,1} se refiere al string 'hola'

• x{2,2} se refiere al arreglo doble [1 3; 7.1 9]
• ¿A qué se refiere x{2,1}(3)?
• ¿A qué se refiere x{2,1}(1,2)?
4) Variables Estructuradas. Una estructura es una colección de data que representa una idea
simple u objeto. Dentro de una estructura tenemos una lista de campos (fields), y cada campo
es una subcolección de data. Entendamos este tipo de variable con el siguiente ejemplo. Primero
creamos una estructura de arreglo 1 × 1 llamada paciente en un script:
paciente . nombre = ' Sancho Panza ';

paciente . cuenta = 250.00;
paciente . test = [80 25 83; 170 150 189.5; 300 311 305];
Ahora digitamos en el Command Window paciente y obtenemos:
paciente =
struct with fields:
nombre: 'Sancho Panza'

cuenta: 250
test: [3x3 double]
Así, paciente es un arreglo que contiene una estructura con tres campos. Para expandier la es-
tructura de arreglo, añadimos índices luego del nombre de la estructura:
paciente (2) . nombre = ' Zoila Vaca ';

paciente (2) . cuenta = 30.50;
paciente (2) . test = [65 70 71; 110 114 109; 120 200 201];
La estructura de arreglo paciente ahora es de tamaño 1 × 2. Observe que una vez que una estruc-
tura contenga más de un simple elemento, M ATLAB ya no exhibe los contenidos de cada campo
cuando invoquemos el nombre del arreglo, más bien muestra un resumen del tipo de información
que la estructura contiene
>> paciente
paciente =
1x2 struct array with fields:
nombre
cuenta
test
Podemos también usar la función fieldnames() para obtener esta información. La diferencia es
que obtendríamos un arreglo celular de strings con los nombres de todos los campos:
>> fieldnames(paciente)
ans =
3x1 cell array
'nombre'
'cuenta'
'test'
Si expandimos la estructura, M ATLAB rellenará cualquier campo no especificado con una matriz
vacía, de manera que todas las estructuras en el arreglo tengan el mismo número de campos.
Por ejemplo, si digitamos paciente(3).nombre = 'Juan Pistolas', expandimos la estructura
paciente a un arreglo 1 × 3, y los campos paciente(3).cuenta y paciente(3).test contendrán
arreglos vacíos.
Tomar en cuenta que el tamaño de los campos no tienen que ser igual para cada elemento en un
arreglo. Por ejemplo, en la estructura paciente el campo nombre tiene diferentes longitudes, y el
campo test puede contener matrices de diferentes tamaños.
5) Clases y Objetos. M ATLAB permite el uso de clases en la programación, esto es, soporta progra-
mación orientada a objetos (OOP, Object Oriented Programming).2
Sección 2.2
Fundamentos de Programación
Un algoritmo3 es un conjunto de reglas que deben seguirse ordenadamente para resolver un proble-
ma dado. Cuando escribimos un algoritmo en papel, comúnmente usamos palabras del habla común
mezcladas con lenguaje técnico. El resultado de este proceso es llamado seudocódigo,4 y por la
estrecha relación entre ambos, usaremos algoritmo y seudocódigo como sinónimos.
Un lenguaje de programación (como C, C++, Java, Python, etc.) es un lenguaje formal que
nos permite controlar el comportamiento de una máquina, en especial, una computadora. Debemos
tener en cuenta que un seudocódigo es escrito independientemente del lenguaje de programación
que uno escoja. Cuando implementamos un seudocódigo en la computadora, vía algún lenguaje de
programación, obtenemos un programa. Un programa es así un conjunto de instrucciones que una
computadora puede interpretar y ejecutar. A diferencia del seudocódigo, un programa debe elaborar-
se cumpliendo ciertas reglas de sintaxis y de formato (determinadas por el lenguaje de programación
que se use), caso contrario, el programa no funcionará (o compilará) correctamente.
Ilustramos los conceptos mencionados en el siguiente ejemplo.
Deseamos encontrar una solución de la ecuación f (x) = 0, donde se sabe que f es continua en el
intervalo [a, b] y f (a) · f (b) < 0.5 De los varios algoritmos que existen presentamos el siguiente:
Algoritmo de Bisección
INPUT extremos a, b; tolerancia TOL; maximo numero de iteraciones N0 .
OUTPUT solucion aproximada r o mensaje de no exito.
Paso 1: Inicializar i = 1;
Poner F A = f (a)
Paso 2: Mientras que i ≤ N0 hacer Pasos 3-6.
Paso 3: Poner r = a + ( b − a)/2;
FR = f ( r )
Paso 4: Si FR = 0 o (b − a)/2 < TOL entonces
OUTPUT (r); (Procedimiento completado exitosamente)
PARAR.
Paso 5: Poner i = i + 1
Paso 6: Si F A · FR > 0 entonces (calculo de nuevo intervalo [a i , b i ])
poner a = r ;
F A = FR ;
sino
poner b = r (F A conserva su valor)
Paso 7: OUTPUT('Metodo no converge despues de N0 iteraciones');
(El procedimiento no tuvo exito).
PARAR.
2 Verhttps://www.cs.ubc.ca/~murphyk/Software/matlabTutorial/html/objectOriented.html
3 Léase http://es.wikipedia.org/wiki/Algoritmo para más información.
4 Léase http://es.wikipedia.org/wiki/Algoritmo#Pseudoc.C3.B3digo.
5 El requerimiento f (a) · f (b) < 0 implica que f evaluada en los extremos del intervalo [a, b] toma valores con signos
opuestos y como f es continua, la gráfica de la función debe pasar necesariamente por el eje X .
2.2. FUNDAMENTOS DE PROGRAMACIÓN 17
Este algoritmo es implementado en el siguiente programa implementado en M ATLAB.6
Programa 2.1: Función biseccion.m

1 function r = biseccion (f ,a ,b , tol , nmax )
2 % Esta funci ó n calcula el valor aproximado de la
3 % ra í z de la ecuaci ó n f ( x ) =0 , localizada en [a , b ]
4 % Input :
5 % f : funci ón , lado izquierdo de la ecuaci ó n
6 % a , b : extremos del intervalo [a , b ]
7 % tol : tolerancia
8 % nmax : n ú mero m á ximo de iteraciones
9 % Output :
10 % r : soluci ó n ( ra í z ) aproximada de f ( x ) =0
11 i =1;
12 FA = f ( a ) ;
13 while i <= nmax
14 r = a + (b - a ) /2;
15 FR = f ( r ) ;
16 if FR ==0 || (b - a ) /2 < tol
17 fprintf ( 'El valor aprox . de la ra í z es %6.3 f ' ,r ) ;
18 return ;
19 end
20 i = i +1;
21 if FA * FR >0
22 a = r;
23 FA = FR ;
24 else
25 b = r;
26 end
27 end
28 disp ( ' El m é todo no converge , n ú mero m á ximo de \ n ') ;
29 disp ( ' iteraciones alcanzado .\ n ') ;
En estas notas usaremos el lenguaje de programación desarrollado en el software M ATLAB (del

acrónimo Matrix laboratory, en inglés). Este lenguaje es un híbrido de los lenguajes conocidos como
C, C++ y Fortran. La conveniencia de usar M ATLAB es que tiene “casi todo en un solo lugar”, esto
es, librerías gráficas, de funciones, de manejo de data, etc. La documentación de M ATLAB puede ser
leída online:
http://www.mathworks.com/help/matlab/index.html
pero en estas notas introduciremos el manejo de este lenguaje poco a poco resolviendo problemas.
2.2.1 Entrada y Salida de Data

En el siguiente problema aprenderemos a cómo utilizar los siguientes comandos de M ATLAB: input,
disp y fprintf; así como las operaciones aritméticas básicas y los tipos de errores que podemos
cometer cuando diseñamos un programa.
6 Cuando profundicemos más en el curso nos daremos cuenta que este programa puede ser mejorado u optimizado, en
un sentido que estableceremos más adelante.

P ROBLEMA 1 (I NCREMENTO DEL A REA DE UNA E SFERA ): El área de la superficie de una es-
fera con radio r está dada por la fórmula
A(r) = 4π r 2 .
Deseamos averiguar cómo se incrementa el área de la superficie cuando el radio se incrementa

por una pequeña cantidad δ r:
∆ A = 4π(r + δ r)2 − 4π r 2 (2.1a)

∆ A = 4π(2r + δ r) · δ r (2.1b)
∆ A ≈ 8π r · δ r (2.1c)
Las dos primeras fórmulas son exactas, mientras que la tercera es una aproximación que ignora
el término (δ r)2 .
Escribir un programa en que solicite (input) el radio de la esfera r (en kilómetros), el incremen-
to δ r (en milímetros), y exhiba (output) en pantalla el incremento en el área de la superficie (en
metros cuadrados). ¿Cuál es el incremento (∆ A) cuando el radio de un planeta esférico (r = 6367
km) se incrementa por unos pocos milímetros (δ r)? Explore la respuesta a esta pregunta usando
cada una de las tres fórmulas dadas arriba.
Solución del Problema

Un número irracional como π no puede ser almacenado en la computadora, así que convengamos
en aproximarlo a cinco cifras: 3.1416. En el command prompt de M ATLAB, >>, asignamos el número
6367 a la variable r:
>> r = 6367;
Ponemos el punto y coma al final de la sentencia y luego presionamos E NTER para llevar a cabo esta
asignación en la memoria del computador (¿qué sucede si no hubiésemos puesto el punto y coma? 7 ).
Podemos ahora hallar el área de la esfera usando la variable r:
>> 4 * 3.1416 * r ^2
lo cual nos da en pantalla (ventana de comandos o “Command Window”)
ans =
5.0943e+08
Si deseamos mostrar más dígitos podemos usar el comando format con el estilo long:8
>> format long

>> 4 * 3.1416 * r^2
ans =
5.094253814496000e+08
7 R PTA : Después de presionar E NTER, aparte de almacenar la asignación en memoria también se mostrará nuevamente
esta línea en el prompt de M ATLAB. Además, obsérvese que no hemos tipeado una ecuación, sino una asignación, el signo
de igualdad en M ATLAB está dado por ==, y veremos más adelante cómo utilizarlo.
8 para conocer más variaciones de estilo, digitar doc format.
Observemos que M ATLAB utiliza internamente 16 dígitos para sus cálculos (en la mantisa de la
representación punto flotante de un número9 ), pero como en la respuesta final es común usar de 3 a
5 dígitos, regresamos al formato original usando format short.
Si deseamos usar el área calculada en operaciones posteriores, es preferible almacenar la res-
puesta en una variable10 , por ejemplo:
>> A = 4 * 3.1416 * r ^2;
Ahora elaboremos el siguiente programa, el cual resolverá parcialmente nuestro problema. En

M ATLAB existen dos clases de programa: scripts (llamados script M-files en inglés) y funciones (o
function M-files). La diferencia entre ellos es que un programa del tipo función necesita al ser invo-
cado de uno o más inputs y da como resultado uno o más outputs11 . Un script no necesita de input al
ser invocado. En nuestro caso utilizaremos un script y lo llamaremos AreaSup.m (la terminación .m
es añadida por defecto por el software). Comenzamos invocando el editor de M ATLAB para tipear el
programa:
>> edit AreaSup
Se abre una nueva ventana12 y tipeamos el siguiente código:
Programa 2.2: Script AreaSup.m

1 % Programa que calcula el á rea superficial de una esfera
2 % dado el radio provisto por el usuario
3 clc ;
4 r = input ( ' Ingrese radio r : ') ;
5 A = 4* pi * r ^2; % f ó rmula del á rea de una esfera
6 fprintf ( ' Para un radio r = %10.3 f el á rea es A = %10.3 e \ n ' ,r , A ) ;
Para usar este programa hacemos clic izquierdo en el triángulo verde de la parte superior, eti-
quetado con R UN, o simplemente tipeamos el nombre del script en el Command Window
>> AreaSup
e inmediatamente la pantalla es limpiada y aparece
Ingrese radio r:
digitamos el valor 5940.32, presionamos E NTER y obtenemos la siguiente línea:
Para un radio r = 5940.320 el área es A = 4.434e+08
En caso quiera usarse otro valor para r, simplemente hay que invocar el programa de nuevo.
Pasamos a describir lo que hace cada línea.
9 más sobre esto en la sección de Aritmética Punto Flotante.

10 En M ATLAB una variable siempre empieza con una letra y luego pueden usarse sólo letras, números y el símbolo de
subrayado, con un máximo de 64 caracteres. M ATLAB también es sensible a las letras mayúsculas, así que las variables
Area, area, arEa, son todas diferentes.
11 Brindaremos varios ejemplos más adelante para que la diferencia entre scripts y funciones quede más claro. Es posible
construir una función que no necesite input y no dé outputs, este tipo de función tiene sus fines, pero no la usaremos aquí.
12 Dependiendo de la configuración de M ATLAB que uno tenga, es posible que la ventana del Editor aparezca incrustada
en la ventana principal.
Líneas 1-2: Cualquier código escrito a la izquierda del símbolo % es un comentario, y no son vistos
por el compilador a la hora de correr el programa. Los comentarios son útiles para hacer el
código más entendible.
Línea 3: El comando clc borra los contenidos del Command Window, es un acrónimo para clear
command window.
Línea 4: El comando input pide entrada de data al usuario. La sentencia entre signos de apóstrofe
es una data del tipo string,13 la cual aparecerá en pantalla al correr el programa. El radio
ingresado se almacenerá en la variable r.
Línea 5: Se calcula el área de la esfera y el resultado se almacena en la variable A. Obsérvese que

debido al punto y coma al final de la línea, no se exhibe el valor de A en pantalla.
Línea 6: Se exhibe en pantalla la sentencia indicada dentro de los apóstrofes, con los reemplazos de
%10.3f por el valor almacenado en r y %10.3e por el valor almacenado en A. El símbolo \n sig-
nifica E NTER o cambio de línea.
Expliquemos ahora el formato %10.3f. Primero recordemos que esta expresión será reemplazada por
el valor que tenga la variable r, la cual está representando un número real. El número 10 significa
que se reservarán 10 espacios en pantalla para el valor de r, incluyendo el punto decimal. El número
3 significa que se reservarán 3 espacios comenzando desde la derecha para la parte decimal del
valor de r. En caso el valor de r contenga más de 3 cifras decimales, M ATLAB utilizará redondeo
para forzar la exhibición a sólo tres cifras. Quedan entonces 6 espacios libres para exhibir la parte
entera de r. En caso la parte entera contenga más de 6 cifras, entonces se exhibirá TODA la parte
entera de r. Si la parte entera tiene menos de 6 cifras, entonces se usarán espacios en blanco (a la
izquierda del número) para completar los espacios que faltan. Por último, la letra f indica que se
usará notación decimal estándar.14 Otras opciones son %d (para valores enteros) y %e (para exhibir
valores en formato exponencial).
Se recomienda correr el programa para diversos valores de r para entender cómo funciona este
formato.
Ahora nos concentramos en cómo implementar en M ATLAB el incremento del área usando las
tres fórmulas dadas en (2.1a)–(2.1c).
delta_A1 = 4* pi *( r + delta_r ) ^2 - 4* pi * r ^2
delta_A2 = 4* pi *(2* r + delta_r ) * delta_r
delta_A3 = 8* pi * r * delta_r
Las operaciones usadas arriba son parte de las operaciones aritméticas usuales predefinidas en M AT-
LAB (Ver Tabla 2.1). ¿Cuál es el orden de operaciones que efectúa la computadora en la expresión?
4*pi*(r + delta_r)^2 - 4*pi*r^2? Para responder esto es importante recordar la jerarquía en el
orden de operaciones aritméticas (ver Tabla 2.2) que está implícito en las computadoras. El nivel 1
es el que se ejecuta primero, y así sucesivamente subiendo de nivel.
En adición a las operaciones aritméticas básicas existen varias funciones predefinidas en M AT-
LAB . La lista de las funciones más usuales puede ser invocada así
>> help elfun
13 En caso necesitemos saber más sobre un comando determinado, por ejemplo input, podemos usar los comandos help o
doc. Por ejemplo, úsese help input.
14 Es útil observar que se pudo usar simplemente %f. El número 10.3, o el que elija el usuario, simplemente nos brinda
mayor dominio sobre la exhibición de data en pantalla.

Tabla 2.1: Operaciones Aritméticas más usuales en M ATLAB
Símbolo Significado
+ Adición
- Sustracción
* Multiplicación
/ División
^ Exponenciación
Tabla 2.2: Jerarquía en el Orden de Operaciones Aritméticas
Nivel Operaciones
1 Las expresiones dentro de todos los paréntesis son evaluadas, comenzando des-
de lo más interno hacia afuera.
2 Todas las expresiones exponenciales son evaluadas, comenzando de izquierda
a derecha.
3 Todas las multiplicaciones y divisiones son evaluadas, de izquierda a derecha.
4 Todas las adiciones y sustracciones son evaluadas, de izquierda a derecha.
El siguiente programa resuelve el problema propuesto.
Programa 2.3: Script IncrementoAreaSup.m

1 % Este programa calcula el incremento del á rea
2 % superficial de una esfera
3
4 % Adquirir y exhibir la data de entrada
5 r = input ( ' Ingrese radio ( en kil ó metros ) : ') ;
6 delta_r = input ( ' Ingrese incremento ( en mil í metros ) : ') ;
7 clc
8 fprintf ( ' Radio de la esfera = %12.6 f kil ó metros \ n ' ,r )
9 fprintf ( ' Incremento radial = %12.6 f mil í metros \ n \ n ', delta_r )
10 disp ( ' Incremento del Area superficial : ')
11
12 % Convertir de mil í metros y kil ó metros a metros
13 dr = delta_r /10^3;
14 r = r *10^3;
15
16 % M é todo 1
17 delta_A1 = 4* pi *( r + dr ) ^2 - 4* pi * r ^2;
18 fprintf ( '\ n M é todo 1: %15.6 f metros cuadrados \ n ', delta_A1 )
19
20 % M é todo 2
21 delta_A2 = 4* pi *(2* r + dr ) * dr ;
22 fprintf ( ' M é todo 2: %15.6 f metros cuadrados \ n ' , delta_A2 )
23
24 % M é todo 3
25 delta_A3 = 8* pi * r * dr ;
26 fprintf ( ' M é todo 3: %15.6 f metros cuadrados \ n ' , delta_A3 )
Para los valores de r = 6367 y δ r = 1.234 obtenemos el siguiente output:
Radio de la esfera = 6367.000000 kilómetros

Incremento radial = 1.234000 milímetros
Incremento del Area superficial:
Método 1: 197464.812500 metros cuadrados

¿Por qué la diferencia en los tres métodos? Es más, los métodos 1 y 2 deberían darnos la misma res-
puesta, ya que representan el mismo valor en aritmética exacta. Resulta que la computadora realiza
otro tipo de aritmética, llamada aritmética punto flotante o aritmética computacional, la cual conlleva
a errores, dado que no todo número real puede representarse exactamente en la computadora.
Tipos de Errores
Cuando usamos la computadora hay que tener en cuenta los tipos de errores que debemos minimizar,
con la finalidad de poder confiar en los resultados que nos dé la computadora. Son dos los tipos de
errores:
Error de redondeo: Ocurren debido a que las computadoras tienen una capacidad limitada para
representar exactamente un número. Por ejemplo, el número irracional π.
Error de truncación: Ocurren cuando formulaciones matemáticas exactas son sustituidas por
aproximaciones. Por ejemplo, en el script que calcula el incremento en el área de la superfi-
cie de una esfera, observamos que si bien el método 3 no es la fórmula exacta, nos brinda una
buena aproximación del incremento de área.
Aprenderemos más adelante cómo evitar la magnificación del error de redondeo en un programa.
En el siguiente problema ilustraremos el uso de una estructura condicional, la cual permite que
un bloque de código sea o no ejecutado, si cierta proposición que lo precede es verdadera o falsa.
P ROBLEMA 2 (U SO DE C ONDICIONALES : IF - ELSEIF ): Una función cuadrática de la forma f (x) =

x2 + bx + c alcanza su mínimo valor en el punto crítico x c = − b/2. Calcule el valor mínimo de f en
un intervalo [L, R] para el cual se requiere cierta verificación: si x c está en el intervalo, entonces
f (x c ) es la respuesta; caso contrario, f es minimizado en x = L o bien en x = R.
Escriba un script que solicite los números reales L, R, b y c, e imprima en pantalla el valor
mínimo de la función cuadrática en el intervalo [L, R]. El script debe también exhibir el valor de
x donde el mínimo ocurre.

Para adquirir un sentido geométrico del problema, describimos las tres alternativas posibles:
R ≤ xc
Mínimo valor = f (R)
L ≤ xc ≤ R
Mínimo valor = f (x c )
xc ≤ L
Mínimo valor = f (L)
Dadas estas tres posibilidades, el siguiente seudocódigo identifica y procesa la alternativa correc-
ta
Si R < x c
Imprimir f (R ) y R
sino si L ≤ x c ≤ R
Imprimir f ( x c ) y x c
sino
Imprimir f (L) y L
fin
Para implementar este seudocódigo en M ATLAB usaremos el condicional if-elseif y expresiones

lógicas o booleanas, esto es, expresiones a las que puede asignarse un valor de verdad (verdadero o
falso). Comenzamos listando los operadores de relación definidos en M ATLAB:
Símbolo Significado
< menor que
<= menor o igual que
== igual
> mayor que
>= mayor o igual que
~= no igual
En M ATLAB una proposición verdadera toma el valor 1, y una falsa el valor cero. Por ejemplo, sabe-
mos que 5<8 es verdadero, si almacenamos esto en una variable y luego invocamos su valor, veremos
que obtenemos 1:
>> a = 5 < 8;
>> a
a =
1
Cuando los operadores de relación son usados entre dos operandos numéricos siempre nos dan como
resultado un valor lógico, verdadero (1) o bien falso (0). Es importante mencionar ya desde ahora que
debemos tener cierto cuidado al usar los operadores == y ~= al comparar dos valores numéricos. Esto
se debe a los errores de redondeo al aproximar números reales. Estudiemos el siguiente ejemplo
>> x = 0;
>> y = sin(pi);
>> x == y
ans =
0
Esto nos dice que sen(π) 6= 0 según la computadora. ¿A qué se debe esto? Si mostramos en pantalla
el valor de sen(π), según M ATLAB, obtendremos 1.2246 × 10−16 , lo cual es muy pequeño y cercano
a cero, pero no es exactamente cero. Esto se debe a que la función seno no es evaluada en el valor
exacto de π sino en una aproximación decimal de este número en la computadora, con un número
máximo de cifras significativas. Recordar que al ser π un número irracional, no puede representarse
con un número finito de cifras decimales, ni como el cociente de dos enteros.
Ahora bien, dado que pedir igualdad exacta entre dos valores reales no es recomendable ¿cómo
verificar entonces una igualdad? Se sugiere mejor verificar una cercanía o distancia mínima. Por
ejemplo, podríamos asumir que x y y son iguales en la computadora si la distancia entre ellos es
menor que 10−14 . Esto lo digitamos así
>> abs(x-y) < 1e-14

ans =
1
A continuación presentamos los operadores booleanos o lógicos. Estos operan sobre uno o dos
operandos lógicos (proposiciones) y conducen a un resultado lógico. Existen cinco operadores binarios
lógicos: Y (& y &&), O inclusiva (| y ||), y O exclusiva (xor); y un operador lógico unitario: NO
(~), véase la Tabla 2.3. Observemos que existen dos versiones para los operadores Y y O inclusivo.
¿Cuándo es usada la versión con atajo en evaluación? Para responder esto debemos recordar la tabla
de verdad de estos operadores y tener en cuenta que este tipo de atajo solo se usa entre expresiones
Tabla 2.3: Operadores Lógicos definidos en M ATLAB
Símbolo Operación
& Y lógico
&& Y lógico con atajo en evaluación
| O lógico inclusivo
|| O lógico inclusivo con atajo en evaluación
xor O lógico exclusivo
∼ NO lógico
Tabla 2.4: Tabla de Verdad de los Operadores Lógicos
Inputs Y O XOR NO
x y x&y x && y x|y x || y xor(x,y) ∼x
0 0 0 0 0 0 0 1
0 1 0 0 1 1 1 1
1 0 0 0 1 1 1 0
1 1 1 1 1 1 0 0
lógicas 1 × 1 ó escalares15 (véase la Tabla 2.4). Es claro que las versiones con atajo tienen la misma
tabla de verdad, pero la conveniencia radica en que a veces no necesitamos verificar el valor de verdad
de una expresión para obtener el resultado lógico de verdadero o falso.16 En la siguiente tabla se ha
puesto un signo de interrogación en los lugares donde la computadora no necesita averiguar el valor
de verdad de la segunda proposición:
x y x && y x y x||y
0 ? 0 0 0 0
0 ? 0 0 1 1
1 0 0 1 ? 1
1 1 1 1 ? 1
¿Qué se gana con esto? TIEMPO, y el tiempo de ejecución de un programa es primordial cuando
resolvemos problemas muy elaborados, sobre todo si nuestro programa tiene que realizar muchos
cálculos intermedios antes de obtener la respuesta final.
En la jerarquía de operaciones, los operadores lógicos son evaluados después de que todas las
operaciones aritméticas y todos los operadores relacionales han sido evaluados. Ver Tabla 2.5. Como
ocurre con las operaciones aritméticas, paréntesis pueden ser usados para cambiar el orden por
defecto en la evaluación.
Ahora mostraremos cómo implementar en Matlab el algoritmo mencionado al principio de la so-
lución de este problema. Primero pedimos al usuario los extremos del intervalo [L, R], los coeficientes
b y c de la función cuadrática, y calculamos el punto crítico x c = − b/2:
b = input ( ' Ingrese b : ') ;

c = input ( ' Ingrese c : ') ;
L = input ( ' Ingrese L : ') ;
R = input ( ' Ingrese R ( L < R ) : ');
15 Cuando expresiones lógicas están dentro de un arreglo (array, en inglés) no escalares, los operadores lógicos solo podrán
usarse componente a componente. Comprenderemos mejor esta aseveración más adelante.
16 Si analizamos la Tabla 2.4 observamos que para el conectivo lógico Y, el resultado lógico sólo es verdadero (1) si ambas
proposiciones x y y son verdaderas, en los demás casos es falso (0), mientras que para el conectivo O inclusivo el resultado
sólo es falso si ambas proposiciones son falsas.
Tabla 2.5: Jerarquía en el Orden de Operadores Lógicos
Nivel Operaciones
1 Todos los operadores aritméticos son evaluados en el orden previamente descrito.
2 Todos los operadores relacionales (==, ∼=, >, >=, <, <=) son evaluados, de
izquierda a derecha.
3 Todas los operadores de negación (∼) son evaluados.
4 Los operadores & y && son evaluados, de izquierda a derecha.
5 Los operadores |, || y xor son evaluados, de izquierda a derecha.
xc = -b /2;
Ahora calculamos f (L), f (R) y f (x c ):
fL = L ^2 + b * L + c ;
fR = R ^2 + b * R + c ;
fxc = xc ^2 + b * xc + c ;
y a continuación traducimos en el lenguaje de M ATLAB el seudocódigo mencionado al principio:
if R < xc
fprintf ( ' Valor m í nimo f ( R ) = %5.2 f \ n ' , fR ) ;
fprintf ( ' alcanzado en R = %5.2 f \ n ',R ) ;
elseif L <= xc && xc <= R
fprintf ( ' Valor m í nimo f ( xc ) = %5.2 f \ n ' , fxc ) ;
fprintf ( ' alcanzado en xc = %5.2 f \ n ' , xc ) ;
else
fprintf ( ' Valor m í nimo f ( L ) = %5.2 f \ n ' , fL ) ;
fprintf ( ' alcanzado en L = %5.2 f \ n ',L ) ;
end
En este último bloque de código hemos visto cómo se utiliza la construcción if-elseif. En general,
toma la siguiente forma:
if " expresi ó n l ó gica 1"

bloque 1 ( a ejecutar si expresi ó n l ó gica 1 es verdadera )
elseif " expresi ó n l ó gica 2"
bloque 2 ( a ejecutar si expresi ó n l ó gica 2 es verdadera )
else
bloque 3 ( a ejecutar si ninguna de las expresiones ...
l ó gicas anteriores son verdaderas )
end
Existen otras variantes de la construcción if-elseif, como if-else o simplemente if.

Antes de presentar en un solo sitio el script que resuelve el problema, sería bueno hacer una pau-
sa aquí e intentar elaborar por si mismo un script. Luego compare con la solución que presentamos
a continuación.
Programa 2.4: Script ValorMinimo.m

1 % Script ValorMinimo
2 % Calcula el m í nimo de una funci ó n cuadr á tica de la forma
3 % x ^2 + bx + c en el intervalo [L , R ].
4
5 % Adquisici ó n y exhibici ó n de la data de entrada ( input )

6 b = input ( ' Ingrese b : ') ;
7 c = input ( ' Ingrese c : ') ;
8 L = input ( ' Ingrese L : ') ;
9 R = input ( ' Ingrese R (L < R ) : ') ;
10 clc
11 fprintf ( ' Cuadr á tica : x ^2 + bx + c , b = %5.2 f , c = %5.2 f \ n ',b , c ) ;
12 fprintf ( ' Intervalo : [L , R ] , L = %5.2 f , R = %5.2 f \ n \ n ',L , R ) ;
13 % C á lculo del punto cr í tico
14 xc = -b /2;
15 if xc < L
16 % M í nimo se encuentra en el extremo izquierdo de [L , R ]
17 fL = L ^2 + b * L + c ;
18 fprintf ( 'M í nimo en x = %5.2 f \ n ' ,L )
19 fprintf ( ' Valor m í nimo f ( x ) = %5.2 f \ n ' , fL )
20 elseif L <= xc && xc <= R
21 % M í nimo se encuentra en el punto cr í tico
22 fxc = c - ( b /2) ^2;
23 fprintf ( 'M í nimo en x = %5.2 f \ n ' , xc )
24 fprintf ( ' Valor m í nimo f ( x ) = %5.2 f \ n ' , fxc )
25 else
26 % Minimo se encuentra en el extremo derecho de [L , R ]
27 fR = R ^2 + b * R + c ;
28 fprintf ( 'M í nimo en x = %5.2 f \ n ' ,R )
29 fprintf ( ' Valor m í nimo f ( x ) = %5.2 f \ n ' , fR )
30 end
Funciones adicionales en MATLAB

Se recomienda usar la palabra help o doc para profundizar en el uso de las siguientes funciones:
• max y min: El valor de
max(expresion1,expresion2)
es el más grande de los dos valores, obtenidos al evaluar las expresiones. La función min trabaja
de manera similar. Ambas funciones se pueden aplicar a una lista de dos o más expresiones.
• floor,ceil,round: Estas son funciones que son aplicadas a números reales para obtener valores
enteros cercanos. Así, para un número decimal x, tenemos que floor(x) nos da el mayor número
entero menor o igual que x, ceil(x) nos da el menor número entero mayor o igual que x, y round←-
(x) nos da el número entero más cercano a x. En caso x esté localizado exactamente entre dos
enteros, round(x) nos dará el mayor de los enteros. Por ejemplo, round(4.5) nos da 5.
• rem: Si x y y son enteros positivos, entonces rem(x,y) es el resto de dividir x por y. Por ejemplo,
rem(23,3) nos da 2. La función rem también puede aplicarse a números reales de cualquier signo.
2.2.2 Problemas Propuestos

1) Una temperatura puede convertirse de la escala Farenheit a la Celsius usando la fórmula
5
c = ( f − 32),
9
donde f es una temperatura dada en ◦ F y c es la temperatura convertida a ◦ C. Escriba un script

que pida una temperatura en ◦ F , haga la conversión a grados Celsius y exhiba en pantalla esta
temperatura.
2) Un elipsoide, tal como nuestro planeta Tierra, se obtiene al rotar una elipse alrededor de su eje
menor. Se sabe que el radio ecuatorial de la Tierra es aproximadamente 21 km más largo que su
radio polar. El área superficial de un elipsoide está dado por17
r 22 ³ cos γ ´¶
µ ¡ ¢
2
A(r 1 , r 2 ) = 2π r 1 + ln
sen(γ) 1 − sen(γ)
donde r 1 es el radio ecuatorial, r 2 es el radio polar, y
³r ´
2
γ = arc cos .
r1
Asumimos r 2 < r 1 . Escriba un script que lea los radios ecuatorial y polar y exhiba en pantalla
A(r 1 , r 2 ) y la aproximación
³ r + r ´2
1 2
4π .
2
Aplique el script a los datos terrestres (r 1 , r 2 ) = (6378.137, 6356.752). Exhiba suficientes dígitos
para mostrar la discrepancia entre la fórmula exacta y la aproximación.
3) Una elipse con semiejes a y b está especificada por

³ x ´2 ³ y ´2
+ = 1.
a b
Si r = a = b entonces esto define una circunferencia cuyo perímetro está dado por P = 2π r. Desafor-
tunadamente, si a 6= b entonces no hay una fórmula simple para el perímetro y debemos recurrir
a una aproximación. Se conocen varias posibilidades:18
³ 3h ´
P1 = π(a + b) P5 = π(a + b) 1 + p
10 + 4 − 3h
p 64 − 3h2
P2 = π 2(a2 + b2 ) P6 = π(a + b)
64 − 16h
s
(a − b)2 256 − 48h − 21h2
P3 = π 2(a2 + b2 ) − P7 = π(a + b)
2 256 − 112h + 3h2
³ h ´2 ³ 3 − p1 − h ´
P4 = π(a + b) 1 + P8 = π(a + b)
8 2
Aquí,
³ a − b ´2
h= .
a+b
Escriba un script que solicite a y b e imprima en pantalla los valores P1 , . . . , P8 en una manera
que facilite la comparación. Los valores de h también deben ser exhibidos. Intente los siguientes
valores de entrada
(a, b) = (1, 1), (1, 0.9), . . . , (1, 0.1).
¿Qué puede afirmar acerca de las diferencias entre las fórmulas del perímetro cuando la elipse se
vuelve más ovalada en forma?
17 Recordar que M ATLAB está en inglés, y por ende, la función seno es representada por sin. Además, la función logaritmo
natural, escrita comúnmente como ln(x), es escrita en M ATLAB como log(x).

18 Véase http://www.mathsisfun.com/geometry/ellipse-perimeter.html
4) Suponga que el valor de x es un entero positivo. Escriba una expresión booleana que es verdadera
si x es divisible por 2, 5 y 7.
5) Escriba un script que solicite números reales L y R e imprima el valor máximo de f (x) = cos(x)
en el intervalo [L, R]. Asuma L < R. Recuerde que el coseno alcanza un valor máximo de uno en
múltiplos enteros de 2π.
6) Considere la función cúbica

q(x) = ax3 + bx2 + cx + d, a 6= 0.
Decimos que q es simple si sus tres raíces son reales y distintas. Decimos que q es monótona si es
siempre creciente o bien siempre decreciente. Como ejemplo, consideremos
q 1 (x) = (x − 1)(x − 2)(x − 3) = x3 − 6x2 + 11x − 6

q 2 (x) = (x − 1)(x − 2)(x − 3) + 100 = x3 − 6x2 + 11x + 94
q 3 (x) = − x(x2 + 1) = − x3 − x
entonces q 1 es simple pero no monótona, q 2 no es simple ni monótona, y q 3 es monótona pero

no simple. En este problema debe escribir un script que puede ser usado para verificar estas
situaciones. Como sugerencia, observe que la derivada q0 (x) y sus raíces r 1 y r 2 son importantes
de tener en cuenta.
El script debe comenzar pidiendo al usuario los cuatro coeficientes a, b, c y d. Si r 1 y r 2 no son
reales y distintos, entonces q es monótona y el mensaje Función Monótona debe ser impreso en
pantalla. Caso contrario, los valores de r 1 , q(r 1 ), r 2 y q(r 2 ) deben ser impresos. En este último
caso q no es monótona y podemos deducir de los cuatro valores calculados si es o no simple. En
particular, q es simple si y sólo si q(r 1 ) q(r 2 ) < 0.
2.2.3 Bucles For y While

Muchas veces es necesario realizar un proceso iterativo varias veces. Si muy bien se podría realizar
esto a mano, resulta aburrido y monótono, y es mucho mejor dejar este trabajo a la computadora, la
cual no se va a quejar.
Estudiaremos dos maneras de iterar en M ATLAB.
Bucle for: Es una construcción que permite repetir un bloque de código un número determinado de
veces. La repetición está controlada por los valores que toma el índice. Su forma más simple
es:
for indice = valor_inicial:valor_final

Bloque de código
end
Es posible usar un bucle for dentro de otro. La necesidad de esto se verá más adelante.
Bucle while: Es una construcción que permite repetir un bloque de código un número indetermi-
nado de veces, hasta que se cumpla cierta condición. Toma la forma
while condición_es_verdadera
Bloque de código
end
Cada vez que se repite el bloque de código el programa verifica si la condición sigue siendo
verdadera. Esto implica que dentro de este bloque debe existir alguna expresión que llegará a
cambiar el valor de verdad de la condición.
Ilustraremos el uso de estos tipos de bucle con los siguientes problemas resueltos.
P ROBLEMA 3 (E MBALDOSADO DE UN D ISCO. U SO DEL B UCLE F OR ): Supongamos que n es

un entero positivo y dibujemos el círculo x2 + y2 = n2 en un papel cuadriculado con cuadrados de
área unitaria y lado 1. La Fig. 2.1 exhibe el caso n = 10. Nótese que el área del disco es π n2 y cada
baldosa no cortada tiene área unitaria. Si hay N baldosas sin cortar, entonces concluimos que
N ≈ π n2
ya que las baldosas sin cortar casi cubren el disco. Escriba un script que solicite al usuario un
entero n y exhiba la aproximación del número π
N
ρn =
n2
junto con el error ¯ρ n − π¯.

¯ ¯
Figura 2.1: Embaldosado del disco x2 + y2 ≤ n2 .

Para simplificar la solución podemos usar simetría. Cada cuadrante tiene exactamente el mismo
número de baldosas sin cortar. Así, solo necesitamos contar el número de baldosas sin cortar en el
primer cuadrante, como exhibimos en la Fig. 2.2. Si N1 es el número de baldosas sin cortar en el
primer cuadrante, entonces ρ n = 4N1 /n2 .
Para calcular N1 sumamos el número de baldosas sin cortar que están localizadas en cada fila
horizontal. Refiriéndonos a la Fig. 2.2, vemos que hay nueve baldosas sin cortar en la fila 1 (row 1,
en inglés), nueve baldosas sin cortar en la fila 2, etc. Para un n genérico, procedemos como sigue:
Inicializar : N1 = 0
Repetir para filas 1 hasta la n
Calcular el n ú mero de baldosas sin cortar en la fila .
A ñ adir este n ú mero a N1 .
Poner ρ n = 4 N1 / n2 .
Figura 2.2: Embaldosado del Primer Cuadrante, caso n = 10.
Dado que sabemos el número de veces que se va a repetir nuestro programa, utilizaremos un
bucle for. Refinamos nuestro seudocódigo así:
Asumir que n est á inicializada y poner N1 a cero

for k =1: n
Calcular el n ú mero de baldosas sin cortar en la fila k .
A ñ adir este n ú mero a N1 .
end
Poner ρ n = 4 N1 / n2
Ahora nos concentramos en cómo calcular el número de baldosas sin cortar en la fila k. En la Fig. 2.2
recordemos que el origen (0, 0) se encuentra en la esquina inferior izquierda. El punto (0, 1) (relativo
a la fila 1), se encuentra obviamente encima del origen. La recta horizontal que pasa por (0, 1), esto
es, y = 1, interseca a la circunferencia en el punto
³p ´
n 2 − 12 , 1 .
Supongamos que (h, 1) son las coordenadas del extremo superior derecho del último cuadradito de la
fila 1 que deseamos saber si es cortado o no por la circunferencia x2 + y2 = n2 . Recuérdese que h es
un número entero positivo. Para que este último cuadradito sea considerado como baldosa sin cortar,
debe suceder que p
h ≤ n 2 − 12 ,
o lo que es mejor, h = floor(sqrt(n^2-1^2)). Podemos generalizar esta idea a la fila k, y obtener

que el número de baldosas sin cortar, en esa fila, debe ser
h = floor ( sqrt ( n ^2 - k ^2) )
Ahora pensemos en cómo calcular N1 . Antes del bucle for debemos inicializar N1 = 0. Para actuali-
zar este valor dentro del bucle for basta escribir
N1 = N1 + h
Esta sentencia se lee “el valor de N1 se actualiza con el valor que resulte de calcular N1 + h”. Es
importante recordar que el símbolo = no significa igualdad cuando programamos en M ATLAB.
Resumiendo, formamos el siguiente pedazo de código
% asumiendo que ya hemos provisto el valor de n

N1 = 0;
for k = 1: n
h = floor ( sqrt ( n ^2 - k ^2) ) ;
N1 = N1 + h ;
end
ro_n = 4* N1 / n ^2;
Aquí, la variable k, llamada índice del bucle, toma valores del 1 al 10:
>> k = 1:10
k =
1 2 3 4 5 6 7 8 9 10
Podemos cambiar el conjunto de valores que toma el índice proveyendo un incremento, por ejemplo
>> k = 1:2:10
k =
1 3 5 7 9
En general, esta construcción, llamado arreglo (unidimensional), puede construirse así:
nombre_arreglo = valor_inicial:incremento:valor_final
Aquí, el incremento puede ser cualquier número real positivo y asumimos que el valor_inicial es
menor o igual que el valor_final. Estos valores pueden ser números decimales. Por ejemplo:
>> k = -2.35:1.4:4.5
k =
-2.3500 -0.9500 0.4500 1.8500 3.2500
Si el valor_inicial es mayor que el valor_final obtenemos un arreglo vacío. Por otro lado, es
posible también usar un decremento si deseamos un arreglo ordenado en forma decreciente:
nombre_arreglo = valor_inicial:decremento:valor_final
El siguiente script resuelve el problema.
Programa 2.5: Script AproximaPi.m

1 % Este programa aproxima el n ú mero irracional pi usando
2 % conteo de á reas de cuadrados inscritos en un disco
3
4 % Ingrese el radio del disco
5 clc
6 n = input ( ' Ingrese un radio con valor entero n : ') ;
7
8 % Embaldosar solo el primer cuadrante y luego multiplicar por cuatro
9 N1 = 0;
10 for k = 1: n
11 % calcular el n ú mero de baldosas " sin cortar " en fila k
12 h = floor ( sqrt ( n ^2 - k ^2) ) ;
13 N1 = N1 + h ;
14 end
15
16 % c á lculo de la aproximaci ó n de pi
17 ro_n = 4* N1 / n ^2;
18
19 % Exhibici ó n del estimado y del error
20 clc
21 fprintf ( 'n = %12 d \ n ',n ) ;
22 fprintf ( ' ro_n = %12.8 f \ n ' , ro_n ) ;
23 fprintf ( ' error = %12.8 f \ n ' , abs (pi - ro_n ) ) ;
Por ejemplo, para n = 104 obtenemos
n = 10000
ro_n = 3.14119052
error = 0.00040213
P ROBLEMA 4 (P OLÍGONOS INSCRITOS Y CIRCUNSCRITOS . U SO DEL B UCLE W HILE ): Con- si-

deremos n puntos igualmente espaciados sobre el círculo unitario x2 + y2 = 1. Si conectamos estos
puntos en orden obtendremos un polígono regular inscrito de n lados. Pero si conectamos las rec-
tas tangentes en estos puntos obtendremos un polígono regular circunscrito de n lados. La Fig. 2.3
exhibe estas construcciones para el caso n = 6. El área del círculo unitario es π. Las áreas de los
polígonos inscrito y circunscrito están dadas por
n ³ 2π ´
An = sen ,
2 n
³π´
B n = n tan ,
n
respectivamente. Es obvio que para todo n
A n < π < Bn.
Usando la Fig. 2.3 para guiarnos, podemos deducir que el promedio de las dos áreas poligonales
A n + Bn
ρn =
2
es una aproximación de π, con error absoluto que satisface
¯ρ n − π¯ < B n − A n .
¯ ¯
Claramente el error converge a cero cuando n → ∞.

Escriba un script que solicite al usuario una tolerancia positiva real δ y exhiba el valor de ρ m ,
donde m es el entero más pequeño tal que
(B m − A m ) ≤ δ.
Esto garantiza que ρ m esté a una distancia menor que δ de π.

Es útil comenzar exhibiendo un script que revele la calidad de la aproximación ρ n cuando n crece:
Figura 2.3: Hexágonos inscritos y circunscritos.
n = 10;
for k = 1:8
A_n = (n /2) * sin (2* pi / n ) ;
B_n = n* tan ( pi / n ) ;
ro_n = ( A_n + B_n ) /2;
n = 10* n ;
end
Los valores de ensayo son n = 10, 102 , 103 , . . . , 108 . Si añadimos más líneas de código para imprimir
en pantalla, conseguimos:
n A_n B_n |ro_n - pi|

---------------------------------------------------------------
1.00e+01 2.9389262614623659 3.2491969623290631 4.75e-02
1.00e+02 3.1395259764656687 3.1426266043351152 5.16e-04
1.00e+03 3.1415719827794755 3.1416029890561563 5.17e-06
1.00e+04 3.1415924468812859 3.1415927569440529 5.17e-08
1.00e+05 3.1415926515227079 3.1415926546233357 5.17e-10
1.00e+06 3.1415926535691225 3.1415926536001288 5.17e-12
1.00e+07 3.1415926535895866 3.1415926535898966 5.15e-14
1.00e+08 3.1415926535897909 3.1415926535897944 4.44e-16
Podemos deducir que para un n dado, el error es aproximadamente 1/n2 .

Nuestra investigación con el bucle for es interesante, pero no nos ayuda directamente con el
problema a mano. Debemos encontrar el mínimo valor de n tal que (B n − A n ) < δ. Supongamos que
este mínimo ocurre en n = m. Como no sabemos cuándo sucederá esto, podemos acercarnos poco a
poco aumentando el valor de n. El siguiente seudocódigo ilustra el proceso a seguir:
Ingrese la tolerancia ( error ) δ y ponga n = 3 .

Calcule A 3 , B3 y la cota del error B3 − A 3 .
Repetir el siguiente bloque mientras la cota del error sea mayor que δ :
Incrementar n por uno .
Actualizar A n , B n y la cota del error B n − A n .
Poner m = n y exhibir ρ m .
Un bucle for no puede ser usado porque debemos saber el número de iteraciones por adelantado
y esto es justo lo que buscamos. Un mecanismo alternativo para iterar es el bucle while. Su uso en el
seudocódigo anterior se implementa en M ATLAB como sigue:
delta = input ( ' Ingrese delta : ') ;

n = 3; A_n = ( n /2) * sin (2* pi / n ) ; B_n = n * tan ( pi / n ) ;
CotaError = B_n - A_n ;

while CotaError > delta
n = n + 1;
A_n = ( n /2) * sin (2* pi / n ) ;
B_n = n * tan ( pi / n ) ;
CotaError = B_n - A_n ;
end
m = n ; ro_m = ( A_n + B_n ) /2;
Si la expresión lógica CotaError > delta es verdadera, entonces el cuerpo del bucle es ejecuta-
do. Esta ejecución se repetirá hasta que la expresión lógica sea falsa. Cuando el bucle termina de
ejecutarse, el programa continuará ejecutando las líneas debajo de la palabra end. Debemos tener
presente que variables que sean utilizadas por el bucle while deben inicializarse antes de comenzar
el bucle.
La solución final del problema es presentada en el siguiente script.
Programa 2.6: Script AproximaPiConBucleWhile.m

1 % Este programa aproxima el n ú mero pi con una tolerancia dada por el
2 % usuario , en base a á reas de pol í gonos inscritos y circunscritos a una
3 % circunferencia .
4
5 % Ingresar los par á metros de iteraci ó n
6 clc
7 delta = input ( ' Ingrese la tolerancia del error : ') ;
8 nMax = input ( ' Ingrese n ú mero m á ximo de iteraciones : ') ;
9
10 % Caso primario : triangular
11 n = 3; % N ú mero de lados del pol í gono
12 A_n = ( n /2) * sin (2* pi / n ) ; % Area inscrita
13 B_n = n * tan ( pi / n ) ; % Area circunscrita
14 CotaError = B_n - A_n ; % Cota del error
15
16 % Iterar mientras que el error sea demasiado grande y n peque ñ o
17 while ( CotaError > delta ) && ( n < nMax )
18 n = n +1;
19 A_n = ( n /2) * sin (2* pi / n ) ;
20 B_n = n * tan ( pi / n ) ;
21 CotaError = B_n - A_n ;
22 end
23
24 % Exhibir la aproximaci ó n final
25 m = n;
26 ro_m = ( A_n + B_n ) /2;
27 clc
28 fprintf ( ' delta = %10.3 e \ n m = %1 d \ n nMax = %1 d \ n \ n ' , delta ,m , nMax ) ;
29 fprintf ( ' ro_m = %20.15 f \ n Pi = %20.15 f \ n ', ro_m , pi ) ;
Observemos que la condición al inicio del bucle while en un principio es verdadera (caso contrario
no se ejecutaría iteración alguna), y basta que una de las proposiciones encerradas entre paréntesis
llegue a ser falsa para que el bucle termine. Si no se hubiese tenido el cuidado de incorporar un
número máximo de iteraciones (nMax) podría haberse incurrido en un bucle infinito (esto ocurriría si
damos, por ejemplo, un δ muy pequeño como 10−20 ). Recordar que en caso el programa tarde mucho
en exhibir el output, debe usarse Ctrl-c para abortar la ejecución.
Para el input δ = 10−6 y nMax=104 obtenemos la siguiente salida
delta = 1.000e-06
m = 5569
nMax = 1000000
ro_m = 3.141592486963389
Pi = 3.141592653589793
Recomendaciones para una Buena Programación

Aunque las siguientes reglas pueden obviarse, recomendamos seguirlas con el fin de evitar posibles
errores dentro de un programa, sobre todo cuando éste se vuelve más largo y complejo.
• Indentar los bloques de código dentro del bucle respectivo. Esto permitirá una lectura más
fácil del programa, no solo para uno mismo, sino para que sea entendido por otras personas.
• Nunca modificar el valor del índice dentro del bloque de un bucle for. Este error común de
programador novato trae consecuencias desastrosas, y cuando son cometidos, a veces son difíciles
de ubicar.
• Decidir si la iteración es mejor manipulada por un bucle for o un bucle while. Si el
número de iteraciones es conocido por adelantado, entonces un bucle for es apropiado.
• Asegurarse que el rango de valores de un bucle for es el correcto. Con frecuencia, probar
primero con un rango de valores pequeño ayuda para verificar si la idea que queremos llevar a
cabo es correcta.
Comandos Adicionales para uso en bucles: BREAK y CONTINUE

Podemos controlar las iteraciones de los bucles for y while con los siguientes dos comandos que se
colocan en el cuerpo del bucle:
break: Cuando es ejecutada la línea que contiene al comando break, el bucle es obligado a terminar
y el control del programa pasa a la línea de código ubicada después del fin (end) del bucle.
continue: Cuando es ejecutado este comando el bucle vuelve a reiniciarse con el próximo valor de
la iteración.
El siguiente programa ilustra el uso del comando break.
Programa 2.7: Script UsoDeBreakConFor.m

1 % uso de break dentro del bucle for
2 clc
3 x = 0;
4 for i =1:5
5 x = x + 1;
6 if i == 3
7 break ;
8 end
9 fprintf ( ' Para iteraci ó n i = %1d , el valor de x es %1 d \ n ',i , x ) ;
10 end
11 disp ( ' Este aviso aparece despu é s de terminado el bucle ') ;
12 fprintf ( ' El valor actual de x en memoria es %1 d \ n ' ,x ) ;
Si compilamos el programa obtenemos
Para iteración i = 1, el valor de x es 1

Este aviso aparece después de terminado el bucle
El valor actual de x en memoria es 3
Queda como ejercicio justificar porqué el valor en memoria para x es 3. El siguiente programa ilustra
el uso del comando continue.
Programa 2.8: Script UsoDeContinueConFor.m

1 % uso de continue dentro del bucle for
2 clc
3 x = 0;
4 for i =1:5
5 if i == 3
6 continue ;
7 end
8 x = x + 1;
9 fprintf ( ' Para iteraci ó n i = %1d , el valor de x es %1 d \ n ',i , x ) ;
10 end
11 disp ( ' Este aviso aparece despu é s de terminado el bucle ') ;
12 fprintf ( ' El valor actual de x en memoria es %1 d \ n ',x ) ;
Si compilamos el programa obtenemos

Este aviso aparece después de terminado el bucle
El valor actual de x en memoria es 4
2.2.4 Problemas Propuestos

1) Usando diversos valores de n en el programa aproximapi.m, ¿en cuánto varía el error si n←-
es incrementado por un factor de 10?, ¿cuánto se incrementaría el tiempo de ejecución si n←-
es incrementado por un factor de 10? Use el comando tic ... toc para calcular el tiempo de
ejecución.
2) Modifique el script aproximapi.m para que exhiba el valor de ¯π − ρ n ¯ para n = 10 j , j = 1:8.

¯ ¯
3) Para n suficientemente grande se saben que
1 (−1)n+1 X n (−1) k+1 π

Rn = 1 − +···− = ≈ ,
3 2n − 1 k=1 2k − 1 4
1 1 n 1
X π2
Tn = 1 + + · · · + = ≈ ,
22 n2 k=1 k2 6
1 1 n 1
X π4
Un = 1 + + · · · + = ≈ ,
24 n4 k=1 k4 90
los cuales nos dan tres maneras diferentes para estimar π:

p p
ρn = 4 Rn, τn = µn =
4
6 Tn, 90Un .
Escriba un script que exhiba los valores de ¯π − ρ n ¯, |π − τn |, y ¯π − µn ¯ para n = 100, 200, 300, . . . , 1000.
¯ ¯ ¯ ¯
4) Escriba un script que verifique las desigualdades
2 p Xn p 4n + 3 p
n n≤ k≤ n
3 k=1 6
para n = 1, . . . , 100.
5) Modifique el programa 2.6 de manera que termine el bucle cuando | A n+1 − A n | ≤ δ ó |B n+1 − B n | ≤
δ.
6) Añada comandos apropiados (fprintf o disp) al fragmento de código (bucle for) de la pág. 34,
para obtener la tabla de valores mostrada en esa página, para n = 10, 102 , . . . , 108 .
7) La implementación del programa 2.6 requiere la constante predefinida pi. Usando repetidamente
las fórmulas de ángulo mitad
s s
1 + cos(θ ) 1 − cos(θ )
cos(θ /2) = , sen(θ /2) = ,
2 2
podemos evitar esta dependencia y producir un método genuino que “descubra” π. Observe que el
caso n = 3 involucra
p p
cos(π/3) = 1/2, sen(π/3) = 3/2, y sen(2π/3) = 3/2.
Las evaluaciones necesarias para coseno y seno para A 6 y B6 pueden ser obtenidas de estos va-
lores vía las fórmulas de ángulo mitad. Usando estas evaluaciones trigonométricas actualizadas,
podemos obtener los valores necesarios de coseno y seno para A 12 y B12 , etc. Use esta idea para
escribir un script que pida un número positivo δ > 10−12 y calcule un estimado de π con un error
menor que δ.
2.2.5 Ploteo en Dos Dimensiones

En esta sección aprenderemos a cómo manipular arreglos, tener acceso a sus componentes, y plotear
funciones reales de variable real.
P ROBLEMA 5: Escriba un script que exhiba la gráfica de la función
sen(5x) e− x/2
f (x) =
1 + x2
en el intervalo [−2, 3].

Primero ilustraremos cómo plotear una función más simple g(x) = sen(x) en el intervalo [0, 2π]. Debe-
mos entender que la manera cómo funciona el ploteo en M ATLAB (y en todos los software) es dibujar
en el plano un conjunto finito de puntos y luego unirlos con segmentos de recta. En otras palabras,
lo que se plotea en una computadora es realmente un camino poligonal. Sin embargo, podemos “en-
gañar” nuestra vista si ploteamos una cantidad suficiente de puntos. Por ejemplo, grafiquemos la
función seno con una cantidad insuficiente de puntos. Primero construimos los puntos (x, g(x)) alma-
cenando las abscisas y ordenadas en arreglos separados:
>> x = 0:2*pi
x =
0 1 2 3 4 5 6
>> y = sin(x);
Observar que la función seno se aplica componente a componente, y por ende el arreglo y es del
mismo tamaño que x. Para plotear los puntos y unirlos con segmentos de recta usamos el comando
plot
plot(x,y);
shg; % show graph (mostrar gráfica)
El comando shg es usado para que la ventana de ploteo (titulada Figure 1), aparezca encima de
todas las ventanas abiertas (lo cual a veces no es necesario si no hay otras figuras abiertas y no hay
aplicaciones que ocupen toda la pantalla del computador). El resultado de este ploteo se muestra en
la Fig. 2.4. Respondemos el problema con el siguiente script
Figura 2.4: Ploteo de y = sen( x) usando 7 puntos.
Programa 2.9: Script miploteo2d.m

1 clc , close all
2 x = linspace ( -2 ,3 ,80) ;
3 y = sin (5* x ) .* exp ( - x /2) ./(1+ x .^2) ;
4 plot (x , y ) , grid on
5 shg
lo cual nos da la siguiente gráfica
Podemos “embellecer” la gráfica agregando un título (comando title), etiquetas a los ejes (co-
mandos xlabel y ylabel), y eliminando franjas verticales u horizontales en donde la función no
es graficada (comando axis tight). Además, podemos hacer la línea más ancha con el comando
linewidth, así como agregar comentarios para especificar lo que hace el script y cada bloque del
código.
Programa 2.10: Script miploteo2dver2.m

1 % Este programa grafica la funci ó n y = ( sen (5 x ) * e ^{ - x /2}) / (1+ x ^2)
2 % en [ -2 ,3]
3
4 clc ; clear all ;
5
6 % Creaci ó n de puntos de la gr á fica
7 x = linspace ( -2 ,3) ;
8 y = ( sin (5* x ) .* exp ( - x /2) ) ./ (1 + x .^2) ;
9
10 % Ploteo de la funci ó n
11 plot (x ,y , ' linewidth ' ,2) , grid on
12 title ( ' Gr á fica de y = ( sen (5 x ) e ^{ - x /2}) / (1+ x ^2) ') ;
13 xlabel ( 'x '); ylabel ( 'y ') ;
14 axis tight
El resultado de compilar este script es el siguiente:

Gráfica de y = (sen(5x) e -x/2 ) / (1+x 2 )

0.8
0.6
0.4
0.2
0
y
-0.2
-0.4
-0.6
-0.8
-1
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3
x
Si deseamos concentrarnos en un pequeño pedazo de nuestra gráfica (sin cambiar el dominio de

graficación), basta reemplazar la linea 14 por el comando axis con la sintaxis:
axis([xmin xmax ymin ymax]);
donde [xmin,xmax] es el rango de valores para el eje horizontal y [ymin, ymax], el rango de valores
para el eje vertical.
2.2.5.1 Ploteos más elaborados en 2D
Es posible graficar dos o más funciones en una misma ventana
Programa 2.11: Script ploteo3funciones.m

1 % Ploteo de seno , coseno y la exponencial en una misma ventana
2 clc
3 x = linspace (0 ,2* pi ,85) ;
4 y1 = sin ( x ) ;
5 y2 = cos ( x ) ;
6 y3 = exp ( x ) ;
7 plot (x , y1 ,x , y2 ,x , y3 ) , grid on
8 %axis ([0 2* pi -1 10])
9 title ( ' Gr á ficas de 3 funciones en una misma ventana ')
10 xlabel ( 'x ( en cm ) ')
11 ylabel ( 'y ( en cm ) ')
12 shg
En esta gráfica vemos que hay un problema. Las oscilaciones de seno y coseno no se aprecian
en la gráfica. Esto se debe al rango de valores en el eje vertical. Si activamos la línea 8 (borrar el
símbolo %), obtenemos un mejor resultado. Sin embargo, si deseamos que el ancho de línea sea más
grueso, nos damos cuenta que no basta con agregar el atributo 'linewidth',2 a cada par de data
(x, yi ), más bien conseguimos un error. Para aumentar el ancho de línea debemos usar el comando
plot para cada data y los comandos hold on ... hold off, como se muestra a continuación:
Programa 2.12: Script ploteo3funcionesV2.m

1 % Ploteo de seno , coseno y la exponencial en una misma ventana
2 clc
3 x = linspace (0 ,2* pi ,85) ;
4 y1 = sin ( x ); y2 = cos ( x ) ; y3 = exp ( x ) ;
5 plot (x , y1 , ' linewidth ' ,2) , hold on
6 plot (x , y2 , ' linewidth ' ,2)
7 plot (x , y3 , ' linewidth ' ,2) , grid on
8 hold off
9 axis ([0 2* pi -1 10])
10 title ( ' Gr á ficas de 3 funciones en una misma ventana ')
11 xlabel ( 'x ( en cm ) ') , ylabel ( 'y ( en cm ) ')
12 shg
Existe aún un problema, suponiendo que no sabemos cómo lucen las gráficas de estas funciones,
¿cómo sabe una persona que no ha hecho el programa, cuál linea le corresponde a cada función? A
nuestra ayuda viene el comando legend. Basta agregar, después del comando hold off (línea 8), la
siguiente sentencia
legend ( ' sen ( x ) ' , ' cos ( x ) ' ,' exp ( x ) ')
Hasta aquí podemos decir que se entiende qué línea le corresponde a cada función. Pero si imprimi-
mos esta gráfica en una impresora a blanco y negro se presentaría un problema. Podemos resolver
esto si cambiamos el estilo de cada línea que es graficada. La tabla 2.6 nos brinda los comandos dis-
Tabla 2.6: Opciones en M ATLAB al plotear en 2D
ponibles en M ATLAB. La primera columna hace referencia al estilo de línea, la segunda al color de
la línea, y la tercera a los marcadores (markers), los cuales son útiles de cambiar en caso estemos
ploteando data discontinua. Por ejemplo, basta cambiar las líneas 5–7 del programa anterior por
plot (x , y1 , '-k ', ' linewidth ' ,2) , hold on

plot (x , y2 , ' --r ' , ' linewidth ' ,2)
plot (x , y3 , ' -. b ' , ' linewidth ' ,2) , grid on
y así obtener la siguiente figura

Gráficas de 3 funciones en una misma ventana

10
sen(x)
9
cos(x)
exp(x)
8
6
(en cm)
4
y
-1
0 1 2 3 4 5 6
x (en cm)
Bibliografía
[1] D Xue, Y Chen, Scientific Computing with MATLAB, CRC Press, Taylor & Francis Group, Flori-
da, 2da Ed, 2016. Disponible en biblioteca.
45

Clase 02

Cargado por

Copyright:

Formatos disponibles

Clase 02

Cargado por

Información del documento

Título 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

Clase 02

Cargado por

Copyright:

Formatos disponibles

Computación Científica

1 Justificación del Curso 1

2 Fundamentos de Programación en MATLAB y Visualización Científica 9

Soluciones en Computadora a Problemas Matemáticos

Una vez compilado obtenemos

en el Command Window. Es fácil ver que esta expresión significa que

diary ( ' mivar . txt ')

Esto creará el archivo mivar.txt en donde se almacenará completamente el contenido de la varia-

En algunas aplicaciones el cálculo de las raíces de un polinomio es requerido. El teorema de Abel-

134 4 135 3 1215 2 729 729

embargo, en la segunda línea intentamos resolver el problema simbólicamente:

p = [1 9 135/4 135/2 1215/16 729/16 729/64]; roots ( p ) % sol . numerica

La respuesta numérica se muestra a continuación:

¿Recuerda Ud. cómo encontrar el determinante de una matriz n × n?

no invertible). Por ejemplo, consideremos la matriz de Hilbert dada por

>> n = 80; H = hilb ( n ) ; det ( H )

Por otro lado, si cálculo simbólico es usado,

>> H = sym ( hilb ( n ) ) ; d = det ( H )

obtenemos (en aproximadamente 4 segundos)

lo cual es más confiable, y nos indica que el determinante no es exactamente cero.

Considere la ecuación no lineal de Van der Pol

y00 + µ(y2 − 1)y0 + y = 0, y(0) = −1, y0 (0) = 1.

Cuando µ = 1000, muestre cómo resolver numéricamente la ecuación en el intervalo de tiempo

Obtenemos la siguiente gráfica:

Resuelva el siguiente problema de programación lineal

mı́n − 2x1 − x2 − 4x3 − 3x4 − x5 ,

Solución. Como la minimización es un problema de optimización lineal restringida, el método

>> P . f = [ -2 -1 -4 -3 -1]; P. Aineq = [0 2 1 4 2; 3 4 5 -1 -1];

Optimal solution found.

De aquí la solución numérica encontrada se traduce a x1 = 19.7850, x2 = 0, x3 = 3.3200, x4 =

1.1.2 Soluciones analíticas versus soluciones numéricas

exitosamente por décadas. Para simulación de experimentos, soluciones numéricas de ecuaciones en

5. Poderoso sistema de simulación. La poderosa técnica de modelación basada en diagrama de

Elementos Esenciales de MATLAB

2.1.1 Variables y Constantes en MATLAB

y tenemos algunas palabras reservadas usuales muy útiles en programación:

for end while if else elseif fprintf case disp

rado, ver https://www.mathworks.com/help/matlab/matlab_prog/create-a-structure-array.html

2.1.2 Estructuras de Datos

II. Tipo de dato Simbólico

>> syms t x ; assume ( t <= 4) ; assume ( x ~= 2) ;

>> syms x real ; assume (x >= -1) ; assumeAlso (x <5)

Alternativamente, podríamos haber usado una forma más compacta

>> syms x real ; assume (x >= -1 & x <5)

Exhiba los primeros 300 dígitos decimales de π.

>> vpa ( pi ,300)

III. Otros tipos de datos

>> s = ' Hola mundo '

nos da una variable tipo string almacenada en s, de tamaño 1 × 10.

2) Arreglos multidimensionales. Arreglos tridimensionales y de dimensiones mayores son la ex-

>> a = [10 20 30 25 13]

Un vector columna puede ser ingresado de varias maneras

>> b = [2; 4; 6; 8; 10]

o tomando la transpuesta a un vector fila, con el operador ',

>> b = [2 4 6 8 10] '

Una matriz 3 × 3 puede ser asignada como sigue

>> A = [[1 4 7] ' [2 5 8] ' [3 6 9] ']

Your variables are:

o si deseamos más detalle, usamos el comando whos: