DEPARTAMENTO DE

Matem´atica Aplicada UNIVERSIDAD COMPLUTENSE DE MADRID

´ PRACTICAS. ´ CALCULO COMPUTACIONAL EN QU´IMICAS Licenciatura de Ciencias Qu´ımicas Prof.: Rosa Pardo San Gil David Usero

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) Matrices Pr´ actica 1

1.

Introducci´ on

En esta pr´actica vamos a profundizar un poco en las capacidades de Matlab para trabajar con matrices. Veremos en primer lugar algunas operaciones y comandos b´asicos y no tan b´asicos que tiene el programa para trabajar con matrices.

2.

Matrices en Matlab Para introducir una matriz en Matlab se procede de la forma siguiente. Si por ejemplo tenemos la matriz µ ¶ 1 2 3 4 A= 5 6 7 8

se introduce como: >>A=[1 2 3 4; 5 6 7 8] A = 1 5

2 6

3 7

4 8

O bien, >>A=[1,2,3,4;5,6,7,8]; Observemos que unas matrices especiales son los vectores, de esta forma, el vector fila v = (1.0, 1.1,1.2,1.3, . . . , 1.9,2.0), se escribe en Matlab como >>v=[1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]

3.

Operaciones y comandos para Matrices

Hemos visto c´omo se introducen las matrices en Matlab. Veamos un ejemplo para introducir algunos de los comandos b´asicos: Ejemplo 1 Operaciones Elementales Definimos dos matrices: >>A=[2 1;3 2] A = 2 3

1 2

>>B=[3 4;-1 5] B = 3 -1

4 5

• Para sumar las dos matrices: >>A+B ans = 5 2

5 7

• Para multiplicar una matriz por un escalar:

>>3*A ans = 6 9

3 6

• Producto de matrices: >>C=A*B C = 5 7

13 22

Siempre que los tama˜ nos de las matrices sean los adecuados. Para saber cu´ al es el tama˜ no de una matriz con la que estamos trabajando, >>size(A) ans = 2

2

Que quiere decir, evidentemente, 2 filas y 2 columnas. • Para calcular la matriz transpuesta: >>A’ ans = 2 1

3 2

Ejercicio 1 Utilizando las matrices definidas en el ejemplo anterior, comprobar que (AB)t = B t At . (At es la transpuesta de A). Ejemplo 2 Operaciones t´ermino a t´ermino: .* ./ .^ Matlab tiene tres operaciones, que las llamaremos operaciones con punto, que permiten i) multiplicar matrices t´ermino a t´ermino: .* ii) dividir matrices t´ermino a t´ermino: ./ ii) elevar los t´erminos de una matriz a una cierta potencia: .^ Si v es el vector definido en la Secci´ on 2, explorar qu´e hace la orden >>v.^2 Por otra parte, si A y B son las matrices definidas anteriormente, explorar qu´e hacen las ´ ordenes >>A.*B >>A./B Estas operaciones con punto son esenciales en el c´alculo num´erico y se utilizan para representar funciones num´ericamente. Ejemplo 3 Matrices especiales con Matlab • Para generar la matriz identidad cuadrada, >>eye(3) ans = 1 0 0

0 1 0

0 0 1

¿Por qu´e habr´ an elegido el nombre eye? • Una matriz 3 × 2 llena de unos, >>ones(3,2) • Y si queremos que est´e llena de ceros, >>zeros(3,2) • Para generar una matriz con n´ umeros aleatorios uniformemente distribuidos entre 0 y 1, >>rand(3,2) Si se usa el comando randn los n´ umeros aleatorios son normalmente distribuidos, siguiendo la Normal Estandar N (0, 1).

Ejemplo 4 Rango, Inversa y Determinante • Definimos la matriz, >>X=[2 3 4; 1 -1 0] X = 2 3 4 1 -1 0 Para calcular su rango, >>rank(X) ans = 2 • Supongamos que tenemos definida la siguiente matriz, H = 8 3 4

1 5 9

6 7 2

Para calcular su inversa, >>inv(H) ans = 0.1472 -0.0611 -0.0194

-0.1444 0.0222 0.1889

0.0639 0.1056 -0.1028

Y si queremos ver el resultado en forma racional, >>format rational >>inv(H) ans = 53/360 -13/90 -11/180 1/45 -7/360 17/90

23/360 19/180 -37/360

(Para ver todas las opciones del comando format hacer help format) • Para calcular el determinante de la matriz anterior H, >>det(H) ans = -360 Ejercicio 2 Generar una matriz cualquiera, por ejemplo 25×25, y calcular su inversa, su rango y su determinante. (¡No imprimirla!) ¿Qu´e ocurre con el determinante de la matriz y el de su inversa? Ejemplo 5 Los comandos especiales rref y rrefmovie • El comando rref produce la forma reducida escalonada por filas de una matriz usando la eliminaci´ on de Gauss-Jordan, es decir, haciendo ceros por debajo y por encima de la diagonal principal sin mover las columnas. Por ejemplo, definimos la matriz, >>A=[-1 2 -1;2 1 2;2 4 2] A = -1 2 2

2 1 4

-1 2 2

Ahora escribimos el comando aplicado a la matriz,

>>R=rref(A) R = 1 0 0

0 1 0

1 0 0

Ejercicio 3 a) Calcular el rango de la matriz siguiente utilizando el comando rref   16 2 3 13  5 11 10 8   A=  9 7 6 12 4 14 15 1 b) Si una matriz H es cuadrada y no singular, es decir det(H) 6= 0, ¿cu´ al ser´ a la matriz R = rref(H)? c) ¿C´ omo podemos utilizar estos comandos para calcular la inversa de una matriz invertible? Aplicarlo a la matriz,   8 1 6 B = 3 5 7 4 9 2 Para verificar el resultado se puede calcular la inversa directamente con inv(B).

4.

Matrices “dispersas”

Ejemplo 6 A veces usamos matrices con muchos ceros. MatLab tiene una forma de trabajar con ellas usando menos bytes con el comando sparse. Ve´ amoslo con un ejemplo: Introducimos una matriz: >>A=[0 0 0 3;0 0 -1 2;3 0 0 1;0 0 0 -2]; - Para convertirla a matriz dispersa >>s=sparse(A) Si preguntamos ‘‘whos" vemos que s ocupa menos que A. - Para recuperar la matriz inicial >>full(s) Se pueden generar directamente matrices “sparse”: >>sparse(i,j,s,m,n) donde: i,j son los sub´ındices de los elementos no nulos (i,j son vectores) s es un vector con los valores de los elementos no nulos (m,n) es el tama˜ no de la matriz. De modo que, en el ejemplo anterior, para generar s deber´ıamos escribir: >>i=[1 2 2 3 3 4]; >>j=[4 3 4 1 4 4]; >>s=[3 -1 2 3 1 -2]; >>m=4;n=4; >>sparse(i,j,s,m,n)# Y obtenemos s. Para obtener la matriz inicial >>full(s) Ejercicio 4 Utilizando el comando sparse, generar la  0 1 0 −1 0 1   0 −1 0   ..  .  0 0 0 0 0 0 (Visualizarla para comprobar que est´a bien.)

matriz 20 × 20 0 0 1

... ... ... .. .

0 0

... ...

 0 0  0  ..  .  0 1 −1 0 0 0 0

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) Sistemas de ecuaciones lineales Pr´ actica 2 En esta pr´actica vamos a ver c´omo se pueden resolver sistemas de ecuaciones lineales utilizando Matlab.

1.

Sistemas de ecuaciones lineales Un sistema de ecuaciones lineales,   a11 x1 + a12 x2 + . . . + a1n xn = b1   a21 x1 + a22 x2 + . . . + a2n xn = b2 ..   .    am1 x1 + am2 x2 + . . . + amn xn = bm con m ecuaciones y n inc´ognitas se puede escribir en forma matricial, Ax = b donde, 

a11  a21  A= .  ..

a12 a22

... ...

am1

am2

...

 a1n a2n   ;  amn

 x1  x2    x= .   ..  

 b1  b2    b= .   ..  

y

xn

bm

Vamos a ver mediante algunos ejemplos y ejercicios c´omo se pueden resolver los sistemas de ecuaciones lineales utilizando algunos de los comandos de Matlab descritos anteriormente. Ejemplo 1 Consideremos el sistema,   2x − y + z = 3 x+y =3   y − 3z = −7 entonces, siguiendo la notaci´ on anterior,   2 −1 1 0 , A = 1 1 0 1 −3

  x x = y  z



y

 3 b= 3  −7

Como se trata de un sistema con soluci´ on u ´nica, ya que el determinante de A es distinto de cero, >>det(A) ans = -8 Una forma de resolver el sistema es escribir la matriz orlada (o ampliada) >>Ab=[A b] y hacer rref(Ab) con lo que obtenemos

 1 0 0

0 1 0

0 0 1

 1 2 3

es decir, la soluci´ on es x = 1, y = 2, z = 3. Otra forma de resolver el sistema consiste en despejar x, x = A−1 b, sin m´ as que escribir

>>x=inv(A)*b x = 1 2 3 Hay otra forma de hacerlo, utilizando lo que en Matlab se denomina como divisi´ on matricial a la izquierda: >>x=A\b x = 1 2 3 En este caso, el resultado es el mismo, pero es diferente la forma en la que trabaja el ordenador. En este segundo caso el m´etodo que utiliza es el de la factorizaci´ on LU, que es una modificaci´ on de la eliminaci´ on gaussiana. Vamos a ver c´ omo resuelve Matlab, internamente, el sistema cuando se utiliza la opci´ on: >>x=A\b. El proceso se puede dividir en tres etapas: 1) Calcula una matriz triangular inferior L, una matriz triangular superior U y una matriz de permutaci´ on P tales que P A = LU . P es simplemente la matriz identidad I con sus filas cambiadas de orden. 2) Resuelve Ly = P b. 3) Por u ´ltimo, se resuelve U x = y. La primera etapa es lo que se conoce con el nombre de factorizaci´ on LU y es el paso m´ as importante. Por lo tanto, en Matlab ser´ıa equivalente utilizar: >>x=A\b que utilizar el siguiente proceso: >>[L,U,P]=lu(A); % Este es el comando que calcula las matrices L, U, P >>B=P*b; >>y=L\B; >>x=U\y Ejercicio 1 Resolver el siguiente sistema utilizando los tres procedimientos anteriormente descritos y comprobar que se obtiene la misma soluci´ on.      1 1 0 3 x1 4 2  x2   1  1 −1 1       3 −1 −1 2  x3  = −3 −1 2 3 −1 x4 4 Ejercicio 2 El m´etodo >>x=A\b funciona si hay distinto n´ umero de ecuaciones que inc´ ognitas utilizando el m´etodo de los m´ınimos cuadrados. Sea A una matriz m × n, el m´etodo busca x que minimiza   2 X X X X   2 f (x) = kAx − bk2 = ((Ax)i − bi ) = aij xj  − 2bi aij xj + b2i  ,  i

i

j

para ello calculamos sus derivadas parciales   X X X ∂f  =2 aij xj  aik − 2 bi aik , ∂xk i j i y

∂f ∂xk

= 0,

j

k = 1, · · · , n

∀k = 1, · · · , n si y s´ olo si X ij

aij aik xj =

X

bi aik ,

∀k = 1, · · · , n

i

es decir >>x=A\b resuelve el sistema lineal AT Ax = AT b, que ahora tiene el mismo n´ umero de ecuaciones n que de inc´ ognitas n y minimiza la distancia kAx − bk.

Cuando el rango de A es k < n, o hay menos ecuaciones que inc´ ognitas, m < n, tenemos un sistema con infinitas soluciones (o ninguna), Matlab nos da una soluci´ on que tiene ceros para algunas de las inc´ ognitas, al menos n − k. Aplicando lo anterior encontrar las soluciones de los sistemas  (  x + y = 1 x+y+z =3 a) x − y = 1 b)  x−y+z =2  2x + y = 2 Ejercicio 3 Con otro comando, >>x=pinv(A)*b se obtiene una soluci´ on del sistema por el m´etodo de los m´ınimos cuadrados tal que la norma eucl´ıdea de x es la m´ as peque˜ na de todas las posibles. ¿Es siempre u ´nica esta soluci´ on con longitud m´ınima? Investigar qu´e hace el comando >>pinv(A) para los ejemplos anteriores Ejercicio 4 Si resolvemos el siguiente sistema con Matlab  (  x + y = 1 x + 2y + 3z = 1 2x + 2y = 2  x + 2y + 3z = 5  3x + 3y = 3 1. Utiliza los m´etodos A\b y pinv(A)*b para resolver los sistemas 2. Resuelve el sistema AT Ax = AT b con los dos m´etodos anteriores 3. Utiliza el comando rref para resolver los dos sistemas anteriores 4. ¿Qu´e hemos obtenido con los comandos A\b y pinv(A)*b? Ejercicio 5 Estudiar el sistema

 x+y+z =1    x − y − z = 2  3x + 3y + 3z = 1    −x − y + z = 1

¿Qu´e obtenemos con los comandos A\b y pinv(A)*b? ¿Es correcto? ¿Que responde Matlab al comando rref([A,b])? ¿Y a los comandos x=inv(A)*b, lu(A)? ¿Qu´e conclusi´ on sacas de esto? Ejercicio 6 Crear las matrices A=rand(200) y b=rand(200,1). Vamos a resolver el sistema Ax = b con MatLab de varias formas, y vamos a contar el tiempo de CPU que se utiliza para resolver el sistema utilizando cuatro m´etodos. Para ello vamos a utilizar el comando cputime. Primero poner el contador a cero. Segundo resolver el sistema de una de las formas y tercero contar el tiempo, escribiendo >>t=cputime; rref([A b]); trref=cputime-t y repetir la operaci´ on resolvi´endolo de otra forma. >>t=cputime; metodo; tmet=cputime-t ¿En cu´ al de los tres casos se utilizan menos tiempo de cpu? i) ii)

>>t=cputime; x=rref([A b]); trref=cputime-t >>t=cputime; x=pinv(A)*b; tpinv=cputime-t

iii)

>>t=cputime; x=inv(A)*b; tinv=cputime-t

iv)

>>t=cputime; x=A\b; t=cputime-t

Comprobar que t>poly2sym(p,’x’) Para evaluarlo en un punto: >>polyval(p,4) Tambi´en es posible evaluar el polinomio en una matriz: >>A=[4 4;1 0]; >>polyval(p,A) >>polyvalm(p,A) Ejercicio 1 Hemos visto que un polinomio puede evaluarse en una matriz utilizando dos comandos: polyval y polyvalm. ¿Qu´e hace exactamente cada uno de estos comandos? ¿Se pueden aplicar a cualquier matriz sea o no cuadrada? ¿Por qu´e?

2.

Operaciones

Ejemplo 1 Sean los polinomios P (x) = 4x5 − 8x2 − 10x + 2 y Q(x) = x2 − x + 1 >>p=[4 0 0 -8 -10 2];q=[1 -1 1]; Para multiplicarlos: >>m=conv(p,q),poly2sym(m,’x’) Para dividir P (x) entre Q(x): >>[c,r]=deconv(p,q) >>cociente=poly2sym(c,’x’) >>resto=poly2sym(r,’x’) Ejercicio 2 ¿Qu´e habr´ a que hacer para sumar dos polinomios? Sumar los polinomios P (x) = x3 − 3x2 − 5x + 2 y 5 2 Q(x) = 4x − 3x + 8x + 2. Sugerencia: Sea >>p=[1 -3 -5 2]; reescribe p para igualar el tama˜ no de los vectores p q y as´ı poder sumarlos >>p=[zeros(...) p]; Ejercicio 3 Dividir el polinomio P (x) = x5 − 8x4 − 10x3 + 10x2 − x + 6 entre Q(x) = x + 3. Comprobar que Dividendo = Divisor·Cociente + Resto Verificar que se cumple el teorema del resto, es decir, que el resto de dividir P (x) entre (x − a) es igual a P (a).

3.

Ra´ıces Para calcular las ra´ıces del polinomio P (x) = x5 − 2x4 − 10x3 + 20x2 + 9x − 18

introducimos el polinomio,

>>p=[1 -2 -10 20 9 -18]; y escribimos: >>r=roots(p) Tambi´en es posible reconstruir el polinomio a partir de sus ra´ıces: >>p=poly(r);poly2sym(p,’x’) Ejemplo 2 Ra´ıces en´ esimas de la unidad. Se llaman ra´ıces en´esimas de la unidad a las soluciones (complejas) de la ecuaci´ on zn = 1 Para calcular las ra´ıces quintas de la unidad, tenemos que calcular las ra´ıces del polinomio P (x) = x5 − 1. Introducimos el polinomio: >>p=[1 0 0 0 0 -1] Calculamos las ra´ıces: r=roots(p) Estas ra´ıces se pueden representar utilizando el comando >>compass(r) Ejercicio 4 Calcular y dibujar las ra´ıces sextas de la unidad. Ejercicio 5 Escribir un programa de tipo function que, tomando como dato de entrada el valor de n, calcule y dibuje las ra´ıces en´esimas de la unidad. Con el programa anterior, calcular y dibujar el caso n = 8.

4.

Gr´ aficas de polinomios

Ejemplo 3 Dibujar la gr´ afica del polinomio

y = x3 − x

en el intervalo [−2, 2]. Introducimos el polinomio: >>p=[1 0 -1 0]; Generamos una tabla de valores para x: >>x=linspace(-2,2,1000); sustituimos en el polinomio: >>y=polyval(p,x); y representamos los puntos con el comando >>plot(x,y) Ejercicio 6 Escribir un programa de tipo function que , tomando como dato de entrada el valor de n, dibuje la gr´ afica de la curva y = xn en el intervalo [−1, 1]. Con el programa anterior, dibujar y = x10 y el caso y = x13 .

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) Gr´ aficas 2D. Pr´ actica 4

1.

Introducci´ on

Con el programa MATLAB podemos dibujar gr´aficas de curvas y funciones en el plano en m´ ultiples formatos y con diferentes presentaciones. Funciones en coordenadas cartesianas, dadas en forma expl´ıcita, es decir, de la forma y = f (x); sin embargo, no se pueden dibujar directamente curvas en forma impl´ıcita, es decir, de la forma g(x, y) = 0. Tambi´en curvas en forma param´etrica, es decir, de la forma ~r(t) = (x(t), y(t)) con a ≤ t ≤ b. Y tambi´en curvas en coordenadas polares, de la forma r = r(θ) con θ1 ≤ θ ≤ θ2 .

2.

Comandos b´ asicos para gr´ aficos 2D

Para dibujar una gr´afica 2D con MATLAB se puede utilizar el comando ezplot, por ejemplo, >>ezplot(’sin(x)) Tambi´en se puede generar una tabla de valores para la x y para la y de la funci´on a dibujar, por ejemplo, >>x=linspace(0,2*pi,30); >>y=sin(x); A continuaci´on utilizar un comando para dibujar, que puede ser, >>plot(x,y) El color y el estilo de las l´ıneas que se utilizan para hacer las gr´aficas se pueden modificar, por ejemplo, con el comando >>plot(x,y,’r:’) obtenemos la gr´afica en color rojo y punteada, en lugar de con l´ınea continua. Para ver los colores y estilos disponibles consultar el comando plot: >>help plot

3.

Coordenadas cartesianas

Ejemplo 1 Dibujar la gr´ afica de la funci´ on  2  x 0 y=  p

3(x − 1)

si x < 0 si 0 < x < 1 si x > 1

Creamos una tabla de valores: >>x=linspace(-3,3,1000); >>y=x.^2.*(x0)&(x1); Y ahora utilizamos alguno de los comandos de dibujo, por ejemplo, >>plot(x,y,’m’) que producir´ a una gr´ afica en color magenta. Ejercicio 1 Dibujar las gr´ aficas de las siguientes funciones eligiendo, en cada caso, una tabla de valores adecuada para que aparezcan los aspectos m´ as representativos de la funci´ on: ( x2 si x < 0 a) f (x) = −1 si x ≥ 0   −x si x < −1 b) f (x) = 1 si 0 < x < 2   2 −x si x > 2 √   1 − x si x < −1 c) f (x) = 1 − x2 si −1 < x < 1  √ x − 1 si x > 1

4.

Ecuaciones param´ etricas

Ejemplo 2 Dibujar la gr´ afica de la curva ~r(t) = (cos(t), sen(t)) ;

−π ≤ t ≤ π

En primer lugar generamos los valores de t en el intervalo indicado, t=linspace(-pi,pi,100); Y ahora lo podemos dibujar: >>plot(cos(t),sin(t)) Ejercicio 2 Dibujar las curvas en param´etricas siguientes, a)~r(t) = (2 cos3 t, 2 sen3 t); b)~r(t) = (3 sen t, 2 sen(2t)); −π ≤ t ≤ π µ µ ¶ µ ¶¶ t t t t c)~r(t) = 12( )2 − 9 , (( )2 − 1)16( )2 + 2 ; π π π π µ ¶ 3 d)~r(t) = cos t(cos t + 1), 2 sen(2t) ; −π ≤ t ≤ π 2 e)~r(t) = (sen(2t) + sen t, − cos(2t) − cos t);

5.

−π ≤ t ≤ π

−3 ≤ t ≤ 3

−π ≤ t ≤ π

Coordenadas polares

Ejemplo 3 Dibujar la gr´ afica de r = 2 − 4 cos(θ),

−π ≤ θ ≤ π

En primer lugar generamos los valores del ´ angulo tetha: >>tetha=linspace(-pi,pi,100); Calculamos los valores de r: >>r=2-4*cos(tetha); Y dibujamos la gr´ afica, >>polar(tetha,r) Ejercicio 3 Dibujar las gr´ aficas de las siguientes funciones, dadas en coordenadas polares: a)r = 7 − 7 sen(θ); −π ≤ θ ≤ π b)r = 3 − 6 sin(θ);

−π ≤ θ ≤ π

c)r = sen(6θ);

−π ≤ θ ≤ π

d)r = cos(8θ);

−π ≤ θ ≤ π

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) Gr´ aficas 3D Pr´ actica 5 Con MATLAB se pueden hacer gr´aficas de funciones de dos variables en el espacio, es decir, funciones de la forma z = f (x, y) (forma expl´ıcita). Tambi´en se pueden dibujar, en el espacio, curvas y ciertas superficies.

1.

Curvas en el espacio

Se generan de una manera similar a las curvas en el plano, con la diferencia de que aqu´ı se utiliza el comando plot3 Ejemplo 1 Dibujar la curva ~r(t) = (sen(t), cos(t), t)

0 ≤ t ≤ 4π

y sobre ella los vectores velocidad. Generamos los valores de t: >>t=linspace(0,4*pi,500); Y ahora podemos utilizar el comando plot3 que nos da el dibujo completo: >>plot3(sin(t),cos(t),t) Ejercicio 1 Representar las curvas siguientes en los intervalos indicados: a)~r(t) = (2 cos3 t, 2 sen3 t, t) − 4 ≤ t ≤ 3. t t t 4 4 d)~r(t) = (e sen(2t), e cos(2t), ) − 10 ≤ t ≤ 4,8. 4

2.

Gr´ aficos de funciones z = f (x, y)

Para dibujar gr´aficos de funciones de dos variables z = f (x, y), al igual que para funciones de una variable, en primer lugar hay que generar tablas de valores para las variables x e y. En realidad ahora lo que tenemos que hacer es generar un mallado sobre un rect´angulo del plano XY . Para eso se utiliza el comando meshgrid. afica de la funci´ on Ejemplo 2 Dibuja la gr´ 2

z = e−x

−y 2

en la regi´ on del plano D = {(x, y)/ − 2 ≤ x ≤ 2, − 3 ≤ y ≤ 3} Habr´ a que efectuar los pasos siguientes: Generamos el mallado >>[x,y]=meshgrid(-2:.2:2,-3:.2:3); Sustituimos en la funci´ on para calcular los valores de z: >>z=exp(-x.^2-y.^2); Y ahora podemos dibujar el gr´ afico con alguno de los siguientes comandos: >>plot3(x,y,z) >>mesh(x,y,z) >>surf(x,y,z) Para conseguir efectos de sombreados y colores diferentes se pueden consultar todas las posibilidades de los comandos colormap y shading. Algo que resulta tambi´en interesante es a˜ nadir una escala de colores al dibujo afica; esto se consigue con el que nos permite conocer las alturas (coordenada z) de los diferentes puntos de la gr´ comando colorbar (despu´es de dibujada la gr´ afica). Las longitudes de los ejes coordenados tambi´en se pueden modificar con el comando: >>axis([xmin xmax ymin ymax zmin zmax]) Ejemplo 3 Representar la gr´ afica de la funci´ on z = x2 + y 2 Dibujando algunas curvas de nivel. Creamos el mallado: >>[x,y]=meshgrid(-2:.1:2);

Sustitu´ımos en la funci´ on, para calcular los valores de z: >>z=x.^2+y.^2; Ahora podemos dibujar la gr´ afica utilizando alguna de los comandos descritos en el punto relativo a los comandos b´ asicos de 3D. Las curvas de nivel se pueden hacer utilizando alguno de los comandos siguientes: >>contour(x,y,z,10) % dibuja 10 curvas de nivel >>contour(x,y,z,1.38) % dibuja la curva de nivel z=1.38 >>contour3(x,y,z,10) % dibuja 10 curvas de nivel en el espacio Ejercicio 2 Escribir un fichero tipo script con nombre missuperficies.m que represente las gr´ aficas de las siguientes funciones de 2 variables, utilizando alguno de los comandos descritos anteriormente. p a)z = − |xy| 2 2 b)z = e−(x +y ) Ejercicio 3 Escribir un fichero tipo script con nombre miscurvasdenivel2D.m que dibuje las curvas de nivel de las funciones del ejercicio anterior, utilizando alguno de los comandos descritos anteriormente.

3.

Superficies pVT

Las ecuaciones de estado de los gases son, desde un punto de vista matem´atico, ecuaciones tipo z = f (x, y) en las que las variables cartesianas (x, y, z) son sustituidas por las variables termodin´amicas temperatura (T ), volumen molar (v = V /n) y presi´on (p). Ejemplo 4 La ecuaci´ on de los gases ideales

RT v puede dibujarse de manera an´ aloga a como lo hemos hecho con otras funciones anteriormente. Primero deberemos generar un mallado: >>[T,v]=meshgrid(0:4:200,2:2:20); A continuaci´ on sustituimos los valores en la funci´ on para calcular los valores de la presi´ on: >>p=0.082*T./v; Y por u ´ltimo pintamos la gr´ afica resultante con el comando mesh o con cualquier otro: >>mesh(v,T,p) p=

3.1.

Isotermas del Gas Ideal

En muchas ocasiones las superficies p v T del gas no son las m´as interesantes debido a la dificultad de ver el comportamiento seg´ un las variables. Por este motivo se suelen utilizar diagramas en dos dimensiones en los que se muestran las variaciones de la presi´on (eje Y ) respecto al volumen molar (eje X) para valores constantes de la temperatura. A estas curvas se las denomina isotermas. Obtenidos los valores de las variables termodin´amicas p, v y T con la ecuaci´on de estado resulta muy sencillo en MATLAB dibujar las isotermas. Basta recordar que el comando contour(x,y,z,N) dibujaba N curvas en el plano (x, y) para valores de z constantes. As´ı se procede de la siguiente manera: >>contour(v,p,T,10) El resultado ser´a una gr´afica con 10 isotermas en las que se observa la variaci´on de la presi´on (eje Y ) respecto a la variaci´on del volumen molar (eje X) para valores constantes de la temperatura. Variando los par´ametros del mallado podremos cambiar la regi´on de la gr´afica dibujada, siempre manteni´endolos dentro de los l´ımites admisibles. (Nunca usar temperaturas negativas, etc.) Ejercicio 4 Cambia el orden en el que se escriben los argumentos en el comando contour() . ¿Que obtienes si escribes contour(v,T,p,10)? ¿Y que obtienes si escribes contour(T,p,v,10)? ¿Cu´ al es su significado termodin´ amico?

3.2.

Otras ecuaciones de Estado

La ecuaci´on de estado de los gases ideales, si bien resulta muy f´acil de manejar, no es una ecuaci´on v´alida para ning´ un gas. Por lo general las mol´eculas de cualquier gas siempre van a sufrir interacci´on entre ellas, y esto originar´a que su comportamiento se desv´ıe del correspondiente a un gas ideal. Para obtener comportamientos m´as parecidos al de los gases reales se han propuesto numerosas ecuaciones de estado que en su formulaci´on incluyen, en cierta medida, este tipo de interacciones.

La m´as conocida fue introducida por J. D. Van der Waals en 1873: p=

RT a − , v − b v2

donde a y b son dos constantes positivas cuyo valor se puede ajustar y dependen del gas. Los mejores valores se obtienen tomando: 27R2 Tc2 a= 64pc b=

RTc 8pc

donde Tc y pc son los valores de la temperatura y la presi´on en el punto cr´ıtico. El par´ametro b se denomina par´ametro de repulsi´on y representa un volumen m´ınimo del gas, por debajo del cual resulta imposible comprimirlo debido al tama˜ no de sus mol´eculas, mientras que el par´ametro a se denomina par´ametro de atracci´on y nos da una correcci´on de la ecuaci´on del gas ideal debido a la atracci´on entre sus mol´eculas. Otra ecuaci´on muy utilizada fue introducida en 1949 por O. Redlich y N. S. Kwong p=

RT a − v − b v(v + b)T 1/2

donde de nuevo a y b son constantes cuyo valor depende del gas y cuyo mejor valor es: 5/2

a=

0,42748R2 Tc pc

b=

0,08664RTc pc

siendo Tc y pc los valores de la temperatura y la presi´on en el punto cr´ıtico. La principal diferencia con la ecuaci´on anterior est´a en que el t´ermino que incluye la interacci´on con las dem´as part´ıculas del gas tiene adem´as una dependencia con la temperatura. Ejercicio 5 En la tabla siguiente se muestran los valores cr´ıticos de la temperatura, presi´ on y volumen molar para diferentes sustancias. Sustancia Dioxido de Carbono (CO2 ) Ox´ıgeno (O2 ) Agua (H2 O) Helio (He)

Tc (K) 304.2 154.4 647.4 5.2

pc (atm) 72.9 49.7 218.3 2.26

vc (l/mol) 0.094 0.074 0.056 0.060

Crea un programa function que realice los sigientes c´ alculos: para cada una de las sustancias anteriores, obt´en los par´ ametros a y b de la ecuaci´ on de Van der Waals y dibuja una gr´ afica con 8 isotermas. Compara las isotermas obtenidas con las correspondientes a la ecuaci´ on de los gases ideales, ¿qu´e diferencias hay?. Ejercicio 6 Repite el ejercicio anterior para la ecuaci´ on de Redlich-Kwong. Ejercicio 7 Tomando como valores los de vc y Tc de la tabla anterior se puede evaluar la presi´ on en cada una de las tres ecuaciones de estado anteriores. De esta manera es posible obtener una estimaci´ on te´ orica para pc en cada ecuaci´ on. Haz los c´ alculos para el CO2 , el O2 , el H2 O y el He. Rellena la tabla siguiente compar´ andolos con el experimental. ¿Porqu´e no se obtiene el valor de pc experimental? (a pesar de que ha sido usado para calcular a y b). ¿Qu´e ecuaci´ on estima mejor cada uno de los valores experimentales?. ¿Cu´ al de las ecuaciones lo estima peor?. ¿Sabes decir porqu´e?. Sustancia CO2 O2 H2 O He

pc Exp. 72.9 49.7 218.3 2.26

pc Gas Id.

pc V-W

pc R-K

4.

Orbitales del ´ atomo de Hidr´ ogeno

El movimiento de un electr´on alrededor de un n´ ucleo de Hidr´ogeno viene descrito, siguiendo un planteamiento cu´antico, por unas funciones de onda Ψn,l,m (x, y, z) cuya norma al cuadrado describe la probabilidad de encontrar dicho electr´on en un punto del espacio (x, y, z) y dentro de un determinado orbital determinado por los n´ umeros cu´anticos n, l, m. Al n´ umero cu´antico n se le llama n´ umero cu´ antico principal y determina el nivel de energ´ıa correspondiente al electr´on En = −1/(2n2 ). Conviene recordar que a diferencia de la descripci´on cl´asica, un electr´on no puede tener cualquier energ´ıa dentro de ese potencial sino solo algunos valores puntuales que forman el espectro de energ´ıa. Los distintos valores de n = 1, 2, 3, ... determinan la capa en la que se encuentra el electr´on, lo que a su vez estar´a relacionado con un sinf´ın de propiedades qu´ımicas del elemento. El n´ umero cu´antico l se le conoce como n´ umero cu´ antico azimutal y determina en momento angular del electr´on. umero azimutal Los valores permitidos para el n´ umero l van desde l = 0 hasta l = n − 1. Los distintos valores del n´ van a determinar los diferentes orbitales existentes dentro de cada capa. As´ı l = 0 se corresponde con el orbital s, l = 1 con el orbital p, l = 2 con el d, etc. El n´ umero m llamado n´ umero cu´ antico magn´etico determina los distintos estados para un valor de l. Determinan, en cierto modo, la orientaci´on y simetr´ıa del orbital en cuesti´on. Sus valores van desde m = −l hasta m = l. Los diferentes valores de m originan los diferentes suborbitales. As´ı para l = 1, (orbital p) tenemos que m = −1 corresponde al orbital py , m = 0 al orbital pz y m = 1 al orbital px . Para l = 2 tenemos los suborbitales dxy (m = −2), dyz (m = −1), dz2 (m = 0), dxz (m = 1), d(x2 −y2 ) (m = 2). Por u ´ltimo existe un cuarto n´ umero cu´antico, que no se tiene en cuenta en la descripci´on de los orbitales, conocido como n´ umero cu´ antico de spin cuyos valores pueden ser de s = ±1/2. Este n´ umero cu´antico hace posible que en cada orbital coexistan dos electrones que s´olo se diferencian en este n´ umero cu´antico de spin.

4.1.

Funciones de Onda

Matem´aticamente hablando, las funciones de onda descritas anteriormente Ψ(x, y, z) son las soluciones de una ecuaci´on diferencial en derivadas parciales conocida como ecuaci´on de Schr¨odinger estacionaria −

¯ 2 h ∇ Ψ(x, y, z) + V (x, y, z)Ψ(x, y, z) = EΨ(x, y, z) 2m

donde V (x, y, z) es el potencial de Coulomb correspondiente a la carga del n´ ucleo, E es la energ´ıa del sistema y ∇ es el operador gradiente y por tanto ∂2 ∂2 ∂ ∇2 = + + 2. ∂x2 ∂y 2 ∂z Las soluciones de esta ecuaci´on se obtienen realizando una transformaci´on a coordenadas polares y posteriormente una separaci´on de variables de manera que para un determinado nivel de energ´ıa En la funci´on de onda se puede descomponer de la forma siguiente Ψn,l,m = Rn,l (r)Yl,m (θ, φ), donde Rn,l es una componente radial de la funci´on y Yl,m es una componente direccional. Escritas en coordenadas cartesianas, las funciones de onda correspondientes a los tres primeros niveles de energ´ıa son (en unidades de a0 ): ³ e−r r ´ e−r/2 (2r2 − 18r + 27) e−r/3 √ √ Ψ1,0,0 = √ Ψ2,0,0 = 1 − Ψ3,0,0 = 2 81 π 8π 3π y e−r/2 z e−r/2 x e−r/2 √ √ √ Ψ2,1,−1 = Ψ2,1,0 = Ψ2,1,1 = 2 2 2 8π ³ 8π ³ 8π ³ ´ ´ r r r ´ −r/3 8 8 8 −r/3 −r/3 y 1− e z 1− e x 1− e Ψ3,1,−1 = √ Ψ3,1,0 = √ Ψ3,1,1 = √ 6 6 6 27 8π 27 8π 27 8π 2 2 1 x y e−r/3 y z e−r/3 Ψ3,2,−2 = √ Ψ3,2,−1 = √ Ψ3,2,0 = √ (3z 2 − r2 ) e−r/3 81 2π 81 2π 81 3π 2 1 x z e−r/3 Ψ3,2,1 = √ Ψ3,2,2 = √ (x2 − y 2 ) e−r/3 81 2π 81 2π p donde r = x2 + y 2 + z 2 es el radio. Las componentes angulares se ponen de manifiesto en la dependencia en cada una de las variables cartesianas de las funciones de onda descritas.

Ejemplo 5 Curvas de nivel del orbital 2px . Es imposible representar gr´ aficamente los orbitales ya que al tratarse de funciones de IR3 sobre IR necesitar´ıamos un espacio de 4 dimensiones para representarlo. Una de las posibilidades para hacernos una idea de los orbitales resulta el dibujar curvas de nivel de los mismos. En el presente ejemplo vamos a representar las curvas de nivel del orbital 2px en el corte con el plano z = 0. Lo mismo se puede intentar para otros valores de z, x o y. Para ello necesitaremos en primer lugar definir un fichero con la correspondiente funci´ on de onda, escribiendo en un fichero: function Psi=FdOn211(x,y,z) % % Funcion de Onda n=2 l=1 m=1 % r=sqrt(x.^2+y.^2+z.^2); Psi=x/2.*exp(-r./2)/sqrt(8*pi);

%

orbital 2px

y lo guardamos con el nombre FdOn211.m. Esto nos permitir´ a evaluar la funci´ on de onda para cada punto (x, y, z) del espacio. A continuaci´ on escribe un fichero de nombre miscurvasdenivel3D.m que 1. defina un mallado de valores de x e y, para valores comprendidos entre −20a0 y 20a0 : >>[x,y]=meshgrid(-20:.1:20); 2. evalue la funci´ on Ψ(x, y, 0) para obtener los puntos sobre los que dibujar las curvas de nivel: >>Psi=FdOn211(x,y,0); 3. dibuje las curvas de nivel con el comando contour: >>contour(x,y,Psi) 4. Si queremos dibujar m´ as curvas de nivel basta a˜ nadir el n´ umero de curvas deseadas >>contour(x,y,Psi,16)

Ejercicio 8 Siguiendo el proceso descrito anteriormente dibuja las curvas de nivel de los orbitales 1. Ψ3,1,1 → Orbital 3px 2. Ψ3,2,−2 → Orbital 3dxy 3. Ψ3,2,2 → Orbital 3d(x2 −y2 ) Ejemplo 6 Escribe en un fichero script de nombre unasuperficiesdenivel.m que dibuje una superficie de nivel para el orbital at´ omico Ψ3,2,0 .

Figura 1: Orbital 3dz2

[x,y,z]=meshgrid(-20:0.5:20); Psi=FdOn320(x,y,z); m=min(Psi(:)); M=max(Psi(:)); clf %clear figure isovalue=(M+m)/2; %valores de Psi=isovalue fv = isosurface(x,y,z,Psi,isovalue); % La superficie de nivel hpatch = patch(fv); isonormals(x,y,z,Psi,hpatch); Alphalevel=0.5; set(hpatch,’FaceAlpha’,Alphalevel,’FaceColor’,[0 0 1],’EdgeColor’,’none’) daspect([1,1,1]) view(3); axis tight camlight left; Ejemplo 7 Escribe en un fichero script de nombre missuperficiesdenivel.m que dibuje una serie de superficies de nivel para el orbital at´ omico Ψ3,2,0 .

Figura 2: Orbitales 3dz2

[x,y,z]=meshgrid(-20:0.5:20); Psi=FdOn320(x,y,z); m=min(Psi(:)); M=max(Psi(:)); n=5 color=[1 1 0;1 0 1; 0 1 1;0 1 0;1 0 0;0 0 1; 0 0 0]; %colores de las superficies clf %clear figure for i=1:n isovalue=m+i*(M-m)/n; %valores de Psi=isovalue fv = isosurface(x,y,z,Psi,isovalue); % La superficie de nivel hpatch = patch(fv); isonormals(x,y,z,Psi,hpatch); Alphalevel=0.5; set(hpatch,’FaceAlpha’,Alphalevel,’FaceColor’,color(i,:),’EdgeColor’,’none’) daspect([1,1,1]) view(3); axis tight

camlight left; end Ejercicio 9 Repite el ejemplo anterior para las funciones de onda citadas en el ejercicio anterior. Ejercicio 10 Seg´ un la teor´ıa cu´ antica, la determinaci´ on seg´ un los ejes cartesianos de los tres orbitales p no es m´ as que una de las muchas soluciones posibles, ya que la colocaci´ on de los ejes es arbitraria. Por lo tanto los tres orbitales p juntos deben tener momento angular total cero, y el aspecto del orbital debe ser esf´erico. Comprueba esta afirmaci´ on pintando en la misma figura los tres orbitales en conjunto.

4.2.

Hibridaci´ on

Seg´ un la teor´ıa elemental de ecuaciones diferenciales, si una ecuaci´on diferencial tiene dos soluciones conocidas x1 y x2 , entonces ser´a soluci´on cualquier combinaci´on lineal de ambas. Esto tambi´en es extensible al caso que nos ocupa, de manera que dadas un conjunto de orbitales soluciones de la ecuaci´on de Schr¨odinger, cualquier combinaci´on lineal de ´estas ser´a un nuevo orbital soluci´on de la ecuaci´on de Schr¨odinger. En el caso que nos ocupa, este fen´omeno se conoce como hibridaci´on y, aunque puede darse entre cualquier ¯ grupo de orbitales, son especialmente conocidos los que se dan entre los orbitales s y p, y en concreto la hibridaci´on 3 sp que se produce especialmente en el carbono originando cuatro orbitales que se distribuyen en el espacio en forma de tetraedro. Los orbitales h´ıbridos resultantes de dicha combinaci´on lineal son: [sp3 ]1 [sp3 ]2 [sp3 ]3 [sp3 ]4

= = = =

1 2 (s 1 2 (s 1 2 (s 1 2 (s

+ px + py + pz ) − px + py − pz ) − px − py + pz ) + px − py − pz )

donde los orbitales s = Ψ2,0,0 , px = Ψ2,1,1 , pz = Ψ2,1,0 y py = Ψ2,1,−1 . Igualmente se pueden obtener hibridaciones sp3 con orbitales de la capa 3, 4, etc. Ejercicio 11 Dibuja los orbitales h´ıbridos [sp3 ]1 y [sp3 ]4 .

4.3.

Orbitales Moleculares

Cuando dos ´atomos se aproximan el uno al otro hasta situarse convenientemente cerca, sus orbitales sufren un solapamiento de manera que se produce algo parecido al fen´omeno de hibridaci´on comentado anteriormente. De esta forma, de dos orbitales que se solapen A y B, se forman otros dos, uno de menor energ´ıa llamado enlazante ΨE y otro de mayor energ´ıa llamado antienlazante ΨA . ΨE = N (A + B)

ΨA = N (A − B)

donde N es una constante de normalizaci´on. Los orbitales A y B son de ´atomos diferentes, y por tanto se evaluar´an en puntos diferentes del espacio. Seg´ un la teor´ıa de los Orbitales Moleculares, los electrones se sit´ uan dentro de estos orbitales que forman las mol´eculas que conocemos. Si el solapamiento se realiza siguiendo la direcci´on del eje de simetr´ıa de los orbitales, al orbital molecular resultante se le llama orbital σ, mientras que si el solapamiento se realiza de forma perpendicular a los ejes de simetr´ıa se le llama orbital π. El solapamiento de muchos orbitales originar´a en los metales dos macro-orbitales que se extienden a todos los ´atomos de la red met´alica. Se llaman bandas de conducci´on y de valencia. Ejemplo 8 Para dibujar los orbitales moleculares enlazante y antienlazante, resultado del solapamiento de dos orun el eje X habr´ a que modificar el script llamado missuperficiesdenivel.m bitales 2pz situados a una distancia de 6a0 seg´ de manera que eval´ ue el orbital enlazante Psi=FdOn210(x,y,z)+FdOn210(x+6,y,z); o el antienlazante Psi=FdOn210(x,y,z)-FdOn210(x+6,y,z);. Nota: Puedes rotar las figuras resultantes para apreciar mejor los detalles. Ejercicio 12 Dibuja cada uno de los siguientes orbitales moleculares: a) σE = 1s + 1s situados a una distancia de 2a0 en el eje X. b) σA = 1s − 1s situados a una distancia de 2a0 en el eje X. c) σE = 2px + 2px situados a una distancia de 8a0 en el eje Y. d) σA = 2px − 2px situados a una distancia de 8a0 en el eje Y.

Figura 3: Orbitales πE y πA .

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) C´ alculo Simb´ olico. Pr´ actica 6

1.

Introducci´ on

En esta pr´actica vamos a a utilizar los comandos del programa Matlab de C´alculo Simb´olico. En particular vamos a ver c´omo se pueden calcular l´ımites, derivadas e integrales.

2.

L´ımites

Ejemplo 1 C´ alculo de l´ımites. Queremos calcular el l´ımite siguiente:

x2 − 3 x→3 3x5 + 5x l´ım

Procedemos de la manera siguiente: Definimos la variable simb´ olica x: >>syms x Escribimos >>l=limit((x^2-3)/(3*x^5+5*x),x,3) Ejecutamos el comando, nos dar´ a el resultado que en este caso es 1 124 Ejercicio 1 Calcular los siguientes l´ımites: ex − 1 a) l´ım x→0 log(1 + x) x + sen(πx) b) l´ım x→0 x − sen(πx) 1 c) l´ım (ex + x) x x→0 xeax − x , (a ∈ IR) d) l´ım x→0 1 − cos(ax) sen(3x) e) l´ımπ x→ 3 1 − 2 cos x 1 4 ln f ) l´ım (x + 1) x x→∞

3.

Derivadas

Ejemplo 2 Queremos calcular la derivada tercera de la funci´ on: y = 4x sen x Escribimos la funci´ on simb´ olica, despu´es de haber definido la variable x.

>>y=4*x*sin(x) Escribimos el comando >>diff(y,3) Y nos debe dar como resultado: −4x cos x − 12 sen x Si hacemos: >>d3=diff(y,3) >>subs(d3,x,pi/2) π Obtenemos como respuesta −12, es decir, el valor de la derivada tercera de la funci´ on y en x = . 2 Ejercicio 2 µ Calcular ¶ las derivadas primeras y segundas de las siguientes funciones en el punto x = 2: x x2 − (7 − 5x6 )2 a)f (x) = 2 3 b)f (x) = (7x − 5)2

c)f (x) =

4x3 + 3x (x − 1)2

dx2 + ex + f x2 x 1 + + a b c

d)f (x) = (x − a)(x − b)(x − c) 1

e)f (x) = (x + a) 3 (x − a) f )f (x) =

4.

−1 3

p

1 + (x2 − 7)1/2

Integrales

Ejemplo 3 Queremos calcular la integral indefinida Z

x dx x2 + 1

Despu´es de definir la variable simb´ olica x, definimos la funci´ on: >>y=x/(x^2+1) Calculamos la integral con el comando >>int(y) Nos dar´ a la soluci´ on: 1 ln(x2 + 1) 2 Ejercicio 3 Calcular las siguientes integrales indefinidas: Z log x √ dx a) x Z √ b) arctg xdx Z c) Z d) Z e) Z f)

a dx x2 + b2 x2

x+1 dx + 4x + 5

x2 dx x2 − x + 1 x2

x2 dx +x+1

Ejemplo 4 Queremos calcular la integral definida: Z 0

1

x dx x2 + 1

Una vez definida la funci´ on simb´ olica, como en el ejemplo anterior, escribimos el comando >>int(y,0,1) Que nos ofrece la soluci´ on simb´ olica: 1/2*log(2) Si queremos el resultado en forma num´erica lo haremos con el comando quad y una funci´ on definida o bien con un fichero function.m o bien del siguiente modo >>F = inline(’x./(x.^2+1)’); >> Q = quad(F,0,1) Q = 0.3466 Ejercicio 4 Calcular las integrales que se calcularon en el ejercicio anterior, pero haci´endolas todas ellas definidas en el intervalo [10, 12], ofreciendo el resultado en forma algebraica y en forma num´erica. Ejemplo 5 Integrales M´ ultiples. Si queremos el resultado de ¶ Z 2 µZ 1 2 (x + y) dx dy −3

0

en forma num´erica (para l´ımites constantes) lo haremos con el comando dblquad y una funci´ on definida o bien con un fichero function.m o bien >> F=inline(’x.^2+y’) o bien >> F=inline(’x.^2+y’,’x’,’y’) >> Q = dblquad(F,0,1,-3,2) Q = -0.8333 Se pueden calcular integrales con l´ımites variables haciendo que el integrando sea cero fuera de la region. Por ejemplo, ! 2 Z ÃZ 2

−3

10−y

(x2 + y) dx

dy

0

o bien >>F = inline(’(x.^2+y).*(x=0)’); >> Q = dblquad(F, 0, 10, -3, 2) Q = 931.0119 En particular, el volumen de una semiesfera se puede calcular >>F = inline(’sqrt(max(1-x.^2-y.^2,0))’,’x’,’y’); Q = dblquad(F, -1, 1, -1, 1) Si queremos calcular integrales triples en forma num´erica (para l´ımites constantes) lo haremos con el comando triplequad y una funci´ on definida o bien con un fichero function.m o bien como en el ejemplo anterior. Ejercicio ultiples: Z 2 Z51 Calcular las siguientes integrales m´ a) (x2 + y)dxdy −3

0

Z b) Z c)

Z

1

1

(1−x2 ) 2

dydx 0

0

Z

1 −1

Z d)

Z

ex+y dydx

−2|x|

Z

1

1

(1−y 2 ) 2

dxdy 0

Z e)

|x|

0

Z

1

Z

1

−1

−1

1

z

Z

1

Z

y

f) 0

Z g)

0

1

Z

0

Z h)

−∞

(xy 2 z 3 )dxdydz

0

y 0

∞

(x2 + y 2 + z 2 )dxdydz

−1

Z

Z

0 ∞

x √ 3

x2 2

e−x

−∞

x dzdxdy + z2

−y 2

dxdy

Z

Z

∞

i)Sabiendo que Comprueba que

0

y

3

xe−y dxdy =

0

1 , 6

calcula una aproximaci´ on nu´erica. ¿Qu´e ocurre?.

>> Q1 = dblquad(F, 0, 1, 0, 10) Q1 = 0.1481 >> Q2 = dblquad(F, 1, 10, 0, 10) Q2 = 0.0186 >> Q1+Q2-1/6 ans = -6.1393e-006

5.

Polinomios de Taylor

on f (x) = sin(x) en x = Ejemplo 6 Calcular el Polinomio de Taylor de grado 6 de la funci´

π 2.

>>taylor(sin(x),pi/2,6) ans= 1-1/2*(x-1/2*pi)^2+1/24*(x-1/2*pi)^4 Pero tambi´en podemos hacer >>taylortool y seguir las indicaciones. Ejercicio 6 Sea la funci´ on f (x) = ex : a) Calcular los polinomios de McLaurin de grados n=2,4,10 y representarlos, utilizando ezplot, junto a la gr´ afica de la funci´ on. b) Calcular los polinomios de Taylor en torno a x = 1 para los grados n=2,4,10 y representarlos junto a la gr´ afica de la funci´ on.

Ejercicio 7 Sea la funci´ on f (x) = sen(x): a) Calcular los polinomios de McLaurin de grados n = 3, 5, 11 y representarlos junto a la funci´ on. b) Calcular los polinomios de Taylor en x = π2 para los grados n = 3, 5, 11 y representarlos junto a la funci´ on. c) Utilizar el polinomio de McLaurin con n = 3 y el polinomio de Taylor, tambi´en con n = 3, en torno a x = π2 para calcular aproximadamente en valor de sen(0,1). Compararlo con los valores reales. Lo mismo con el valor de sen(1,5). ¿En qu´e casos nos se obtiene una mejor aproximaci´ on? ¿Por qu´e? Ejercicio 8 Investigar qu´e hace el comando >>funtool Hacer help y resolver el problema que se plantea en ´el: (“The Demo button poses the following challenge...”)

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) An´ alisis de Datos Pr´ actica 7

1.

M´ınimos Cuadrados.

El M´etodo de ajuste por M´ınimos Cuadrados sirve para encontrar una funci´on y = f (x, α1 , α2 , . . . , αm ), en la que habr´a que calcular los par´ametros α1 , α2 , . . . , αm . Esta funci´on debe ser la que se ajuste lo mejor posible a una tabla de valores que relaciona las dos variables x e y obtenida experimentalmente: x x1 y y1

x2 y2

... ...

xn yn

Para calcular los par´ametros se impone la condici´on de que sea m´ınima la funci´on S(α1 , α2 , . . . , αm ) =

n X

[f (xi , α1 , α2 , . . . , αm ) − yi ]

2

i=1

Como S(α1 , α2 , . . . , αm ) es una funci´on de m variables, una condici´on necesaria para que tenga un valor extremo en un punto es que sus derivadas parciales en ese punto sean todas nulas. De aqu´ı obtenemos un sistema de m ecuaciones, con m inc´ognitas: ∂S = 0, ∂α1

∂S = 0, ∂α2

...,

∂S =0 ∂αm

cuyas soluciones, los par´ametros α1 , α2 , . . . , αm nos indican c´omo es la funci´on que mejor se ajusta a los datos, es decir, f (x, α1 , α2 , . . . , αm ). La funci´on f (x, α1 , α2 , . . . , αm ) puede ser de cualquier tipo, te´oricamente, sin embargo, en la pr´actica, con programas como MATLAB (o Matlab) s´olo se pueden calcular (directamente) cuando la funci´on es un polinomio. En otros casos, que veremos en los ejercicios, habr´a que modificar ligeramente el problema para que se pueda tratar con el ordenador. Para realizar estos c´alculos MATLAB dispone del comando polyfit. Veamos en un ejemplo c´omo se puede ajustar una recta con M´ınimos Cuadrados: Ejemplo 1 Dada la tabla de valores, x 1 2 3 4 y 2,1 4,3 6 7,8 Encontrar la funci´ on de la forma y = ax + b que mejor se ajuste a los datos. Como se trata de una funci´ on polin´ omica se puede hacer directamente. Introducimos primero la tabla de valores en dos variables: >>x=[1 2 3 4] >>y=[2.1 4.3 6 7.8] El comando a utilizar es polyfit(x,y,k), donde k es el grado del polinomio que queremos obtener. Por lo tanto para obtener una recta k = 1: >>c=polyfit(x,y,1) Que nos da como resultado los coeficientes de la recta: c = 1.8800

0.3500

Es decir, que la recta que hemos encontrado es, y = 1,88x + 0,35 Para representar la informaci´ on obtenida gr´ aficamente: Primero dibujamos la tabla de valores, por ejemplo: >>plot(x,y,’*’) De esta forma conseguimos que dibuje s´ olo los puntos, con asteriscos o con cualquier otro formato.

Para dibujar la recta, lo hacemos como para dibujar cualquier funci´ on. Generamos una tabla, (que llamaremos 8 con un nombre diferente de x para no borrar la tabla de los datos del problema): >>xp=linspace(1,4,20); Para para calcular los 7 valores de xp en la recta y = 1,88x + 0,35 podemos utilizar el comando polyval que eval´ ua el polinomio utilizando los coeficientes, que ten´ıamos en la variable c: >>yp=polyval(c,xp);6 >>hold on % para mantener el dibujo anterior Y, por u ´ltimo, 5 plot(xp,yp) y=1.88x+0.35

4

3

2

1

1.5

2

2.5

3

3.5

4

Figura 4: Tabla y Recta y = 1,88x + 0,35 El error que se comete al aproximar la funci´ on emp´ırica (tabla inicial) por la funci´ on te´ orica (recta) se puede cuantificar de varias formas. Una manera es precisamente calcular la suma de las desviaciones en cada punto de la tabla al cuadrado, es decir, el valor de la funci´ on: S(α1 , α2 , . . . , αm ): Primero calculamos los valores de f (xi , α1 , α2 , . . . , αm ): >>fxi=polyval(c,x) fxi = 2.2300

4.1100

5.9900

7.8700

Y ahora sustituimos en la f´ ormula: >>error=sum((fxi-y).^2) error = 0.0580 Otra forma consiste en calcular lo que se llama Error Cuadr´atico Medio, que viene dado por la expresi´ on: v u n u1 X 2 ε=t [f (xi , α1 , α2 , . . . , αm ) − yi ] n i=1 Para calcularlo: >>n=length(x) % calcula el n’umero de t’erminos n = 4 >>errorcm=sqrt(sum((fxi-y).^2)/n)

errorcm = 0.1204 Por u ´ltimo, se pueden utilizar los c´ odigos que llevan incorporados los comandos polyfit y polyval para obtener un intervalo de confianza al 95 % en el que se encuentra cada yi . Veamos en qu´e consiste: Primero volvemos a utilizar el comando polyfit, de una forma diferente: >>[c,s]=polyfit(x,y,1) c = 1.8800 0.3500 s = R: [2x2 double] df: 2 normr: 0.2408 s es la estructura con la que se va a calcular el error. Ahora volvemos a utilizar polyval, pero tambi´en modificado: >>[fxi,delta]=polyval(c,x,s) fxi = 2.2300 4.1100 5.9900 delta = 0.2220 0.1942 0.1942

7.8700 0.2220

Hemos obtenido, en cada xi , un intervalo de confianza al 95 %, [fxi-2*delta,fxi+2*delta], dentro del cual se encuentra el valor yi , i. e. ³ ´ P |f (xi ) − yi | ≤ 2δi = 0,95 donde P denota la probabilidad. Si representamos la gr´ afica siguiente se entender´ an de manera gr´ afica los c´ alculos que acabamos de realizar. >>plot(x,y,’*’,x,fxi,’g-’,x,fxi+2*delta,’r:’,x,fxi-2*delta,’r:’) Ejercicio 1 Mediante el M´etodo de M´ınimos Cuadrados, elegir una funci´ on cuadr´ atica del tipo f (x) = ax2 + bx + c para los valores de x e y dados por la siguiente tabla x y

7 7,4

8 8,4

9 9,1

10 9,4

11 9,5

12 9,5

13 9,4

Calcular la funci´ on f (x). El error, el error cuadr´ atico medio y representar la funci´ on obtenida junto a los datos de la tabla con una banda alrededor de la gr´ afica de f (x), dentro de la cual se encuentren los puntos de la tabla a un nivel del 95 %. Ejercicio 2 De ciertas medidas realizadas se han obtenido los siguientes valores: x y

0 1

0,2 0,833

0,5 0,667

1 0,54

1,5 0,4

2 0,333

2,5 0,286

3 3,5 4 0,25 0,222 0,2

4,5 5 0,182 0,167

a) Representar la tabla de valores gr´ aficamente.

1 . Esta funci´ on no la ax + b podemos obtener directamente con el comando polyfit. Pero la podemos transformar en, b) Se puede observar que se trata aproximadamente de una funci´ on del tipo y = 1 = ax + b y Ahora, haciendo X = x; Y =

1 . Nos queda la expresi´ on y Y = aX + b

Calcular esta funci´ on. Habr´ a que hacer una nueva tabla de valores con las nuevas variables. c) Deshacer el cambio de variable y obtener la funci´ on original y representarla con los valores del apartado a). d) Calcular el error y el error cuadr´ atico medio.

Ejercicio 3 Hallar por el M´etodo de M´ınimos Cuadrados una funci´ on del tipo S = Atq para los datos de la siguiente tabla. Representarla junto a los datos de la tabla: t S

1 7,1

2 27,8

3 62,1

4 110

5 161

(Indicaci´ on: Tomar logaritmos neperianos en la funci´ on para transformarla en una funci´ on de tipo lineal). Ejercicio 4 La tabla siguiente contiene la presi´ on p en kp/cm2 de un vapor saturado, en funci´ on del volumen 3 espec´ıfico, v, medido en m /kp: v p

3,334 0,482

1,63 1,034

0,8657 2,027

0,4323 4,247

0,2646 0,1699 0,1146 7,164 11,48 17,6

Elegir una funci´ on p = f (v) que sea adecuada a la tabla de valores y calcularla con el M´etodo de M´ınimos Cuadrados. Representarla despu´es con los valores. Seg´ un esta funci´ on ¿Cu´ al ser´ıa la presi´ on aproximada que corresponder´ıa a un volumen de 3m3 /kp?

2.

Analisis Multivariante. Regresi´ on

Se utiliza cuando y es una funci´on de m´as de una variable independiente, las ecuaciones matriciales que expresan las relaciones entre las variables se deben ajustar a esos datos. Esto se conoce como regresi´on m´ ultiple. Supongamos que medimos una cantidad y para varios valores de x1 y x2 . >> x1 = [.2 .5 .6 .8 1.0 1.1]’; >> x2 = [.1 .3 .4 .9 1.1 1.4]’; >> y = [.17 .26 .28 .23 .27 .24]’; Un modelo lineal de regresi´on calcula los coeficientes , y = a0 + a1 x1 + a2 x2 , por el m´etodo de los m´ınimos cuadrados. Construye y resuelve el sistema de ecuaciones formando una matrix X y calcula los par´ametros con el operador \ del siguiente modo >> X = [ones(size(x1)) >> a = X\y >> a = 0.1018 0.4844 -0.2847

x1 x2];

i.e. el modelo de ajuste por m´ınimos cuadrados es y = 0,1018 + 0,4844x1 − 0,2847x2 . Para validar el modelo, calcula el m´aximo del valor absoluto de las desviaciones de los datos a partir del modelo >> Y = X*a; >> MaxErr = max(abs(Y - y)) >> MaxErr = 0.0038 el error relativo est´a entre los l´ımites aceptables. Observaci´ on: Se ha desarrollado este m´etodo con el comando \ de matlab. Podr´ıamos hacerlo con el comando pinv. ales son las semejanzas y diferencias entre Ejercicio 5 Repite el ejemplo anterior utilizando el comando pinv. ¿Cu´ ambos m´etodos?. Recuerda lo comentado en la pr´ actica 2. Ejercicio 6 Ajustar los datos siguientes x 1 2 3 4 y 0 1 2 3 z −2,1 1,3 4,2 5,6 a una funci´ on lineal del tipo z = ax + by + c. Ejercicio 7 Repite el ejemplo anterior utilizando el comando pinv. ¿Por qu´e se obtienen dos soluciones distintas?.

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) Ecuaciones Diferenciales Ordinarias. Pr´ actica 8 En esta pr´actica vamos a resolver Ecuaciones Diferenciales Ordinarias (E.D.O.), calculando aproximaciones num´ericas de las soluciones (usando MATLAB). Algunos ejemplos de ecuaciones diferenciales ordinarias en las ciencias f´ısicas son 1.

dx = kx dt

2.

dx = k(a0 − x)(b0 − x) dt

3.

d2 h = −g dt2

4.

d2 x = −ωx2 dt2

5. −

proceso cin´etico de primer orden proceso cin´etico de segundo orden

cuerpo que cae bajo la acci´on de la gravedad oscilador arm´onico cl´asico

h2 d2 ψ 1 2 + kx ψ = Eψ 8π 2 m dx2 2

ecuaci´on de ondas para el oscilador arm´onico

Ejemplo 1 Un proceso cin´ etico de primer orden. Considera la ecuaci´on diferencial x0 = −10x + e−t ,

x(0) = 0,

en el intervalo, [0, 4] on, con el nombre funcorden1.m, para la funci´on de esta ecuaci´on diferencial, 1. Escribir un fichero tipo funci´ es decir la expresi´on f (t, x), en concreto: f u n c t i o n f=f u n c o r d e n 1 ( t , x ) f =−10∗x+exp(− t ) ; 2. Usar, en la lineas de comandos, la rutina ode45 de Matlab para resolver la ecuaci´on num´ericamente. Dibuja la gr´afica de la soluci´on aproximada. En concreto, escribe >> [ t , u]= ode45 ( @funcorden1 , [ 0 , 4 ] , 0 ) ; >> p l o t ( t , u ) Ejercicio 1 Considera una reacci´ on reversible de primer orden en ambas direcciones A * ) B: x0 (t) = k1 (a − x) − k−1 x,

x(0) = 0,

t ∈ [0, 10].

1. Repetir el punto 1 del ejemplo anterior para a = k1 = 1, k−1 = 0,1. El fichero funci´on que contiene la funci´on de la ecuaci´on diferencial se tiene que llamar funcreversible11.m. 2. Escribe en una linea todos los datos de entrada del comando ode45, es decir: >> f=@ f u n c r e v e r s i b l e 1 1 ; x0 =0; i n t e r v a l o = [ 0 , 1 0 ] ; N=1000; 3. Repetir el punto 2 del ejemplo anterior. Opcional: Usa el comando title para poner un t´ıtulo a la ventana y el comando legend para indicar la(s) curva(s). k

Ejercicio 2 Consideramos la reacci´ on de tercer orden completa A+B +C *D. La concentraci´on de D verifica la ecuaci´on diferencial d0 (t) = k(a0 − d(t))(b0 − d(t))(c0 − d(t)), d(0) = 0. 1. Elige valores de los par´ametros tales que a0 > b0 > c0 > 0. 2. Resuelve num´ericamente la ecuaci´on diferencial. 3. ¿Qu´e ocurre con la soluci´on si a0 = b0 o b0 = c0 ?.

4. ¿Y si a0 = b0 = c0 ?. 5. Establece la analog´ıa con las reacciones k

2A + B *D y k

3A*D. 6. ¿Que cambia en lo anterior si suponemos unas relaciones estequiom´etricas en la reaccion de la forma k

nA A + nB B + nC C *nD D ?. Ejercicio 3 Consideremos una reacci´ on de primer orden autocatal´ıtica k

A*B. En dicha situaci´on la constante de proporcionalidad, al aplicar la Ley de Acci´on de Masas, verifica que k = k([B]) = k0 + α[B] siendo α una constante, usualmente peque˜ na. Si α > 0 se dice que B es un activador de la reacci´on mientras que es un inhibidor si α < 0. 1. Obtener la ecuaci´on diferencial para las concentraciones de A y B y resolverlas supuesto que la concentraci´on inicial de B es nula. Determinar cual es la velocidad m´axima que alcanza la reacci´on. Determinar la concentraci´on final de las sustancias. dx = (k0 + αx)(a0 − x) dt con dato inicial x(0) = 0. 2. Integrar la ecuaci´on diferencial de la velocidad de reacci´on de la siguiente reacci´ on incompleta dx = 3(2 − x)(1 − x) − (1 + x)(2 + x) dt con dato inicial x(0) = 0. ¿Para qu´e valor de la variable x cesa ´esta su variaci´on?. Ejemplo 2 Considera el sistema de ecuaciones diferenciales ½ 0 x = x + 7y − 10x2 y y 0 = −3x − 8y + xy

x(0) = 1 y(0) = 0,

en el intervalo, [0, 3] 1. Escribir un fichero tipo funci´ on, con el nombre funcsistema1.m, para la funci´on de este sistema de ecuaciones diferenciales function f=funcsistema1(t,x) f=zeros(2,1); f(1)=x+ 7*y - 10*x.^2.* y; f(2)=-3*x -8*y + x.*y; 2. Usar, en la lineas de comandos, la rutina ode45 de Matlab para resolver la ecuaci´on num´ericamente. Dibuja la gr´afica de la soluci´on aproximada. En concreto, escribe >> [ t , u]= ode45 ( @funcsistema1 , [ 0 , 3 ] , [ 1 ; 0 ] ) ; Ejercicio 4 La enzima [E] reacciona reversiblemente con el sustrato [S] y forma el complejo [C]. A continuaci´on el complejo se descompone, en una reacci´on irreversible, en un producto [P ] y la enzima original k1

k

2 E + S* ) C *E + P.

k−1

La ley de acci´on de masas para las concentraciones e(t) = [E], s(t) = [S], c(t) = [C], p(t) = [P ], ser´a e0 (t) s0 (t) c0 (t) p0 (t)

= = = =

−k1 e(t)s(t) + k−1 c(t) + k2 c(t) −k1 e(t)s(t) + k−1 c(t) k1 e(t)s(t) − k−1 c(t) − k2 c(t) k2 c(t)

con las condiciones iniciales e(0) = e0 ,

s(0) = s0 ,

c(0) = p(0) = 0.

(2) (3) (4) (5)

1. Demostrar que la suma de concentraciones e(t) + c(t) = e0

(6)

es constante. Elimina la variable e(t) en las ecuaciones (3)-(4). 2. Dando valores iniciales, resolver num´ericamente las ecuaciones (3)-(4). 3. Estado pseudo-estacionario. En los experimentos en laboratorio, s0 À e0 . Michaelis y Menten en 1913 observaron que en estas circunstancias se espera que, despu´es de un corto intervalo de tiempo, haya un balance entre la formaci´on del complejo [C] por la uni´on de enzima y sustrato y la degradaci´on del complejo. Debido a que hay una gran cantidad de moleculas de sustrato s0 À e0 , este balance se alcanzar´a antes de sea perceptible la transformaci´on de sustrato en producto. La formaci´on de producto se puede llevar a cabo con la hip´otesis c0 (t) = 0 es decir k−1 + k2 (7) e(t)s(t) = Km c(t), Km := , k1 Km es la constante de Michaelis. Esta ecuaci´on se conoce como la hip´otesis del estado pseudo-estacionario, c(t) est´a ajustada a los valores instant´aneos de e(t) y de s(t) y estos valores cambian lentamente con el tiempo. Bajo la hip´otesis pseudo-estacionaria, sustituyendo (6) en (7), comprueba que se puede escribir c(t) = en consecuencia s0 (t) = − Resuelve num´ericamente este PVI. 4. Compara los resultados con el caso anterior.

e0 s(t) Km + s(t)

k2 e0 s(t) , Km + s(t)

s(0) = s0 .

(8)

Departamento de Matem´atica Aplicada ´ CALCULO COMPUTACIONAL. Licenciatura en Qu´ımica (Curso 2008-09) Archivos de ´ ordenes. Programaci´ on. Pr´ actica 9

1.

Introducci´ on

Hasta ahora, todos los comandos que hemos estudiado en MatLab, los hemos ido escribiendo en el entorno de trabajo del programa, es decir, a continuaci´on de >>. Sin embargo, en muchas ocasiones resulta necesario ejecutar varios o muchos comandos encadenados y, quiz´as, repetirlos varias veces. Para esto resulta u ´til guardarlos en un fichero, de manera que se puedan ejecutar con s´olo una orden, esto es, queremos construir un programa, que no es otra cosa que un conjunto de comandos encadenados. La manera de guardar comandos y, por tanto, de construir programas que ejecuten determinadas rutinas consiste en crear lo que se llaman ficheros-M. Para crear un fichero-M necesitamos abrir un editor de texto (como el Bloc de Notas de Windows) y escribir los comandos. Despu´es, ese fichero de texto debe guardarse con la extensi´on .m. La versi´on 5 (Student Edition) incorpora un editor de Ficheros-M. Basta abrir la opci´on File New M-File de la barra superior. Vamos a ver algunos ejemplos de ficheros-M: Ejemplo 1 .

1 1 1 a) Vamos a elaborar un programa para calcular S(n) = 1 + 2 + 2 + · · · + 2 . 2 3 n Creamos el siguiente fichero:

function S=suma(n) % Fichero para la funcion suma % de los inversos de los cuadrados N=1:n; NN=N.^2; Ninv=1./NN; S=sum(Ninv); Y lo guardamos con el nombre suma.m. Ahora lo podemos utilizar para calcular valores de la suma >>suma(100) ans = 1.6350 Z Ejercicio 1 Elaborar un programa para calcular la integral definida de una funci´ on, punto medio: si se divide [a, b] en n subintervalos de longitud ∆x =

b

f (x)dx, por el m´etodo del a

b−a , N

la integral aproximada es I=

N X

f (xi )∆x

i=1

donde xi es el punto medio de cada subintervalo. Los datos de entrada ser´ an la funci´ on f, introducida como un fichero f.m, los extremos a, b del intervalo y el n´ umero N de sumandos. El nombre del fichero es miptomedio.m Ejemplo 2 a) Elabora el fichero func1.m para la funci´ on sen(x)/x, es decir: function y=func1(x) % Fichero para la funcion sen(x)/x y=sin(x)./x; b) Teniendo en cuenta los datos de entrada del programa miptomedio, es decir: >> f=@func1 ; a =0; b=1; N=1000;

Z calcula: 0

1

sen(x) dx, es decir: x

>> I=miptomedio ( @func1 , 0 , 1 , 1 0 0 0 ) Ejercicio 2 Calcula: Z π/4 tan(x) a) dx x 0

2.

Z

π

b)

sen(x2 ) dx

0

1 c) √ π

Z

N −N

Z 2

e−x dx

d) 0

1

5x3 + 2 dx 1 + x + x2

Control de flujo. Bucles

Cuando los comandos los ejecutamos en el entorno de trabajo, las instrucciones se van realizando en el orden en el que se van introduciendo. Pero, dentro de un programa, el orden de su ejecuci´on, el flujo, se puede alterar utilizando algunas instrucciones. Estas instrucciones pueden ser: condicionales, en las que las secuencias de ´ordenes deben evaluarse bas´andose en alguna condici´on y, bucles o instrucciones iterativas, mediante las cuales una o un grupo de ´ordenes se ejecutan varias veces.

Condicionales Vamos a ir viendo las diferentes posibilidades mediante ejemplos sencillos: Estructura if simple. Su sintaxis es if condicional comando end El comando se ejecuta si todos los elementos en condicional son verdaderos. Por ejemplo, >>a=3; if a>2 b=a+5 end Mientras que si escribimos >>a=3; if a==2 b=a+5 end no obtenemos respuesta. Estructura if compuesta, if condicional comando primero else comando segundo end Sirve para los casos en los que existen dos alternativas. Si se cumple condicional, se ejecuta comando primero, pero si es de otra forma, se ejecutar´a comando segundo. Por ejemplo, >>a=n; if a>5 b=a+5 else b=a-1 end

En la secuencia anterior de ´ordenes, si n es mayor que 5, el valor de b ser´a a+5, de otra forma, ser´a a-1. Estructura if-else-if. Se puede utilizar cuando hay m´as de dos condiciones que pueden cumplirse. Se escribe if condicional primero comando primero elseif condicional segundo comando segundo elseif condiccional tercero comando tercero ...... else comando final end En este caso el comando final se ejecuta, si no se han realizado ninguno de los anteriores. Por ejemplo, if n>a a = 3

4

5

6

7

8

9

10

11

12

N´otese que ahora hay que poner ; detr´as del comando que se ha de ejecutar si no queremos ver todos los c´alculos en la pantalla. Por supuesto, lo anterior se podr´ıa haber obtenido simplemente escribiendo >>a=3:12 Los bucles, dentro de un programa, consumen tiempo y memoria del ordenador, en la medida de lo posible, hay que tender a evitarlos. Aunque en ocasiones puede que no haya, o no se nos ocurra, otra posibilidad. El bucle de tipo for ejecuta un conjunto de instrucciones un n´ umero de veces que nosotros conocemos a priori. Pero esto puede no ser as´ı, en este caso habr´ıa que utilizar un bucle de tipo while. Su sintaxis es while condicional comando end Que significa que, mientras se cumpla condicional, se ejecutar´a comando, cuando deje de cumplirse condicional, comando dejar´a de ejecutarse.

3.

M´ etodo de Newton

El m´etodo de Newton sirve para aproximar ceros de funciones f (x) = 0. Las funciones f pueden ser de varias variables f : IRk → IRk . Teniendo en cuenta el desarrollo en serie de Taylor de orden 1 podemos escribir f (xn+1 ) = f (xn ) + Df (xn )(xn+1 − xn ) + t.o.s. donde t.o.s. representa a los t´erminos de orden superior, y suponiendo xn+1 ) ≈ f (xn ). El m´etodo de Newton-Raphson es un m´etodo iterativo Dado un valor inicial x = x0 ½ resolver Df (xn )w = f (xn ), calcular xn+1 = xn − w,

que f (xn+1 ) ≈ 0 tendr´ıamos Df (xn )(xn − que se puede describir como: n≥0 n ≥ 0.

De esta forma si xn converge cuando n → ∞ a x,µel l´ımite¶es una soluci´on del sistema no lineal f (x) = 0. Obs´ervese ∂fi (x) que la matriz Df (xn ) viene dada por Df (x) = ∂xj i,j Ejemplo 3 Vamos a escribir en un fichero instrucciones de MATLAB de forma que, utilizando el Algoritmo de Newton, calculemos aproximaciones a un cero de la funci´ on f (x) = x − x3 . function x=minewton(f,df,x0,TOL) % metodo de Newton para ecuaciones, dadas f

y su derivada df

% variables de entrada: f,df ficheros .m % x0 dato inicial vector columna % TOL tolerancia error % variables de salida: x aproximacion al punto fijo en ’i=contador’ iteraciones ERR=TOL+1; x=x0; contador=0; while ERR>TOL contador=contador+1; DF=feval(df,x); F=feval(f,x); w= DF\F; ERR=norm(w,inf); x=x-w; if contador>1000 break disp(’excede numero maximo de iteraciones’) end end Llamad a este fichero minewton.m. Ahora podemos ejecutarlo. Para ello necesitaremos en primer lugar definir un fichero .m con la correspondiente funci´ on f y otro con su derivada, df, como una matriz si f fuera de varias variables. Para ello escribimos en dos ficheros: function F=fun1(x) % F=x-x.^3; function DF=dfun1(x) % DF=1-3*x.^2; y los guardamos con el nombre fun1.m y dfun1.m respectivamente. Esto nos permitir´ a evaluar la funci´ on y su on escribimos las variables de entrada en la linea de comandos derivada para cada x. A continuaci´ >>TOL=1e-10; >>x0=-2; y ahora llamamos al programa minewton.m desde la linea de comandos >> x=minewton ( @fun1 , @dfun1 , x0 ,TOL) ;

Ejercicio 3 Calcule los ceros de las siguientes funciones utilizando el m´etodo de Newton i) f (x) = x − x2 , ii) f (x) = xe−x + 1, iii) f (x1 , x2 ) = (x1 − x22 , sin(x1 ) − 0,7x2 ).