¡Descarga Laboratorio de sistemas de control tutoriales y más Guías, Proyectos, Investigaciones en PDF de Sistemas de Control solo en Docsity!
0.1 Ejercicios de Scilab
- Escriba una funci´on [A, info] = Heron(a, b, c) que calcula, usando la f´ormula de Her´on de Ale- jandr´ıa, el ´area de un tri´angulo cuyos lados miden a, b, c. El par´ametro de salida info valdr´a 1 si los datos est´an bien. En esta y en todas las funciones que usted escriba: - Si los datos pueden ser inadecuados o si el proceso puede no acabar satisfactoriamente debe haber un par´ametro de salida que lo indique, por ejemplo info. - En los comentarios iniciales debe quedar claro lo que hace la funci´on, el m´etodo utilizado, el significado de los par´ametros de entrada (datos) y de salida (resultados). - Cuando se presente alg´un error no previsto, adem´as de influir en el par´ametro info, la funci´on debe escribir un aviso indicando el tipo de error, printf(... ).
- Escriba una funci´on d = dist(x, y) que calcula la distancia entre dos puntos.
- Escriba una funci´on a = area(x, y, z) que calcula el ´area de un tri´angulo cuyos v´ertices son x, y, y z.
- ¿Cual es el ´area del tri´angulo cuyos v´ertices est´an dados por los 6 ´ultimos d´ıgitos de su documento de identidad?
- Escriba una funci´on [xmin, fmin] = minimo(f, a, b, h) que halla el minimizador de una funci´on en un intervalo, simplemente comparando todos los valores f (a), f (a + h), f (a + 2h), ... Muestre el resultado de su utilizaci´on en un ejemplo espec´ıfico.
- Escriba una funci´on [r, info] = raiz(f, a, b, h, eps) que halla una ra´ız de la funci´on en el intervalo, simplemente averiguando si entre los los valores f (a), f (a + h), f (a + 2h), ... hay alguno tal que |f (x)| ≤ ε. Muestre el resultado de su utilizaci´on en un ejemplo espec´ıfico.
- Escriba una funci´on [a, info] = ordenaSegun(a, j) que ordena de manera ascendente las filas de una matriz seg´un los elementos de la columna j. Muestre el resultado de su utilizaci´on en un ejemplo espec´ıfico.
0.2 Procesos inestables
- En un archivo .sci escriba una funci´on function y = redDec(x, n) que redondea a x con n cifras decimales.
- Verifique si la siguiente funci´on redondea un n´umero con n cifras significativas.
function y = fun_xyz(x, n) y = 0. if x == 0, return, end
if x >= 0 signo = 1 else signo = - x = -x end
m = ceil( log10(x) ) k = 10^m t = x/k c = 10^n y = round(tc)/c y = yk*signo endfunction
Pru´ebela varias veces, con n´umero positivos, negativos, con cero, con n´umeros cercanos a cero, con n´umero grandes, con n´umeros negativos de valor absoluto grande, ...
- Escriba una funci´on function z = sumaRed(x, y, n) que redondea a x y a y, efect´ua la suma x + y y redondea el resultado con n cifras significativas.
- Escriba funciones function z = prodRed(x, y, n) function z = restaRed(x, y, n) function z = divRed(x, y, n) function z = raiz2Red(x n) que efect´uan el producto, la resta, la divisi´on y la raiz cuadrada con n cifras significativas.
- Considere la sucesi´on definida por
x 0 = 1 x 1 = 1/ 3
xn =
xn− 1 −
xn− 2 , n ≥ 2. (1)
Demuestre por inducci´on que
xn =
3 n^ , n = 0, 1 , 2 , ... (2)
- Haga un progama, archivo .sce, que produzca una tabla con cinco columnas:
0.3 M´etodo de Gauss sin y con pivoteo
El objetivo es escribir una funci´on que implemente el m´etodo de Gauss sin pivoteo, una funci´on que im- plemente el m´etodo de Gauss con pivoteo parcial y comparar los resultados. En ambos casos, en todas las operaciones se redondea con un n´umero fijo de cifras. Puede ser con un n´umero fijo de cifras significativas o con un n´umero fijo de cifras decimales.
- Escriba funciones, redond(x, n), sumaRed(x, y, n), multRed(x, y, n), divRed(x, y, n), restaRed(x, y, n), redMatr(a, n), para redondear un n´umero con n cifras, sumar (redondea los n´umeros, suma y redondea el resultado), multiplicar, dividir, restar y redondear una matriz.
- Escriba funciones [x, info] = solTriSupR(a, b, nc) [a, b, info] = triangR(a0, b0, nc) [x, info] = GaussR(a, b, nc) que, respectivamente, resuelve un sistema triangular superior, triangulariza un sistema de ecuaciones y resuelve un sistema de ecuaciones por el m´etodo de Gauss. En las tres funciones las operaciones se hacen utilizando redondeo con nc cifras. Siempre, por convenci´on, el par´ametro de salida o resultado info valdr´a 1 si y solamente si el proceso se desarroll´o satisfactoriamente. No olvide incluir en la primera funci´on (en el sentido de su utilizaci´on en el tiempo) el redondeo de los datos. Puede ensayar la funci´on con solTriSupR con ejemplos construidos “a mano” o con datos aleatorios. Por ejemplo, puede usar ´ordenes semejantes a
n = 4; a = (rand(n,n) - 0.5)20; a = triu(a); x0 = (rand(n,1) - 0.5)20; b = a*x0; xTS = solTriSup0(a, b); errTS = norm(xTS-x0)
El valor de errTS debe ser cero o casi cero. Para GaussR puede utilizar ´ordenes semejantes a las anteriores, excluyendo a = triu(a). Para escribir las 3 funciones, pueden ser ´utiles pero no indispensables, los siguientes pasos:
- escribir funciones x = solTriSup0(a, b), [a, b] = triang0(a0, b0), x = Gauss0(a, b) que no controlan la existencia de pivotes nulos ni la posiblidad de que la matriz no sea invertible.
- escribir funciones [x, info] = solTriSup(a, b), [a, b, info] = triang(a0, b0), [x, info] = Gauss(a, b) que controlan la existencia de pivotes nulos, que hacen intercambio de filas cuando es necesario y que consideran la posibilidad de que la matriz sea no invertible. En estas funciones se pueden utilizar las operaciones matriciales de Scilab.
- escribir funciones [x, info] = solTriSup1(a, b), [a, b, info] = triang1(a0, b0), [x, info] = Gauss1(a, b), semejantes a las 3 anteriores pero que no utilizan las operaciones matriciales de Scilab. Por ejemplo, la orden a(i,i+1:n)x(i+1:n) se puede remplazar por t = 0 for j = i+1:n t = t + a(i,j)x(j) end
- Escriba funciones [a, b, info] = triangPP_R(a0, b0, nc) [x, info] = GaussPP_R(a, b, nc) que, respectivamente, triangulariza un sistema de ecuaciones con pivoteo parcial y resuelve un sistema
de ecuaciones por el m´etodo de Gauss con pivoteo parcial. En estas funciones las operaciones se hacen utilizando redondeo con nc cifras. Puede ser ´util escribir funciones previas an´alogas a las indicadas para las funciones triangR, etc.
- Construya un ejemplo 2 × 2 o 3 × 3 o 4 × 4, donde:
- no haya diferencia (o muy peque˜na) al utilizar GaussR y GaussPP_R con 15 cifras significativas.
- haya diferencia clara al utilizar GaussR y GaussPP_R con un n´umero de cifras entre 4 y 12.
- Para resolver una sistema de ecuaciones Ax = b, con matriz definida positiva, es preferible, en lugar del m´etodo de Gauss, utilizar el m´etodo de Cholesky, es decir - Obtener U tal que U TU = A. - Resolver U Ty = b. - Resolver U x = y.
¿Cual es, aproximadamente, el n´umero de operaciones de punto flotante requeridas para el proceso anterior?
- Obtenga una estimaci´on del tiempo requerido para la soluci´on de un sistema con matriz definida positiva mediante las tres ´ordenes
u = chol(a) y = (u’)\b x = u\y
- Compare el tiempo anterior con el tiempo de x = a\b. Concluya.
- Escriba una funci´on x = solTriSup(a, b) que resuelve un sistema de ecuaciones con matriz triangular superior. Suponga provisionalmente que los tama˜nos son compatibles, que la matriz s´ı es triangular superior y que no hay elementos diagonales nulos.
- Escriba una funci´on x = solTriInf(A, b) que resuelve un sistema de ecuaciones con matriz triangular inferior.
- Escriba una funci´on x = sol_Ut_b(U, b) que resuelve el sistema triangular inferior U Tx = b donde U T^ es triangular inferior, es decir, U es triangular superior. Esta funci´on es muy parecida a la anterior, pero la informaci´on sobre la matriz est´a en diferente sitio.
- Escriba la funci´on x = solDP(a, b) que resuelve un sistema de ecuaciones con matriz definida positiva mediante las ´ordenes
u = chol(a) y = sol_Ut_b(U, b) x = solTris(u, y)
- Obtenga una estimaci´on del tiempo de la funci´on anterior. Compare con el tiempo de x = a\b. Concluya.
0.5 M´etodo de Gauss-Seidel
Se desea resolver un sistema de ecuaciones lineales Ax = b. El algortimo para el m´etodo de GS se puede escribir as´ı:
datos: A, b, x^0 , ε, maxit x = x^0 para k = 1, ...,maxit nrmD← 0 para i = 1, ..., n δi = (bi − Ai· x)/aii xi ← xi + δi nrmd←nrmD+|δi| fin-para i si nrmD ≤ ε ent x∗^ ≈ x, salir fin-para k
- Escriba una funci´on [x, k, info] = GS(a, b, x0, eps, maxit) que implementa el m´etodo de GS. Los resultados ser´an: x aproximaci´on de la soluci´on (posiblemente), k el n´umero de iteraciones com- pletas realizadas e info que valdr´a 1 si se obtuvo una buena aproximaci´on de la soluci´on.
- Considere una matriz A de la forma
A =
d 1 0 v 1 0 0 0 · · · 0 d 2 0 v 2 0 0 u 1 0 d 3 0 v 3 0 0 u 2 0 d 4 0 v 4
dn− 2 0 vn− 2 0 dn− 1 0 un− 2 0 dn
En palabras, los elementos diagonales de A son d 1 , d 2 , ..., dn; en la segunda superdiagonal est´an los valores v 1 , v 2 , ..., vn− 2. en la segunda subdiagonal est´an los valores u 1 , u 2 , ..., un− 2. Los otros elementos (entradas) de A son nulos. Escriba una funcion A = matriz3(d, u, v) que, dados los tres arreglos d, u y v, construye expl´ıcitamente la matriz A.
- Escriba una funci´on [x, k, info] = GS3(d, u, v, b, x0, eps, maxit) que implementa el m´etodo de GS para resolver Ax = b, donde A es una matriz con la forma descrita anteriormente. Esta funci´on no construye expl´ıcitamente la matriz A. Los resultados parciales (puede agregar temporalmente disp o printf ) y finales deben ser los mismos por cualquiera de los dos formas siguientes: - utilizar GS. - con matriz3 construir expl´ıcitamente A y utilizar GS.
Es necesario adaptar el algortimo para tener en cuenta las especificidades de cada fila de A:
0.6 M´etodo de Muller
Este m´etodo sirve para hallar ra´ıces reales o complejas de polinomios reales p(x) = a 0 +a 1 x+a 2 x^2 +...+anxn. El polinomio p se puede expresar en funci´on de sus ra´ıces:
p(x) = an(x − r 1 )(x − r 2 ) · · · (x − rn).
Las ra´ıces complejas, no reales, siempre vienen por parejas, es decir si r = a + ib, b 6 = 0, es una ra´ız, entonces ¯r = a − ib, el conjugado de r, tambi´en es ra´ız. Para las ra´ıces complejas q(x) = (x − r)(x − ¯r) divide a p(x).
q(x) = (x − r)(x − ¯r) = (x − a − ib)(x − a + ib) = (x − a)^2 + b^2 = x^2 − 2 ax + (a^2 + b^2 ).
Si se halla una ra´ız real r entonces q(x) = (x − r) divide a p(x).
En general, si q(x) divide a p(x), entonces existe un polinomio s(x) tal que
p(x) = q(x)s(x), grado(p) = grado(q) + grado(s).
Entonces para sequir obteniendo las ra´ıces de p(x) basta con obtener las ra´ıces de s(x), polinomio m´as sencillo.
En el m´etodo de la secante, dados dos valores x 0 y x 1 se busca la recta que pasa por los puntos (x 0 , f (x 0 )), (x 1 , f (x 1 )); el siguiente valor x 2 est´a dado por el punto donde la recta corta el eje x. En el m´etodo de Muller, en lugar de una recta, se utiliza una par´abola. Dados tres valores x 0 , x 1 y x 2 , se construye la par´abola P (x) que pasa por los puntos (x 0 , f (x 0 )), (x 1 , f (x 1 )) y (x 2 , f (x 2 )); el siguiente valor x 3 est´a dado por el (un) punto tal que P (x 3 ) = 0.
La par´abola se puede escribir de la forma P (x) = a(x − x 2 )^2 + b(x − x 2 ) + c. Las f´ormulas que permiten calcular a, b y c son:
d = (x 0 − x 1 )(x 0 − x 2 )(x 1 − x 2 ),
a =
−(x 0 − x 2 )(f (x 1 ) − f (x 2 )) + (x 1 − x 2 )(f (x 0 ) − f (x 2 ) d
b =
(x 0 − x 2 )^2 (f (x 1 ) − f (x 2 )) − (x 1 − x 2 )^2 (f (x 0 ) − f (x 2 ) d
c = f (x 2 ).
Entonces
x 3 − x 2 = −b ±
b^2 − 4 ac 2 a Para reducir los errores de redondeo se “racionaliza” el numerador y buscando que el denominador resultante sea grande (en valor absoluto)
D = b^2 − 4 ac R =
D
δ =
b + R si |b + R| ≥ |b − R| b − R si |b + R| < |b − R|
x 3 = x 2 −
2 c δ
Si en una iteraci´on D = b^2 − 4 ac < 0 es necesario utilizar, a partir de ah´ı, aritm´etica compleja (Scilab lo hace autom´aticamente). Eso hace que los siguientes valores a, b y c no sean necesariamente reales. Muy posiblemente b^2 − 4 ac tampoco es real.
METODO DE MULLER PARA UNA RA´ ´IZ
datos: p, x 0 , x 1 , x 2 , εf , ε 0 , maxit f 0 = p(x 0 ), f 1 = p(x 1 ), f 2 = p(x 2 ) info= 0 para k = 1, ..., maxit si |f 2 | ≤ εf ent r = x 2 , info= 1, parar d = (x 0 − x 1 )(x 0 − x 2 )(x 1 − x 2 ) si |d| ≤ ε 0 ent parar calcular a, b y c seg´un (3) D = b^2 − 4 ac R =
D
δ 1 = b + R, δ 2 = b − R si |δ 1 | ≥ |δ 2 | ent δ = δ 1 , sino δ = δ 2 si |δ| ≤ ε 0 ent parar x 3 = x 2 − 2 c/δ x 0 = x 1 , x 1 = x 2 , x 2 = x 3 , f 0 = f 1 , f 1 = f 2 f 2 = p(x 2 ) fin-para k
Si el algoritmo anterior acaba normalmente, info valdr´a 1 y r ser´a una ra´ız, real o compleja. En la anterior descripci´on del algoritmo, |t| indica valor absoluto si t es real, e indica m´odulo si t es complejo.
- Averig¨ue c´omo efectuar divisi´on de polinomios en Scilab.
- Escriba una funci´on [r, info] = Muller1(coef, x0, x1, x2, epsf, eps0, maxit) que obtiene (si al final info vale 1) una ra´ız de un polinomio cuyos coeficientes reales est´an, en orden creciente de la potencia, en el vector coef. En el algoritmo presentado, esta funci´on va desde tomar o redefinir x 0 , x 1 , x 2 hasta sino q(x) = (x − r)(x − r¯).
- Escriba una funci´on [rr, info] = Muller(coef, x0, epsf, eps0, maxit) que obtiene (si al final info vale 1) todas las ra´ıces de un polinomio cuyos coeficientes reales est´an, en orden creciente de la potencia, en el vector coef. Las ra´ıces quedar´an en el vector rr. Esta funci´on utiliza varias veces Muller.
METODO DE MULLER´ datos: p, x 0 , εf , ε 0 , maxit r = x 0 , h = 0. 5 mientras grado(p) ≥ 3 x 0 = r, x 1 = x 0 + h, x 2 = x 1 + h (r, inf o) = M uller1(p, xo, x 1 , x 2 , εf , ε 0 , maxit) si inf o = 0, ent parar si |imag(r)| ≤ ε 0 ent q(x) = (x − r) sino q(x) = (x − r)(x − ¯r) p(x) = p(x)/q(x) fin-mientras calcular ra´ıces de p (de grado no superior a 2)
- Sean k 1 , k 2 , ..., k 7 los primeros siete d´ıgitos del n´umero de su documento de identidad. Halle las seis ra´ıces del polinomio cuyos coeficientes, en orden creciente de la potencia, son k 7 , k 6 , ..., k 1. Muestre los resultados intermedios para la primera utilizaci´on de Muller.
- La f´ormula para obtener p utilizando polinomios de Lagrange es
p(x) =
∑^ n
k=
ykLk(x).
En seudolenguaje,
p(x) ← 0 para k = 1, ..., n construir Lk(x) p(x) ← p(x) + yk Lk(x) fin-para
Escriba una funci´on p = polInterpPL(x, y) que obtiene el polinomio de interpolaci´on utilizando polinomios de Lagrange. Las coordenadas de los puntos est´an en los vectores columna x y y.
- Escriba una funci´on que calcula el factorial de un entero no negativo.
- Sean x , y dos vectores columna con las coordenadas de n puntos: ( 0, (−1)n−^1 (n − 1)! ), (1, 0), (2, 0), (3, 0), ..., (n − 1 , 0). Considere n = 5 y construya expl´ıcitamente los vectores columna x, y. Obtenga p, el polinomio de interpolaci´on, usando polInterpSE. Construya un vector columna yse con los valores del polinomio p evaluado en los valores del vector columna x. Te´oricamente los dos vectores, y y yse son iguales. Obtenga la norma de y-yse.
- Obtenga p, el polinomio de interpolaci´on, usando polInterpPL. Construya un vector columna ypl con los valores del polinomio evaluado en los valores del vector columna x. Obtenga la norma de y-ypl.
- Repita los dos pasos anteriores para n = 3, 4 , ..., 20. Construya una tabla de tres columnas, en la primera los valores de n, en la segunda columna los valores norm(y-yse) y en la tercera columna los valores norm(y-ypl). Concluya.
0.8 Error local y global del m´etodo del
trapecio y de Simpson
- Escriba una funci´on [a, b] = rectaMinC(X, Y) que obtiene los dos coeficientes de la recta de aproxi- maci´on por m´ınimos cuadrados y = ax + b para los puntos cuyas coordenadas est´an en los dos vectores columna X y Y.
m = size(X,1) A = [ones(m,1) X] aa = (A’A)(A’Y) b = aa(1) a = aa(2)
- Escriba una funci´on s = trapecio(f, a, b, n) que obtiene una aproximaci´on de la integral definida por medio de la f´ormula del trapecio
∫ (^) b
a
f (x) dx ≈ h
f (a) + f (b) 2
n∑− 1
i=
f (xi)
h = b − a n xi = a + ih
Las ´ordenes importantes son semejantes a:
h = (b-a)/n s = ( f(a) + f(b) )/ for i=1:n- s = s + f( a+ih) end s = sh
- Escriba una funci´on s = Simpson(f, a, b, n) que obtiene una aproximaci´on de la integral definida por medio de la f´ormula de Simpson (con n par)
∫ (^) b
a
f (x) dx ≈ h 3
f (a) + f (b) + 4
n∑− 1
i=1, 3 ,...
f (xi) + 2
n∑− 2
i=2, 4 ,...
f (xi)
h =
b − a n xi = a + ih
- Escoja una funci´on (diferente de un polinomio de grado inferior a 7) de la cual conoce la (una) antiderivada, es decir, para la cual se puede obtener f´acilmente el valor exacto de una integral definida. Por ejemplo f (x) = ex. Con a = 0, h = 0.2, b = a + h, n = 1, obtenga el valor aproximado por la f´ormula del trapecio, el valor exacto de
∫ (^) b a f^ (x)^ dx^ y^ e^ el valor absoluto del error cometido.^ Estos valores de h y e deben ser almacenados en arreglos H y E.
- Repita el paso anterior para h = 0. 21 , 0. 22 , ..., 0 .25 (de nuevo b = a + h, n = 1).
0.9 Comparaci´on de m´etodos de Simpson y Gauss-Legendre
Se desea comparar para una funci´on diferente de un polinomio, con antiderivada conocida (se puede obtener el valor exacto de una integral definida), el m´etodo de Simpson y el m´etodo de cuadratura de Gauss-Legendre, utilizando el mismo n´umero de evaluaciones de la funci´on.
- Escoja una funci´on, diferente de un polinomio, de la cual conoce la (una) antiderivada. Por ejem- plo f (x) = ex. Escoja un intervalo “razonable”, por ejemplo, a = 1, b = 2. Escriba una funci´on y = f08(x) donde se calcula la funci´on f.
- Escoja dos valores a < b y obtenga el valor exacto de
∫ (^) b
a
f (x) dx
- Obtenga una aproximaci´on de
∫ (^) b
a
f (x) dx utilizando intg.
- Obtenga una aproximaci´on de
∫ (^) b
a
f (x) dx utilizando la f´ormula de Simpson con 6 subintervalos (hay 7 evaluaciones de f ).
- Escriba una funci´on s = Simpson38(f, a, b, n) que obtiene una aproximaci´on de la integral definida por medio de la f´ormula de Simpson 3/8 con n m´ultiplo de 3. ∫ (^) x 3
x 0
f (x) dx ≈
h(y 0 + 3y 1 + 3y 2 + y 3 ).
Para el intervalo [a, b], con n (m´ultiplo de 3) subintervalos ∫ (^) b
a
f (x) dx ≈
h(y 0 + 3y 1 + 3y 2 + 2y 3 + 3y 4 + 3y 5 + 2y 6 + 3y 7 + · · · + 3yn− 2 + 3yn− 1 + yn).
Los pasos principales pueden ser:
h = (b − a)/n s ← f (a) + f (b)
s ← s + 3
n∑− 2
i=1, 4 ,...
f (xi)
s ← s + 3
n∑− 1
i=2, 5 ,...
f (xi)
s ← s + 2
n∑− 3
i=3, 6 ,...
f (xi)
s ← 3 hs/ 8
donde xi = a + ih. O tambi´en,
h = (b − a)/n s ← f (a) + f (b) para i = 1, ..., n − 1 y ← f (a + ih) si 3 |i ent s ← s + 2y sino s ← s + 3y fin-para s ← 3 hs/ 8
La notaci´on p|q significa “p divide a q”. En Scilab se puede escribir if modulo(q, p) == 0
- Obtenga una aproximaci´on de
∫ (^) b
a
f (x) dx mediante h = (b-a)/5, Simpson(f08, a, a+2h, 2) + Simpson38(f08, a+2h, b, 3) (se utilizan 6 evaluaciones de f ).
- En el m´etodo de cuadratura de Gauss-Legendre de orden m el c´alculo num´erico de la integral se hace mediante la siguiente aproximaci´on: ∫ (^1)
− 1
f (x) dx ≈ wm 1 f (xm 1 ) + wm 2 f (xm 2 ) + · · · wmmf (xmm).
Los valores wmj , los pesos o ponderaciones, y los valores xmj , las abcisas, son fijos y conocidos.
Cuadratura de Gauss-Legendre m ±xmj wmj 1 0.000000000000000 2. 2 0.577350269189626 1. 3 0.000000000000000 0. 0.774596669241483 0. 4 0.339981043584856 0. 0.861136311594053 0. 5 0.000000000000000 0. 0.538469310105683 0. 0.906179845938664 0. 6 0.238619186083197 0. 0.661209386466265 0. 0.932469514203152 0.
En Abramowitz y Stegun hay una tabla muy completa. Para integrar en otro intervalo es necesario hacer un cambio de variable, ∫ (^) b
a
f (x) dx ≈
b − a 2
∑^ m
j=
wmj f
b − a 2 xmj +
a + b 2
Escriba una funci´on s = GaussLeg1(f, a, b, m) que calcula de manera aproximada
∫ (^) b a f^ (x)^ dx^ por el m´etodo de cuadratura de Gauss-Legendre de orden m. En esta funci´on se definen y se utilizan dos matrices X y W. La matriz X es triangular inferior de tama˜no 6 × 6: en las primeras m posiciones de la fila m est´an las m abcisas y, de manera an´aloga, la matriz W tiene los pesos correspondientes. La parte principal de esta funci´on puede ser semejante a
X = [ 0 0 0 0 0 0; // -0.577350269189626 0.577350269189626 0 0 0 0 // 0 -0.774596669241483 0.774596669241483 0 0 0 // -0.339981043584856 0. .. -0.861136311594053 0.861136311594053 0 0 // 0 -0.538469310105683 0. ..
0.10 M´etodo de Newton en Rn
y c´alculo num´erico de la matriz jacobiana
- Considere una funci´on f : Rn^ → Rn, por ejemplo,
f (x 1 , x 2 ) = ( 2x^21 + 3x 1 x 2 + 4x 2 − 38 , − 5 x 1 x 2 + 6x 1 − 6 x 2 + 36 ) Escriba una funci´on fx = f09(x) que, para un vector columna x, construye el vector columna fx con las dos componentes de f (x).
- Encuentre, usando fsolve, un vector x tal que f (x) = 0 para la funci´on definida en f09.
- Dada una funci´on f : Rn^ → Rn, se denotar´a por f ′(x) la matriz jacobiana de f evaluada en x. En esta matriz n × n, la fila i est´a formada por las n derivadas parciales de fi (la i-´esima componente de f ). Escriba una funci´on J = jac09(x) que para un vector columna x, construye la matriz jacobiana J de la funci´on definida en f09.
- En el m´etodo de Newton en Rn^ (para hallar un x tal que f (x) = 0), en cada iteraci´on hay dos pasos fundamentales
resolver f ′(xk) dk^ = −f (xk) xk+1^ = xk^ + dk^.
Estos pasos presuponen que para caulquier xk^ se puede obtener f (xk) y f ′(xk). El m´etodo termina satisfactoriamente si
||f (xk)|| ≤ ε.
Las terminaciones no deseadas, pero posibles, son
- demasiadas iteraciones sin encontrar un x adecuado.
- f ′(xk) no es invertible (su determinante es nulo o casi nulo).
Escriba una funci´on [x, info, k] = Newton1(f, jac, x0, eps, maxit) que implementa el m´etodo de Newton. El resultado info valdr´a 1 si se obtuvo un x adecuado, 2 si en maxit iteraciones no se pudo obtener una soluci´on adecuada y valdr´a 0 si en alguna iteraci´on la matriz jacobiana es singular o casi singular. En la funci´on f se calcula f (x), en la funci´on jac se calcula f ′(x), el vector columna x0 es la aproximaci´on inicial. El n´umero de iteraciones realizadas estar´a en k.
- Verifique la funci´on anterior mediante [x, info, k] = Newton1(f09, jac09, [1; 1], 1.0e-4, 20) y con otros ejemplos.
- Averig¨ue mediante qu´e funci´on de Scilab puede obtener num´ericamente la matriz jacobiana de una funci´on. Escriba una funci´on [x, info, k] = Newton2(f, x0, eps, maxit) que implementa el m´etodo de Newton. Aqu´ı no es necesario tener una funci´on externa que eval´ua la matriz jacobiana, basta con tener f.
- Para una funci´on ϕ : R → R la derivada se puede aproximar por
ϕ′(x) ≈ ϕ(x + h) − ϕ(x) h y tambi´en ϕ′(x) ≈ ϕ(x + h) − ϕ(x − h) 2 h
La derivada parcial se puede aproximar por
∂fi ∂xj
(x) ≈ fi(x 1 , x 2 , ..., xj + h, ..., xn− 1 , xn) − fi(x 1 , ..., xn) h con h cercano a 0. Una mejor aproximaci´on es
∂fi ∂xj (x) ≈
fi(x 1 , x 2 , ..., xj + h, ..., xn− 1 , xn) − fi(x 1 , x 2 , ..., xj − h, ..., xn− 1 , xn) 2 h
Escriba una funci´on Dij = der_fi_xj(f, i, j, x) que calcula la anterior aproximaci´on.
- Escriba una funci´on J = jacobiana3(f, x) que calcula aproximadamente la matriz jacobiana.
- Escriba una funci´on [x, info, k] = Newton3(f, x0, eps, maxit) que implementa el m´etodo de Newton aproximando la matriz jacobiana mediante la funci´on anterior.
- Escriba una funci´on Dj = der_f_xj(f, j, x) que construye un vector columna con als n derivadas parciales con respecto a xj.
- Escriba una funci´on J = jacobiana4(f, x) que calcula aproximadamente la matriz jacobiana us- ando al funci´on anterior.
- Escriba una funci´on [x, info, k] = Newton4(f, x0, eps, maxit) que implementa el m´etodo de Newton aproximando la matriz jacobiana mediante la funci´on anterior.
- Con los mismos par´ametros de entrada, observe en n´umero de iteraciones para cada una de las 3 implementaciones.
- Una manera m´as precisa para tratar de medir la eficiencia, es mediante el n´umero de llamados a la funci´on f. Obtenga este n´umero para cada una de las implementaciones (puede utilizar global). Concluya.
