Metodos Numericos En Matlab 5e2749

UNIVERSIDAD JUÁREZ AUTÓNOMA DE TABASCO DIVISIÓN ACADÉMICA DE INGENIERÍA Y ARQUITECTURA Ingeniería Mecánica Eléctrica

MÉTODOS NUMÉRICOS

Profesor: Luis Adolfo Mejía Vázquez Alumna: Sonia Mariel García Díaz Grupo: I Matrícula: 102D18074

ÍNDICE

1. Introducción  Matlab  Error de aproximación  Análisis numérico 2. Solución de ecuaciones algebraicas y trascendentales  Método de aproximaciones sucesivas  Método de bisección  Método de Newton-Raphson 3. Solución de sistemas de ecuaciones lineales  Método de Gauss  Método de Jordan  Método de Jacobi  Método de Seidel 4. Integración numérica  Método de los trapecios  Método de Simpson 1/3 5. Solución numérica de ecuaciones diferenciales  Método de Euler  Método de Heun  Método de Nystrom  Método de Runge-Kutta 6. Aproximación funcional  Mínimos cuadrados

INTRODUCCIÓN

La enseñanza de los métodos numéricos, actualmente se encuentra íntimamente ligada a la construcción de sistemas matemáticos de alta precisión y exactitud. De ahí que ésta haya evolucionado sobre todo para satisfacer las necesidades y requerimientos de los instrumentos de control numérico, diseño de ingeniería, y en muchos otros de los campos científicos. La mayor parte de los ingenieros dependen de herramientas de tipo robótico, que requieren de aproximaciones cada vez más exactas y precisas. Es, por tanto, deseable que a la par de la enseñanza teórica de los métodos numéricos, vaya la enseñanza de un lenguaje de programación de alto nivel, que permita escribir directamente todos los algoritmos matemáticos que se van desarrollando durante el curso. En el caso del presente trabajo se utiliza el Matlab, por ser un sistema matemático, muy flexible y relativamente sencillo de utilizar. Matlab Algunas cosas que uno debe de saber de Matlab son:   

Las minúsculas y las mayúsculas no son equivalentes Un punto y coma al final del comando hará que no se vea en la pantalla el resultado Matlab utiliza paréntesis () y corchetes [], los cuales no son intercambiables

Algunos comandos que utilizaremos son: plot() crea un gráfico a partir de vectores y/o columnas de matrices, con escalas lineales sobre ambos ejes grid cuadrícula function nombre de la función o fichero *.m limits es un vector de 2 o 4 elementos, cuyos valores son [xmin, xmax] o [xmin, xmax, ymin, ymax] tol tolerancia de error relativo. Por defecto 2e-03 eval(‘cadena de caracteres’) evalúa como expresión de Matlab el texto entre apostrofes feval(nombre_función, calcular(a,b,c)

argumentos) feval(calcular, a, b, c) equivale a

inline transforma en función una cadena de caracteres input permite imprimir un mensaje en la línea de comandos de Matlab y recuperar como valor de retorno un valor numérico o el resultado de una expresión tecleada por el disp permite imprimir en pantalla un mensaje de texto o el valor de una matriz, pero sin imprimir su nombre. En realidad, disp siempre imprime vectores y/o matrices: las cadenas de caracteres son un caso particular de vectores norm norma de un vector zeros devuelve un vector o matriz de ceros Cómo introducir un vector de filas o de columnas: u = [2 4 5] tiene una fila con tres elementos (matriz de 1 x 3). v = [2; 4; 5] tiene tres filas separadas por puntos y comas (matriz de 3 x 1). v = [2 4 5]' o v = u' transpone u para generar la misma v. w = 2:5 define el vector de filas w = [2 3 4 5] mediante valores que aumentan sucesivamente en una unidad. u = 1:2:7 asigna valores que aumentan en dos unidades para obtener u = [1 3 5 7] Cómo definir una matriz (introduciendo las filas una por una): A = [1 2 3; 4 5 6] tiene dos filas (el punto y coma siempre separa unas filas de otras). B = [1 2 3; 4 5 6]' es la transpuesta de A. Así pues, AT es A' en MATLAB. Cómo generar matrices especiales diag(v) genera una matriz diagonal con el vector v como diagonal. toeplitz(v) define una matriz simétrica de diagonal constante con v como primera fila y primera columna. toeplitz(w, v) define una matriz simétrica de diagonal constante con w como primera columna y v como primera fila. ones(n) genera una matriz de n × n con todos los valores iguales a uno.

zeros(n) genera una matriz de n × n con todos los valores iguales a cero. eye(n) genera una matriz identidad de n × n. rand(n) genera una matriz de n × n con elementos de valor aleatorio entre 0 y 1 (distribución uniforme). randn(n) genera una matriz de n × n cuyos elementos siguen una distribución normal (media 0 y varianza 1). ones(m,n), zeros(m,n), rand(m,n) generan matrices de m × n. ones(size(A)), zeros(size(A)), eye(size(A)) generan matrices de la misma forma que A. Error de aproximación El error de aproximación o error numérico es una medida del ajuste o cálculo de una magnitud con respecto al valor real o teórico que dicha magnitud tiene. Un aspecto importante de los errores de aproximación es su estabilidad numérica. Dicha estabilidad se refiere a cómo dentro de un algoritmo de análisis numérico el error de aproximación es propagado dentro del propio algoritmo. El concepto de error es consustancial con el cálculo numérico. En todos los problemas es fundamental hacer un seguimiento de los errores cometidos a fin de poder estimar el grado de aproximación de la solución que se obtiene. Los errores asociados a todo cálculo numérico tienen su origen en dos grandes factores:  

Aquellos que son inherentes a la formulación del problema. Los que son consecuencia del método empleado para encontrar la solución del problema.

Dentro del grupo de los primeros, se incluyen aquellos en los que la definición matemática del problema es sólo una aproximación a la situación física real. Estos errores son normalmente despreciables; por ejemplo, el que se comete al obviar los efectos relativistas en la solución de un problema de mecánica clásica. En aquellos casos en que estos errores no son realmente despreciables, nuestra solución será poco precisa independientemente de la precisión empleada para encontrar las soluciones numéricas. Otra fuente de este tipo de errores tiene su origen en la imprecisión de los datos físicos: constantes físicas y datos empíricos. En el caso de errores en la medida de los datos empíricos y teniendo en cuenta su carácter generalmente aleatorio,

su tratamiento analítico es especialmente complejo pero imprescindible para contrastar el resultado obtenido computacionalmente. En lo que se refiere al segundo tipo de error (error computacional), tres son sus fuentes principales: 1. Equivocaciones en la realización de las operaciones (errores de bulto). Esta fuente de error es bien conocida por cualquiera que haya realizado cálculos manualmente o empleando una calculadora. El empleo de computadores ha reducido enormemente la probabilidad de que este tipo de errores se produzcan. Sin embargo, no es despreciable la probabilidad de que el programador cometa uno de estos errores (calculando correctamente el resultado erróneo). Más aún, la presencia de bugs no detectados en el compilador o en el software del sistema no es inusual. Cuando no resulta posible verificar que la solución calculada es razonablemente correcta, la probabilidad de que se haya cometido un error de bulto no puede ser ignorada. Sin embargo, no es esta la fuente de error que más nos va a preocupar. 2. El error causado por resolver el problema no como se ha formulado, sino mediante algún tipo de aproximación. Generalmente está causado por la sustitución de un infinito (sumatorio o integración) o un infinitesimal (diferenciación) por una aproximación finita. Algunos ejemplos son:  El cálculo de una función elemental (por ejemplo, Seno x) empleando sólo n términos de los infinitos que constituyen la expansión en serie de Taylor.  Aproximación de la integral de una función por una suma finita de los valores de la función, como la empleada en la regla del trapezoide.  Resolución de una ecuación diferencial reemplazando las derivadas por una aproximación (diferencias finitas).  Solución de la ecuación f(x) = 0 por el método de Newton-Raphson: proceso iterativo que, en general, converge sólo cuando el número de iteraciones tiende a infinito. Denominaremos a este error, en todas sus formas, como error por truncamiento, ya que resulta de truncar un proceso infinito para obtener un proceso finito. Obviamente, estamos interesados en estimar, o al menos acotar, este error en cualquier procedimiento numérico. 3. Por último, la otra fuente de error de importancia es aquella que tiene su origen en el hecho de que los cálculos aritméticos no pueden realizarse con precisión ilimitada. Muchos números requieren infinitos decimales para ser representados correctamente, sin embargo, para operar con ellos es necesario redondearlos. Incluso en el caso en que un número pueda

representarse exactamente, algunas operaciones aritméticas pueden dar lugar a la aparición de errores (las divisiones pueden producir números que deben ser redondeados y las multiplicaciones dar lugar a más dígitos de los que se pueden almacenar). El error que se introduce al redondear un número se denomina error de redondeo. Análisis numérico El análisis numérico o cálculo numérico es la rama de las matemáticas que se encarga de diseñar algoritmos para, a través de números y reglas matemáticas simples, simular procesos matemáticos más complejos aplicados a procesos del mundo real. El análisis numérico proporcionará todo el andamiaje necesario para llevar a cabo todos aquellos procedimientos matemáticos susceptibles de expresarse algorítmicamente, basándose en algoritmos que permitan su simulación o cálculo en procesos más sencillos empleando números. Definido el error, junto con el error isible, pasamos al concepto de estabilidad de los algoritmos. Muchas de las operaciones matemáticas pueden llevarse adelante a través de la generación de una serie de números que a su vez alimentan de nuevo el algoritmo (). Esto proporciona un poder de cálculo y refinamiento importantísimo a la máquina que a medida que va completando un ciclo va llegando a la solución. El problema ocurre en determinar hasta cuándo deberá continuar con el ciclo, o si nos estamos alejando de la solución del problema. Finalmente, otro concepto paralelo al análisis numérico es el de la representación, tanto de los números como de otros conceptos matemáticos como los vectores, polinomios, etc. Por ejemplo, para la representación en ordenadores de números reales, se emplea el concepto de coma flotante que dista mucho del empleado por la matemática convencional. En general, estos métodos se aplican cuando se necesita un valor numérico como solución a un problema matemático, y los procedimientos "exactos" o "analíticos" (manipulaciones algebraicas, teoría de ecuaciones diferenciales, métodos de integración, etc.) son incapaces de dar una respuesta. Debido a ello, son procedimientos de uso frecuente por físicos e ingenieros, y cuyo desarrollo se ha visto favorecido por la necesidad de éstos de obtener soluciones, aunque la precisión no sea completa.

SOLUCIÓN DE ECUACIONES ALGEBRAICAS Y TRASCENDENTALES

Una ecuación es una expresión matemática, que contiene términos algebraicos, funciones trigonométricas, funciones exponenciales, logarítmicas, etc., números y un signo de igual. Ejemplo: 3x2+6x-2=1 Existen diferentes tipos de ecuaciones, dependiendo de los términos, las hay lineales como por ejemplo: 45y+26=73. En general se puede hablar de ecuaciones polinomiales de la forma:

Las ecuaciones trigonométricas son las que tienen funciones trigonométricas, ejemplo: sen(2x)+5x=2. Las ecuaciones logarítmicas son las que tienen logaritmos, como por ejemplo: 5ln(3x)+10=1. En general se distinguen dos tipos de ecuaciones, las lineales y las no lineales. Una raíz de una ecuación es un número que satisface le ecuación, es decir que si esta igualada a cero, una raíz será el número real o complejo que la hace cero. Ejemplo: la raíz de la ecuación 5x-1=0 será 1/5 Una ecuación lineal siempre tendrá como raíz un solo número real. Las ecuaciones no lineales pueden tener normalmente más de una raíz. Las raíces de las ecuaciones no lineales pueden ser reales o complejas. En este momento es importante mencionar que una de las partes importantes de los métodos numéricos es el cálculo aproximado de las raíces de ecuaciones que no se pueden obtener por métodos algebraicos. Los métodos tienen la característica de que bajo ciertas condiciones, se aproximan a las raíces de las ecuaciones, con un cierto error de aproximación, hay también que considerar que hay ciertos métodos que tienen algunas restricciones que dificultan su uso. La solución de una ecuación es el número, o conjunto de números reales o complejos que, que al substituirlo en lugar de la variable satisface la igualdad Ejemplo la solución de la ecuación:

3x3+6x2+7x=0 es {-1+1.15i, -1-1.15i, 0} Como se podrá observar a partir de este momento, por cada raíz compleja, que es solución de la ecuación, se tiene su conjugado. Hay que observar que el número total de raíces que tiene un polinomio, corresponde al grado del término mayor del polinomio. En el ejemplo anterior, el grado del polinomio es tres, y tiene 3 raíces.

Método de aproximaciones sucesivas El método de aproximaciones sucesivas representa la esencia de los procesos iterativos ya que permite definir una ecuación de recurrencia que, en apariencia, no tiene sentido desde el punto de vista algebraico, pero que resulta muy atinada si se toma un valor inicial y se mejora a través de las iteraciones. Sin embargo, el método como tal no es cien por ciento aplicable para cualquier ecuación algebraica o trascendente; no obstante, se utiliza como base para completar otros métodos abiertos. Definición del método Aproximaciones sucesivas es un método abierto, es decir, no necesita de un intervalo que atrape a una raíz, sino que requiere de un valor x 0 que representa un valor aproximado a la raíz y de cuya cercanía al valor real dependerá la velocidad en que se cumpla con una tolerancia preestablecida. Una forma sencilla de definir un método de aproximaciones sucesivas consiste en despejar de una ecuación a la variable independiente; esto se aplica particularmente en ecuaciones que por su forma no permiten despejar fácilmente a su incógnita. Por ejemplo, en la ecuación x2+7x-ex=0 no puede lograrse un despeje sencillo, algebraicamente hablando. Desde un punto de vista iterativo, la ecuación puede expresarse como:

En efecto, algebraicamente hablando, el despeje anterior no aporta mejora en la solución de la ecuación. Sin embargo, sí se define en forma iterativa:

Donde xi es un valor inicial y xi+1 es un valor corregido que, en un escenario favorable, tendrá una cantidad de error menor con respecto a la raíz de la ecuación. El proceso iterativo se detendrá cuando entre dos aproximaciones sucesivas se satisfaga la tolerancia preestablecida.

La ecuación que representa el método de aproximaciones sucesivas expresada en forma iterativa es:

La aplicación del método consiste en proporcionar una aproximación inicial a la raíz de la ecuación (que puede obtenerse por medios gráficos o al detectar un cambio de signo en la función tabular) y sustituirla en la ecuación, obteniéndose una nueva aproximación. De nuevo deberá sustituirse la última aproximación hasta que la diferencia entre dos aproximaciones sucesivas satisfaga determinada tolerancia preestablecida. La principal aportación del método de aproximaciones sucesivas es la determinación del criterio de convergencia, ya que se debe aplicar a todos los métodos abiertos o de punto fijo. Una manera de explicar los pasos para llegar a la raíz de una ecuación usando este método sería la siguiente. Primero, tenemos nuestra ecuación f(x). Proponemos un valor inicial x 1 y lo sustituimos en la ecuación para encontrar y1. Después le sumamos a x1 un Δx y obtenemos x2. Luego calculamos y2 sustituyendo x2 en la ecuación. Ahora multiplicamos y1 por y2, si el resultado es positivo entonces quiere decir que entre x1 y x2 no hay raíz por lo tanto x2 se vuelve x1 y y2 se vuelve y1 y volvemos a empezar sumando el Δx a x1 para encontrar nuestro siguiente x2. Cuando la multiplicación de y1 y y2 nos dé negativo entonces sabemos que en ese intervalo hay una raíz y procedemos a dividir Δx entre 10. Y seguimos con el proceso iterativo hasta que nuestro Δx sea menor o igual a la tolerancia.

Descripción gráfica

Algoritmo del método  

Grado de la ecuación G Coeficientes de la ecuación ai 0≤i≤G

               

x1←0 Dx←1 Tolerancia←0.01 Contador←0 y1←F(x1) Repetir x2←x1+Dx y2←F(x2) Si (y1*y2)>0 Entonces x1←x2 y1←y2 De lo contrario Dx←Dx/10 Contador←Contador+1 Hasta que (Dx<=Tolerancia) o (Contador>100)

Programación en Matlab Nombre del programa: aproximaciones_sucesivas.m f=input('ingrese la funcion, entre apostrofos: '); x1=-10; %Valor inicial dx=1; tol=input('ingrese el valor de la tolerancia '); contador=0; f=inline(f); y1=feval(f,x1); while dx>tol contador=contador+1; x2=x1+dx; y2=feval(f,x2); z=y1.*y2; if z>0 x1=x2; y1=y2; else dx=dx/10; end end fprintf(' La raiz es: %1.8f\n Encontrada en %1.0f iteraciones\n ',x1,contador); Ejemplo en Matlab con la ecuación 2+3x-5x2+8x3 Los valores de las raíces son {-0.3494168904, 0.4872+0.6914i, 0.4872-0.6914i} Se corre el programa con diferentes valores de tolerancia y vemos que entre más pequeña sea la tolerancia más nos acercamos al valor real. Pero este método es muy inexacto y necesita tener una tolerancia demasiado pequeña para que nos dé el resultado más aproximado. >> aproximaciones_sucesivas ingrese la funcion, entre apostrofos: '2+3*x-5*x.^2+8*x.^3' ingrese el valor de la tolerancia 0.01 La raiz es: -0.40000000 Encontrada en 17 iteraciones >> aproximaciones_sucesivas ingrese la funcion, entre apostrofos: '2+3*x-5*x.^2+8*x.^3' ingrese el valor de la tolerancia 0.0001 La raiz es: -0.35000000 Encontrada en 24 iteraciones

>> aproximaciones_sucesivas ingrese la funcion, entre apostrofos: '2+3*x-5*x.^2+8*x.^3' ingrese el valor de la tolerancia 0.0000001 La raiz es: -0.34941690 Encontrada en 45 iteraciones >> Gráfica de la ecuación 2+3x-5x2+8x3 Nombre del programa: grafica.m f=input('ingrese la funcion, entre apostrofos: '); x=input('ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: '); f=inline(f); y=feval(f,x); close plot(x,y) grid Corremos el programa y nos muestra la gráfica de la ecuación donde podemos ver el intervalo en que se encuentra la raíz. >> grafica ingrese la funcion, entre apostrofos: '2+3*x-5*x.^2+8*x.^3' ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: [-2:1:10] >>

Método de bisección El método de bisección se aplica a funciones algebraicas o trascendentes y proporciona únicamente raíces reales. Tiene su origen en un popular algoritmo de búsqueda de datos en arreglos vectoriales denominado búsqueda binaria. Es un método cerrado, es decir, requiere de un intervalo en el cual esté atrapada una raíz. Básicamente, consiste en cortar el intervalo en dos justo por la mitad (bisectar) considerando a este punto como una aproximación de la raíz de la función. Posteriormente, debe determinarse si la raíz verdadera se encuentra a la derecha o a la izquierda de la aproximación y, según corresponda, cerrar el intervalo con la aproximación y el limite derecho izquierdo, pero siempre manteniendo a la raíz verdadera en el intervalo. Esta operación se repite hasta que la diferencia entre las dos últimas aproximaciones sea menor que una tolerancia preestablecida. Bisección es un método robusto, aunque resulta lento en su proceso por lo oneroso de los cálculos que deben realizarse; por otra parte, su convergencia puede en ocasiones ser inestable. Definición del método A partir de una función algebraica o trascendente y de un intervalo [ a,b] que pertenece al dominio de la función y para el cual f(a)*f(b)<0, lo que implica que en el intervalo [a,b] existe al menos una raíz. El método consiste en bisectar el intervalo [a,b]:

Obteniendo una aproximación a la raíz x0; la función se evalúa en este nuevo valor y de acuerdo al signo de la función valuada en este punto, deberá sustituirse uno de los extremos del intervalo de búsqueda, de tal forma que se conserve que f(a)*f(b)<0. De acuerdo a la geometría de la figura, la sustitución de los intervalos deberá hacerse de la siguiente forma. Sea a tal que f(a)>0 y b tal que f(b)<0: 

Si f(x0)>0, entonces x0 sustituye a a



Si f(x0)<0, entonces x0 sustituye a b

En cada iteración deberá sustituirse alguno de los límites del intervalo que contiene a la raíz. Repitiendo este proceso, el intervalo se reduce paulatinamente hasta que alguna de las aproximaciones coincide razonablemente con la raíz de la función. El proceso se detiene cuando, entre la aproximación xi y la aproximación anterior xi-1 se satisface un nivel de error (absoluto o relativo) preestablecido (tolerancia). Descripción gráfica

Algoritmo del método  

Grado de la ecuación G Coeficientes de la ecuación ci 0≤i≤G



Intervalo de la raíz [a,b]

      

Valor de la tolerancia x0←0 xf←0 Contador←0 Repetir Contador←Contador+1 x0←(a+b)/2

       

y0←F(x0) Si y0<0 entonces a←x0 De lo contrario b←x0 error←|(x0-xf)| x0←xf Hasta que error >= tolerancia

Programación en Matlab Nombre del programa: biseccion.m f=input('ingrese la funcion, entre apostrofos: '); x=input('ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: '); f=inline(f); y=feval(f,x); close plot(x,y) grid Conv = 1; while Conv >= 0 f=input('ingrese la funcion, entre apostrofos: '); x=input('ingrese el intervalo en que se encuentra la raiz, en la forma [a b]: '); f=inline(f); a=x(1); b=x(2); fa=feval(f,a); fb=feval(f,b); Conv = fa.*fb; end tol = input('\nTeclee el valor de la tolerancia tol = '); Error = tol+1; xo = 0; xf = 0; Cont = 0; %Comienzo del algoritmo iter=1; while Error >= tol Cont = Cont +1; xo=(a+b).*0.5; fxo=feval(f,xo); if fxo < 0 a=xo; else b=xo; end Error = abs(xo-xf); xf = xo; end fprintf(' La raiz es: %1.8f\n Encontrada en %1.0f iteraciones\n ',xo,Cont);

Ejemplo en Matlab con la ecuación 3x3+3x2-x-1 Los valores de las raíces son {0.5773502692, -0.5773502692, -1} Se corre el programa con un intervalo aleatorio, y hasta que damos el intervalo correcto, el programa continúa y nos da la raíz. Nótese que en este programa ya se incluye la gráfica de la ecuación y con ésta podemos ver en qué intervalo se encuentran las raíces, en vez hacerlo al tanteo. El programa se corre con distintos valores de tolerancia y distintos intervalos. >> biseccion ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: [0:1:10] ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese el intervalo en que se encuentra la raiz, en la forma [a b]: [0 10] Teclee el valor de la tolerancia tol = 0.01 La raiz es: 0.57617188 Encontrada en 10 iteraciones >> >> biseccion ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: [0:1:10] ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese el intervalo en que se encuentra la raiz, en la forma [a b]: [-5 3] Teclee el valor de la tolerancia tol = 0.0001 La raiz es: -1.00006104 Encontrada en 17 iteraciones >> >> biseccion ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: [0:1:10] ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese el intervalo en que se encuentra la raiz, en la forma [a b]: [-1 5] ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' ingrese el intervalo en que se encuentra la raiz, en la forma [a b]: [-2 0] ingrese la funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1'

ingrese el [a b]: [-1 ingrese la ingrese el [a b]: [-3

intervalo en que se encuentra la raiz, en la forma 0] funcion, entre apostrofos: '3*x.^3+3*x.^2-x-1' intervalo en que se encuentra la raiz, en la forma 8]

Teclee el valor de la tolerancia tol = 0.0000001 La raiz es: 0.57735030 Encontrada en 27 iteraciones >> Gráfica de la ecuación 3x3+3x2-x-1 Como el programa biseccion.m ya incluye la gráfica de la función, automáticamente nos grafica la ecuación que le introducimos con el intervalo que nosotros escogemos. Esta es la gráfica de la ecuación del ejemplo donde podemos observar que efectivamente hay una raíz en el intervalo [0,1].

Método de Newton-Raphson El método Newton-Raphson es, junto con bisección uno de los más populares. Su preferencia radica en su robustez y velocidad al encontrar la raíz. Por ser un método abierto o de punto fijo, debe verificarse su criterio de equivalencia. Quizá el único punto que pudiera tener en contra radica en la necesidad de contar con la primera y segunda derivadas de la ecuación a resolver. Se aplica a ecuaciones algebraicas y trascendentes y proporciona raíces reales. Definición del método El nombre original del método es de las tangentes. Una tangente es una recta que intersecta a una curva en un solo punto; en consecuencia, es perpendicular a su radio. A partir de la figura 1, se plantea que en un valor x 0 que represente una aproximación a la raíz de la ecuación, se trace una tangente en el punto f(x0). Esta recta tangente deberá cortar al eje horizontal. El punto donde lo corte será la nueva aproximación x1, de tal forma que en el punto f(x1) se trace una nueva tangente. Este proceso se repetirá hasta que el corte de la tangente en el eje horizontal coincida con la raíz de la ecuación, o bien, cuando la diferencia entre dos aproximaciones sucesivas cumpla con una tolerancia preestablecida. La ecuación del método de Newton-Raphson es la siguiente:

Descripción gráfica

Algoritmo del método  

Grado de la ecuación G Coeficientes de la ecuación ai 0≤i≤G

             

Valor de la tolerancia x1←0 Contador←0 Repetir y1←F(x1) dy1←F’(x1) x2←x1-(y1/dy1) Si |x1-x2|>tolerancia Entonces x1←x2 Contador←Contador+1 Hasta que (|x1-x2|<=tolerancia) o (contador>100) Si contador<100 Entonces imprimir ‘x2’ De lo contrario imprimir ‘no se encontró la raíz’

Programación en Matlab Nombre del programa: newton_raphson.m f=input('ingrese la funcion a graficar, entre apostrofos: '); p=input('ingrese un intervalo en el que se graficara la ecuacion,\nen la forma [a:incremento:b]: '); f=inline(f); s=feval(f,p); close plot(p,s) grid syms x y=input('ingrese la funcion sin apostrofos F(x) = '); %funcion a derivar dy1=diff(y,1) %derivada de primer orden dy2=diff(y,2) %derivada de segundo orden z=(y.*dy2)./dy1.^2; %Criterio de convergencia disp(z) z=inline(z); x=input('ingrese un valor de x = '); fz=feval(z,x); disp(fz) if abs(fz) < 1 tol = input('\nTeclee el valor de la tolerancia tol = '); Error = tol+1; Cont = 0; xi=x; while Error >= tol Cont = Cont +1; xf=xi-y./dy1; disp(xf) xf=inline(xf) fxf=feval(xf,xi); Error = abs(fxf-xi); disp(Error) xi = fxf; end fprintf(' \nLa raiz es %1.8f,\nEncontrada en %1.0f Iteraciones ',fxf,Cont); else fprintf(' \nLa ecuacion no converge '); end

Ejemplo en Matlab con la ecuación x3-7x2+14x-6 Los valores de las raíces son {0.5857864376, 3.414213562, 3} >> newton_raphson ingrese la funcion a graficar, entre apostrofos: 'x.^37*x.^2+14*x-6' ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: [0:1:10] ingrese la funcion sin apostrofos F(x) = x.^3-7*x.^2+14*x-6 ingrese un valor de x = 0 Teclee el valor de la tolerancia tol = 0.0001 La raiz es 0.58578644, Encontrada en 5 Iteraciones >> Gráfica de la ecuación x3-7x2+14x-6

Ejemplo en Matlab con la ecuación e-x-x >> newton_raphson ingrese la funcion a graficar, entre apostrofos: 'exp(-x)-x' ingrese un intervalo en el que se graficara la ecuacion, en la forma [a:incremento:b]: [-4:1:4] ingrese la funcion sin apostrofos F(x) = exp(-x)-x ingrese un valor de x = 0 Teclee el valor de la tolerancia tol = 0.00001 La raiz es 0.56714329, Encontrada en 4 Iteraciones >> Gráfica de la ecuación e-x-x

SOLUCIÓN DE SISTEMAS DE ECUACIONES LINEALES

En general, un sistema con m ecuaciones lineales y n incógnitas puede ser escrito en forma normal como:

Donde son las incógnitas y los números son los coeficientes del sistema. Es posible reescribir el sistema separando con coeficientes con notación matricial:

Si representamos cada matriz con una única letra obtenemos:

Donde A es una matriz m por n, x es un vector columna de longitud n y b es otro vector columna de longitud m. El sistema de eliminación de Gauss-Jordan se aplica a este tipo de sistemas, sea cual sea el cuerpo del que provengan los coeficientes. La matriz A se llama matriz de coeficientes de este sistema lineal. A b se le llama vector de términos independientes del sistema y a x se le llama vector de incógnitas. Lo que se pretende hacer es desarrollar métodos numéricos y sus algoritmos, que permitan resolver los sistemas lineales que se sabe tienen solución, lo cual ha sido establecido, por un análisis con las herramientas del álgebra lineal.

Haciendo un resumen de los sistemas lineales y las características, para que estos tengan solución, se puede enunciar: Sea Ax=b un sistema lineal de ecuaciones  Cuando un sistema tiene al menos una solución se dice que, el sistema es consistente, en caso contrario es inconsistente.  Se dice que el sistema es homogéneo cuando b=0, es decir el vector b es el vector nulo. Estos sistemas siempre tienen solución, por lo tanto estos sistemas son consistentes. Las solución trivial de estos sistemas es también X=0  Si b≠0 se dice que el sistema es no homogéneo.  Un sistema no homogéneo es consistente (tiene al menos una solución ) si rango(A) = rango(A,b)  Los sistemas no homogéneos son consistentes con una solución si rango(A)=n Es consistente con más de una solución si rango(A)=rango(A,b)= r rango(A)

Método de Gauss A través de las operaciones fundamentales aplicadas a la matriz de coeficientes A y al vector de términos independientes b (para mantener la igualdad en las ecuaciones), se convierte a la matriz A en una matriz triangular superior donde la diagonal principal la constituyen números. Esta transformación implica que en la ecuación n aparezca una sola incógnita, dos en la ecuación n-1, tres en la n-2 y así consecutivamente para después realizar una sustitución hacia atrás. Definición del método Sea el sistema de ecuaciones lineales Ax=b, la matriz ampliada del sistema es:

Por medio de las operaciones fundamentales debe transformarse en una matriz ampliada triangular superior de la forma:

Donde los valores a’ij y b’i son los coeficientes modificados por la aplicación de las operaciones fundamentales. Este procedimiento se logra tomando cada uno de los elementos de la diagonal principal de la matriz A original, al cual se le denomina pivote y después normalizando la ecuación con dicho valor. Posteriormente, a través de las operaciones fundamentales, se hace la eliminación de los elementos inferiores de la columna correspondiente al pivote. Una vez realizada la transformación se realiza la sustitución hacia atrás obteniéndose el valor de cada una de las incógnitas.

Algoritmo del método  

Número de ecuaciones NE Valores de la matriz aumentada ai,I 1≤i≤NE 1≤I≤NE+1

   

Hacer I de 1 hasta NE-1 Hacer R de I+1 hasta NE Hacer C de I+1 hasta NE+1 aR,C←aR,C-(aI,C/aI,I)*aR,I



xNE←aNE,NE+1/aNE,NE

   

Hacer I de NE-1 hasta 1 Suma←0 Hacer K de I+1 hasta NE Suma←Suma+aI,K*xK



xI←(1/aI,I)*(aI,NE+1-Suma)

Programación en Matlab El siguiente programa está hecho para trabajar como una función de Matlab, por lo que antes de usarlo tenemos que introducir la matriz de coeficientes A y la matriz de términos independientes b. Nombre de la función: gauss.m function [x]=gauss(A,b) % Solución de Ax=b mediante eliminación de Gauss n=size(A,1); x=zeros(n,1); for i=1:n-1 % Transformación matriz A [p maxk]=max(abs(A(i:n,i))); maxk=maxk+i-1; if i~=maxk A([i maxk],:)=A([maxk i],:); b([i maxk]) =b([maxk i]); end j=i+1:n; A(j,i)=A(j,i)/A(i,i); A(j,j)=A(j,j)-A(j,i)*A(i,j); b(j)=b(j)-b(i)*A(j,i); end for i=n:-1:1 % Sustitución inversa x(i)=(b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i); end

Ejemplo en Matlab con el sistema de ecuaciones:

Las soluciones del sistema son: x1=9/13, x2=-2/13, x3=-10/13 >> A=[3 2 1;2 2 4;-1 0.5 -1] A = 3.0000 2.0000 -1.0000

2.0000 2.0000 0.5000

>> b=[1;-2;0] b = 1 -2 0 >> [x]=gauss(A,b) x = 0.6923 -0.1538 -0.7692 >>

1.0000 4.0000 -1.0000

Método de Jordan Es una ampliación del método de Gauss con la diferencia que la matriz de coeficientes A se transforma en la matriz identidad de tal forma que se conserva una única incógnita por ecuación eliminando la sustitución hacia atrás. Definición del método Sea el sistema de ecuaciones lineales Ax=b, por medio de las operaciones fundamentales debe transformarse en la matriz identidad de la forma:

Donde los valores

b’i

son los coeficientes modificados por la aplicación de las

operaciones fundamentales. De nuevo, el procedimiento utilizado es la elección del pivote, la normalización de las ecuaciones y la eliminación de los elementos superiores e inferiores de la columna del pivote correspondiente.

Algoritmo del método  

Número de ecuaciones NE Valores de la matriz aumentada ai,k 1≤i≤NE 1≤k≤NE+1

    

Hacer I de 1 hasta NE Hacer R de 1 hasta NE Si R<>I entonces Hacer C de I+1 hasta NE+1 aR,C←aR,C-(aI,C/aI,I)*aR,I

 

Hacer C de I+1 hasta NE+1 aI,C←aI,C/aI,I

 

Hacer I de 1 hasta NE xI←aI,NE+1

Programación en Matlab El siguiente programa también está hecho para trabajar como una función de Matlab, por lo que antes de usarlo tenemos que introducir la matriz aumentada a. Nombre de la función: Jordan.m function [X] = Jordan(a) [n,t] = size(a); for k = 1:n for j = n + 1:-1:k a(k, j) = a(k, j) / a(k, k); end for i = 1:n if i ~= k for j = n + 1:-1:k a(i, j) = a(i, j) - a(i, k) * a(k, j); end end end end for i=1:n X(i) = a(i,n+1); end end

Ejemplo en Matlab con el sistema de ecuaciones:

Las soluciones del sistema son: x=1, y=1, z=1 >> a=[5 -3 -1 1;1 4 -6 -1;2 3 4 9] a = 5 1 2

-3 4 3

-1 -6 4

>> [X]=Jordan(a) X = 1 >>

1

1

1 -1 9

Método de Jácobi Los métodos de Jácobi y de Seidel son los equivalentes en la solución de sistemas de ecuaciones lineales al método de aproximaciones sucesivas en la solución de ecuaciones algebraicas y trascendentes. Consiste básicamente en obtener una ecuación de recurrencia (matricial en este caso) y proponer un vector solución inicial; posteriormente, se deberán realizar las iteraciones necesarias hasta que la diferencia entre dos vectores consecutivos cumpla con una tolerancia predefinida. En realidad estos métodos representan una adaptación vectorial de un proceso escalar, lo que implica la necesidad de adaptar los conceptos necesarios: los procesos iterativos se detienen cuando entre dos aproximaciones consecutivas se cumple con determinado error preestablecido. En este caso, deberá medirse la norma entre dos vectores para reconocer el momento en que se satisface la cota de error. Por otra parte, resta el hecho de tener que evaluar un criterio de equivalencia el cual, naturalmente, tendrá carácter vectorial. Definición del método Sea el sistema de ecuaciones lineales AX=b, donde A es la matriz de coeficientes, X es el vector de incógnitas y b el vector de términos independientes. En esta ecuación se puede sustituir a la matriz A por la suma de dos matrices: A=D+R. En donde la matriz D es una matriz cuyos elementos son cero excepto los elementos de la diagonal que corresponden a los elementos de la matriz A y R que es una matriz con ceros en la diagonal y sus restantes elementos coinciden con los respectivos de A.

Al sustituir D y R en la ecuación del sistema se llega a la siguiente ecuación:

Esta ecuación no aporta una solución por sí misma, si se observa desde la óptica del algebra matricial. Sin embargo, si se aplica desde una forma recursiva:

Para k=0,1,2,…,n y donde X(k) representa un vector solución inicial y X(k+1) representa una aproximación posterior a la inicial X(k). Se puede constatar claramente que la ecuación obtenida es totalmente representativa de un método de aproximaciones sucesivas. Esta ecuación requiere de un breve análisis para su aplicación práctica. En principio, la matriz D, solo posee elementos diferentes de cero (que corresponden a los propios de A) en su diagonal principal. Es fácilmente comprobable que la matriz inversa D-1 también posee únicamente valores diferentes de 0 en su diagonal principal y que estos valores corresponden a los recíprocos de sus valores en la matriz A, es decir, serán 1/ aii. Por otra parte, el resto de los

elementos de cada renglón de la matriz A se encuentran en la matriz R y son restados del vector de términos independientes. En contexto la ecuación equivale, a partir del sistema de ecuaciones lineales, a despejar a la incógnita de ubicada en la diagonal principal de cada una de las ecuaciones que conforman el sistema, de la siguiente forma:

El método Jácobi propone que el vector inicial X0 sea igual a cero. A partir de esta propuesta, el vector siguiente será X1=bi/aii, es decir, el elemento independiente entre el coeficiente de la diagonal principal para cada ecuación.

Este vector X(1) se sustituye en las ecuaciones de los despejes anteriores obteniéndose el siguiente vector X(2). El proceso se realiza consecutivamente hasta que la norma entre dos vectores consecutivos es menor que cierta tolerancia preestablecida. La norma θ se calcula como:

El método de Jácobi es susceptible de los efectos del pivoteo. En consecuencia, su criterio de convergencia lo conforman los criterios de la diagonal pesada, mismo que posee dos condiciones: 1. Condición necesaria: Es condición necesaria que el elemento ubicado en la diagonal principal de cada ecuación sea mayor en valor absoluto que el resto de los elementos de la misma ecuación. 2. Condición suficiente: Es condición suficiente que el elemento ubicado en la diagonal principal de cada ecuación sea mayor en valor absoluto que la suma del resto de los elementos de la misma ecuación.

Algoritmo del método  

Número de ecuaciones NE Valores de la matriz aumentada ai,k 1≤i≤NE 1≤k≤NE+1

         

Valor de la tolerancia Contador←0 Hacer I de 1 hasta NE xI←0 Repetir Hacer I de 1 hasta NE Suma←0 Hacer K de 1 hasta NE Si K<>I entonces Suma←Suma+aI,K*xK



zI←(1/aI,I)*(aI,NE+1-Suma)

           

Hacer I de 1 hasta NE Bandera←cierto Si |zI-xI|>tolerancia entonces Bandera←falso Contador←Contador+1 Hacer I de 1 hasta NE xI←zI Hasta que (Bandera=cierto) o (Contador>100) Si Bandera=cierto entonces Imprimir xI De lo contrario Imprimir ‘No se encontraron las soluciones’

Programación en Matlab El siguiente programa también está hecho para trabajar como una función de Matlab, por lo que antes de usarlo tenemos que introducir la matriz de coeficientes A, la matriz de términos independientes b, el vector inicial x, la tolerancia tol y el número máximo de iteraciones maxit. Nombre de la función: jacobi.m function [x,iter]=jacobi(A,b,x,tol,maxit) N=diag(diag(A)); P=N-A; corr=1; errest=1; iter=0; while abs(errest)>tol & iter<maxit iter=iter+1; x0=x; corr0=corr; x=N\(P*x0+b); corr=norm(x-x0,inf); normest=corr/corr0; if normest>=1 & iter>=2 error('norma de la matriz de iteración > 1') end errest=normest/(1-normest)*corr; end iter

Ejemplo en Matlab con el sistema de ecuaciones:

Las soluciones del sistema son: x1=-1, x2=3, x3=5 >> A=[10 1 2;4 6 -1;-2 3 8] A = 10 4 -2

1 6 3

2 -1 8

>> b=[3;9;51] b = 3 9 51 >> x=[0;0;0] x = 0 0 0 >> tol=0.00001; >> maxit=100; >> [x,iter]=jacobi(A,b,x,tol,maxit) x = -1.0000 3.0000 5.0000 iter = 15 >>

Método de Seidel Este método es una versión acelerada del método de Jácobi. En el método de Jácobi es necesario contar con un vector aproximado completo para proceder a la sustitución en las ecuaciones de recurrencia y obtener una nueva aproximación. Definición del método En el método de Seidel se propone ir sustituyendo los nuevos valores de la aproximación siguiente conforme se vayan obteniendo sin esperar a tener un vector completo. De esta forma se acelera la convergencia. A partir de las ecuaciones de recurrencia del método de Jácobi:

El criterio de convergencia del método de Seidel corresponde totalmente al criterio de la diagonal pesada cuyas condiciones ya se expresaron anteriormente.

Algoritmo del método  

Número de ecuaciones NE Valores de los coeficientes ai,k 1≤i≤NE 1≤k≤NE+1

          

Valor de la tolerancia Contador←0 Hacer I de 1 hasta NE xI←0 Repetir Bandera←cierto Hacer I de 1 hasta NE Suma←0 Hacer K de 1 hasta NE Si K<>I entonces Suma←Suma+aI,K*xK



b←(1/aI,I)*(aI,NE+1-Suma)

       

Si |b-xI|>tolerancia Entonces Bandera←falso xI←b Contador←Contador+1 Hasta que (Bandera=cierto) o (Contador>100) Si Bandera=cierto Entonces imprimir xI De lo contrario imprimir ‘No se encontraron soluciones’

Programación en Matlab El siguiente programa también está hecho para trabajar como una función de Matlab, por lo que antes de usarlo tenemos que introducir la matriz de coeficientes A, la matriz de términos independientes b, el vector inicial x, la tolerancia tol y el número máximo de iteraciones maxit. Nombre de la función: seidel.m function [x,iter]=seidel(A,b,x,tol,maxit) N=tril(A); P=N-A; corr=1; errest=1; iter=0; while abs(errest)>tol & iter<maxit iter=iter+1; x0=x; corr0=corr; x=N\(P*x0+b); corr=norm(x-x0,inf); normest=corr/corr0; if normest>=1 & iter>=2 error('norma de la matriz de iteración > 1') end errest=normest/(1-normest)*corr; end iter

Ejemplo en Matlab con el sistema de ecuaciones: A manera de comparación usaremos el mismo ejemplo que usamos en el método de Jácobi para ver la diferencia entre los dos métodos.

Las soluciones del sistema son: x1=-1, x2=3, x3=5

>> A=[10 1 2;4 6 -1;-2 3 8] A = 10 4 -2

1 6 3

2 -1 8

>> b=[3;9;51] b = 3 9 51 >> x=[0;0;0] x = 0 0 0 >> tol=0.00001; >> maxit=100; >> [x,iter]=seidel(A,b,x,tol,maxit) x = -1.0000 3.0000 5.0000 iter = 6 >> Como pudimos observar, el método de Jácobi resolvió el sistema de ecuaciones en 15 iteraciones, mientras que con el método de Seidel se resolvió en solo 6 iteraciones. Por lo que vemos que realmente el método de Seidel es una versión acelerada del método de Jácobi.

INTEGRACIÓN NUMÉRICA

En la práctica profesional, es común encontrarse con problemas que requieren del conocimiento del área bajo una curva de una función que representa el comportamiento de un proceso, que depende de un conjunto de variables de control; las funciones, no siempre se encuentran en forma analítica, muchas veces se reportan en forma tabular, o varían drásticamente entre intervalos, por este motivo, no siempre es posible integrarlas por métodos analíticos, y se requiere de un tipo de método apropiado para cada caso. La integración definida es el procedimiento sustentado por el teorema fundamental del cálculo, por medio del cual se puede calcular de forma exacta una sumatoria de Riemman. La notación de Leibnitz permite escribir la integración indefinida de una función como:

Los métodos de integración numérica se basan prácticamente en usar una forma aproximada (obtenida por interpolación) de una función y aplicarle métodos simples para obtener un valor numérico como aproximación al valor real de la integral. Desde el punto de vista numérico, existen dos tipos generales de métodos de integración, estos son los métodos de integración cerrada, y los métodos de integración abierta. La integración numérica es un procedimiento algebraico y aritmético, que permite el cálculo de integrales definidas, que o bien no se pueden resolver analíticamente o su solución es muy complicada. Los métodos cerrados de integración, son aquellos que tienen un alcance exclusivamente dentro del intervalo de definición de la función (intervalo en el cual se reportan los valores de la función). A estos se les conoce como los métodos de integración cerrada de Newton-Cotes.

Dentro de los métodos cerrados de Newton-Cotes se tienen:  La regla del trapecio  La regla de Simpson 1/3 Los métodos de Newton-Cotes funcionan en general en dos pasos: 1. Se divide el intervalo [a,b] en n intervalos de igual amplitud, cuyos valores extremos son sucesivamente:

Para quedar en la nueva notación x0=a y xn=b 2. Se aproxima f(x) por un polinomio de grado n, P n(x) y se integra, para obtener la aproximación de la integral.

Método de los trapecios Este es un método cerrado, parte de los métodos de Newton-Cotes y tiene como objetivo hallar la integral aproximada de una función a través de un polinomio de primer grado, es decir uniendo mediante una recta los puntos en donde se evaluará la función. Definición del método La regla del trapecio consiste en: el intervalo original queda igual, la función se aproxima por una línea recta, y el área bajo la curva corresponde al área de un trapecio. De acuerdo a los cursos de geometría, el área de un trapecio, es igual a la suma de la base mayor, más la base menor, por la altura sobre dos. Así que cada trapecio que se forma dentro de una curva tiene como área:

Al sumar todas las áreas se tiene que:

Considerando los incrementos constantes la expresión anterior queda:

Descripción gráfica

Algoritmo del método  

Grado de la ecuación G Valores de los coeficientes ai 0≤i≤G

           

Límite inferior LI Límite superior LS Número de segmentos NS Dx←(LS-LI)/NS y1←F(LI) yNS+1←F(LS) x←LI Suma←0 Hacer I de 2 hasta NS x←x+Dx Suma←Suma+F(x) A←(Dx/2)*(y1+2*Suma+yNS+1)

Programación en Matlab Nombre del programa: trapecios.m clear all; clc; fprintf('Calculo de la integral por el metodo trapecial\n\n'); f=input('introduce la funcion:','s'); a=input('limite inferior:'); b=input('limite superior:'); c=input('numero de segmentos a dividir:'); h=(b-a)/c; z=0; for x=a:h:b k=eval(f); if x==a,d=k; end if x==b,e=k; end z=z+k; end z=z-d-e; z=z*2; z=z+d+e; z=z/(2*c); z=z*(b-a) fprintf('Resultado ');

Ejemplo en Matlab con la integral:

Nota: el ángulo es en radianes. Valor real del área: 8.183479208 Calculo de la integral por el metodo trapecial introduce la funcion:2+sin(2*sqrt(x)) limite inferior:1 limite superior:6 numero de segmentos a dividir:10000 z = 8.1835 Resultado >> Gráfica de la función 2+sen(2√x)

Método de Simpson 1/3 Este es otro método cerrado, parte de los métodos de Newton-Cotes. Tiene condición que el numero n de puntos que conforman la función tabular sea par. Definición del método La regla de Simpson consiste en: dividir el intervalo [a,b] en 2M sub intervalos [xk,xk+1] de la misma anchura:

mediante una partición de nodos equi espaciados xk=a+kh, y con esta

que corresponde a la integral:

y se le llama regla de Simpson 1/3.

Descripción gráfica

Algoritmo del método  

Grado de la ecuación G Valores de los coeficientes ai 0≤i≤G

                 

Límite inferior LI Límite superior LS Número de segmentos NS par Dx←(LS-LI)/NS y1←F(LI) yNS+1←F(LS) N←NS/2 PAR←0 x←LI-Dx Hacer I de 1 hasta N x←x+2*Dx PAR←PAR+F(x) IMPAR←0 x←LI Hacer I de 1 hasta N-1 x←x+2*Dx IMPAR←IMPAR+F(x) A←(Dx/3)*(y1+4*PAR+2*IMPAR+yNS+1)

Programación en Matlab clear all; clc; fprintf('Calculo de la integral por el metodo de Simpson de 1/3\n\n'); f=input('introduce la funcion:','s'); a=input('limite inferior:'); b=input('limite superior:'); c=input('numero de segmentos a dividir (numero par):'); h=(b-a)/c; z=0; x=a; for i=1:c; if (-1)^i==1 k=eval(f); z=z+k; end x=h*i; end zz=0; x=a; for i=2:c; if (-1)^i==-1 k=eval(f); zz=zz+k; end x=h*i; end x=a; if x==a d=eval(f); end x=b; if x==b e=eval(f); end z=z*4; v=zz*2; z=z+v+d+e; z=z/(3*c); z=z*(b-a) fprintf('Resultado ');

Ejemplo en Matlab con la integral:

Valor real del área: 1.718281828 Calculo de la integral por el metodo de Simpson de 1/3 introduce la funcion:exp(x) limite inferior:0 limite superior:1 numero de segmentos a dividir (numero par):10 z = 1.7183 Resultado >> Gráfica de la función

SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

Una ecuación diferencial es una ecuación en la que intervienen derivadas de una o más funciones desconocidas. Dependiendo del número de variables independientes respecto de las que se deriva, las ecuaciones diferenciales se dividen en: Ecuaciones diferenciales ordinarias: aquellas que contienen derivadas respecto a una sola variable independiente. Ecuaciones en derivadas parciales: aquellas que contienen derivadas respecto a dos o más variables. Una ecuación diferencial es una ecuación que incluye expresiones o términos que involucran a una función matemática incógnita y sus derivadas. Algunos ejemplos de ecuaciones diferenciales son:

Es una ecuación diferencial ordinaria, donde especificada de la variable independiente derivada de con respecto a .

representa una función no

, es decir

,

es la

La expresión:

Es una ecuación en derivadas parciales. A la variable dependiente también se le llama función incógnita (desconocida). La resolución de ecuaciones diferenciales es un tipo de problema matemático que consiste en buscar una función que cumpla una determinada ecuación diferencial. Se puede llevar a cabo mediante un método específico para la ecuación diferencial en cuestión o mediante una transformada. Ecuación diferencial lineal Se dice que una ecuación es lineal si tiene la forma

Es decir:   

Ni la función ni sus derivadas están elevadas a ninguna potencia distinta de uno o cero. En cada coeficiente que aparece multiplicándolas sólo interviene la variable independiente. Una combinación lineal de sus soluciones es también solución de la ecuación.

Ecuaciones semilineales y cuasilineales No existe un procedimiento general para resolver ecuaciones diferenciales no lineales. Sin embargo, algunos casos particulares de no linealidad sí pueden ser resueltos. Son de interés el caso semilineal y el caso cuasilineal. Una ecuación diferencial ordinaria de orden n se llama cuasilineal si es "lineal" en la derivada de orden n. Más específicamente, si la ecuación diferencial ordinaria para la función puede escribirse en la forma:

Se dice que dicha ecuación es cuasilineal si .

es una función afín, es decir

Una ecuación diferencial ordinaria de orden n se llama semilineal si puede escribirse como suma de una función "lineal" de la derivada de orden n más una función cualquiera del resto de derivadas. Formalmente, si la ecuación diferencial ordinaria para la función puede escribirse en la forma:

Se dice que dicha ecuación es semilineal si

es una función lineal.

Solución de una ecuación diferencial Tipos de soluciones Una solución de una ecuación diferencial es una función que al reemplazar a la función incógnita, en cada caso con las derivaciones correspondientes, verifica la ecuación, es decir, la convierte en una identidad. Hay tres tipos de soluciones:  Solución general: una solución de tipo genérico, expresada con una o más constantes. Es un haz de curvas. Tiene un orden de infinitud de acuerdo a su cantidad de constantes (una constante corresponde a una familia simplemente infinita, dos

constantes a una familia doblemente infinita, etc). En caso de que la ecuación sea lineal, la solución general se logra como combinación lineal de las soluciones (tantas como el orden de la ecuación) de la ecuación homogénea (que resulta de hacer el término no dependiente de ni de sus derivadas igual a 0) más una solución particular de la ecuación completa.

 Solución particular: Si fijando cualquier punto por donde debe pasar necesariamente la solución de la ecuación diferencial, existe un único valor de C, y por lo tanto de la curva integral que satisface la ecuación, éste recibirá el nombre de solución particular de la ecuación en el punto , que recibe el nombre de condición inicial. Es un caso particular de la solución general, en donde la constante (o constantes) recibe un valor específico.

 Solución singular: una función que verifica la ecuación, pero que no se obtiene particularizando la solución general. Solución de la ecuación no consistente en una particular de la general.

Método de Euler Las herramientas numéricas para la solución de ecuaciones diferenciales aportan una aproximación muy buena en comparación con los resultados obtenidos en forma analítica. Son varios los métodos numéricos que ofrecen la solución de estas ecuaciones diferenciales. Los métodos paso a paso se fundamentan en los polinomios interpolantes que provienen de tablas de diferencias; en consecuencia será necesario disponer de un intervalo de solución de la variable independiente cuyo primer punto debe coincidir con la condición inicial. El método de Euler es un método predictivo-correctivo o de paso a paso ya que resuelve la ecuación en intervalos entre dos pares de puntos y repitiendo esta solución en todos los puntos que conforman a la función tabular. Definición del método El planteamiento inicial es resolver una ecuación como la siguiente:

Integramos ambos de la ecuación entre dos puntos consecutivos.

De la que se obtiene:

Resolviendo la integral a partir del polinomio interpolante de Newton-Gregory obtenemos la siguiente ecuación:

A esta ecuación se le denomina método de Euler para resolver ecuaciones diferenciales de primer orden con condiciones iniciales en un intervalo equiespaciado. Descripción gráfica

Algoritmo del método            

Valor inicial de yn Valor inicial de to Valor inicial de tu Número de incrementos NI Dt←( tu- to)/NI tn← to Hacer I de 1 hasta NI dyn←f(yn,tn) yn+1←yn+dyn*Dt Imprimir yn, Dt, dyn, yn+1 yn←yn+1 tn←tn+Dt

Programación en Matlab Nombre del programa: euler.m fprintf('\n \tRESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO METODO DE EULER\n') f=input('\nIngrese la ecuacion diferencial de la forma: dy/dx=f(x,y)\n','s'); x0=input('\nIngrese el primer punto x0:\n'); x1=input('\nIngrese el segundo punto x1:\n'); y0=input('\nIngrese la condicion inicial y(x0):\n'); n=input('\nIngrese el numero de pasos n:\n'); h=(x1-x0)/n; xs=x0:h:x1; y1=y0; fprintf('\n''it x0 x1 y1'); for i=1:n it=i-1; x0=xs(i); x=x0; x1=xs(i+1); y=y0; y1=y0+h*eval(f); fprintf('\n%2.0f%10.6f%10.6f%10.6f\n',it,x0,x1,y1); y0=y1; end fprintf('\n El punto aproximado y(x1) es = %10.6f\n',y1);

Ejemplo en Matlab

>> euler RESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO METODO DE EULER Ingrese la ecuacion diferencial de la forma: dy/dx=f(x,y) sqrt(x^2+y^2) Ingrese el primer punto x0: 2 Ingrese el segundo punto x1: 2.3 Ingrese la condicion inicial y(x0): 0.5 Ingrese el numero de pasos n: 3 'it x0 x1 y1 0 2.000000

2.100000

0.706155

1

2.100000

2.200000

0.927710

2

2.200000

2.300000

1.166470

El punto aproximado y(x1) es = >>

1.166470

Método de Heun La mejora del método de Heun consiste en la aproximación a la pendiente mediante la aplicación de dos derivadas del intervalo, una en el punto inicial y otra en el punto final. La aproximación mejorada de la pendiente será el promedio de las dos derivadas. Definición del método Recordando el método de Euler, la pendiente al principio del intervalo es:

Que se emplea en extrapolar linealmente el valor de ‘y’ al final del intervalo:

Este valor no es la solución final sino una predicción intermedia, por lo que se ha distinguido a ésta con el superíndice 0. Esta ecuación se denomina ecuación predictora y proporciona una aproximación del valor de ‘y’ al final del intervalo. Este valor nos permite a su vez calcular la pendiente aproximada en dicho punto.

Combinando las dos pendientes obtenemos el promedio del intervalo:

Ésta pendiente promedio se usa para extrapolar linealmente el valor de la función en el siguiente punto usando el método de Euler:

Esta ecuación se denomina ecuación correctora. El método de Heun sigue un esquema predictor-corrector que es el mismo de los métodos de pasos múltiples y se expresa como:

Descripción gráfica

Algoritmo del método             

Valor inicial de yn Valor inicial de to Valor inicial de tu Número de incrementos NI Dt←( tu- to)/NI tn← to Hacer I de 1 hasta NI dyn←f(yn,tn) yn+1←yn+dyn*Dt dyn+1←yn+((dyn+dyn+1)2)*Dt Imprimir yn, tn, dyn, dyn+1, yn+1 yn←yn+1 tn←tn+Dt

Programación en Matlab El siguiente programa está hecho para trabajar como una función de Matlab, por lo que antes de usarlo tenemos que introducir la función, los valores iniciales y el número de incrementos. Nombre de la función: heun.m function [tvals,yvals]=heun(f,t0,b,y0,h) n=floor((b-t0)/h)+1; m=length(y0); k1=zeros(1,m); tvals=zeros(n+1,1); yvals=zeros(n+1,m); index=[0:1:n]'; tvals=t0+h*index; yvals(1,:)=y0; h2=h/2; for i=2:n+1 k1=feval(f,tvals(i-1), yvals(i-1,:)); yvals(i,:)=yvals(i-1,:)+h*k1; yvals(i,:)=yvals(i-1,:)+h2*(k1+feval(f,tvals(i),yvals(i,:)); end Ejemplo en Matlab Consideremos aquí las ecuaciones diferenciales que se obtienen de las leyes de Newton aplicadas al problema de dos cuerpos. Suponemos aquí que uno de los cuerpos es mucho más masivo que el otro de modo que su movimiento es descartable, ej., la tierra y un satélite. Suponemos también que el movimiento es en un plano. Como la fuerza gravitacional es inversamente proporcional a la distancia entre los cuerpos, tenemos tomando todas las constantes envueltas como uno, que la posición (x(t),y(t)) del cuerpo pequeño está dada por el sistema de ecuaciones diferenciales:

Tomamos como condiciones iniciales

Debido a que este es un sistema de orden dos, tenemos que hacer la sustitución:

Lo cual transforma el sistema de arriba al siguiente sistema de orden uno:

Definimos ahora la siguiente subrutina en MATLAB que evalúa el lado derecho de este sistema: function f=satelite(t,u) f=zeros(1,4); denom=(u(1)^2+u(3)^2)^1.5; f(1)=u(2); f(2)=-u(1)/denom; f(3)=u(4); f(4)=-u(3)/denom; Ahora calculamos y graficamos la solución del problema de valor inicial con la siguiente secuencia de instrucciones en MATLAB. Note que graficamos el conjunto de puntos (x(t),y(t)) para los t's generados en lugar de (t,x(t)) y (t,y(t)). Tenemos pues: >> >> >> >> >>

[t,y]=heun('satelite',0,10,[0.4,0,0.1,2],0.01); plot(y(:,1),y(:,3)) xlabel('X');ylabel('Y'); title('Solucion particular del problema de dos cuerpos')

Resultado:

Note que la curva es efectivamente una elipse aunque esto se acentúa en la gráfica por que los ejes tienen unidades de largo distintas.

Método de Nystrom En este método, que es de tipo predictor-corrector, se cuenta con una fórmula para determinar el valor de yi+1 aproximado, y otra para corregirlo hasta satisfacer una condición de error. El método de Nystrom es de segundo orden, por lo que al utilizarlo para resolver una ecuación diferencial, es necesario contar con los primeros dos pares de datos x,y. Definición del método La fórmula de predicción es:

La fórmula de corrección es:

Una vez satisfecha la condición de error, se calcula el siguiente valor de yi+1 con la fórmula de predicción y ese valor se trata con la fórmula de corrección (nuevamente hasta que la condición de error se satisfaga), y de esa forma se calculan los siguientes valores de yi+1, hasta que el valor de y sea yf.

Algoritmo del método                

Valor inicial de yn-1 Valor inicial de to Valor inicial de tu Número de intervalos NI Dt←( tu- to)/NI tn-1← to Hacer I de 1 hasta NI dyn+1←f(yn-1,tn-1) yn←yn-1+dyn+1*Dt dyn←f(yn,tn-1+Dt) yn+1←yn-1+2*dyn*Dt Imprimir yn, tn, dyn, dyn+1, yn+1 yn-1←yn tn-1←tn-1+Dt yn←yn+1 tn←tn-1+Dt

Programación en Matlab % Metodo de Nystrom % Uso: [t,y,dy]=nystrom(funcion,inicio,fin,alfa,paso) function [ts,ys,fs]=nystrom(f,inicio,fin,alfa,h) if nargin<5 disp('Faltan argumentos de entrada'); return; else ts=inicio:h:fin; ys(:,1)=alfa; n=(fin-inicio)/h; for i=1:n k1=feval(f,ts(i),ys(:,i)); k2=feval(f,ts(i)+h*(2/3),ys(:,i)+h*(2/3)*k1); k3=feval(f,ts(i)+h*(2/3),ys(:,i)+h*(2/3)*k2); fs(:,i)=k1; ys(:,i+1)=ys(:,i)+(2*k1+3*k2+3*k3)*h/8; end; end;

Ejemplo del método Obténgase la solución de la siguiente ecuación diferencial por el método de Nystrom, si x0=0, y0=2, h=0.1, xf=0.5 y e=1×10-5.

Recuérdese que para resolver ecuaciones diferenciales por este método es necesario conocer los valores de x1 y y1. Estos valores se han determinado a partir de los datos de x0, y0 y h por métodos de un solo paso.

Estos valores nos permiten conocer (por medio de las fórmulas predictora y correctora), los valores de yi+1 siguientes, resultando:

Método de Runge-Kutta El procedimiento de Runge-Kutta se aplica de forma directa a una ecuación diferencial de primer orden, pero también se extiende a un sistema de ecuaciones de primer orden, a una ecuación diferencial de segundo orden y a un sistema de ecuaciones diferenciales de segundo orden. El procedimiento de Runge-Kutta se puede programar fácilmente en los ordenadores y además, se emplea mucho en la práctica, debido a su exactitud relativamente elevada de la solución aproximada de la ecuación diferencial. Definición del método Sea una ecuación diferencial de primer orden, con la condición inicial de que en el instante t0 el valor inicial de x es x0. Se elige una anchura de paso h y se calculan cuatro números k1, k2, k3, k4 de acuerdo con el procedimiento esquematizado en la tabla adjunta. Según el procedimiento ordinario de Runge-Kutta, a partir del valor de x en el instante t se determina el valor de x en el instante t+h mediante la última fórmula en la tabla.

Algoritmo del método                

Valor inicial de yn Valor inicial de to Valor inicial de tu Número de intervalos NI Dt←( tu- to)/NI yn←yo tn← to Hacer I de 1 hasta NI k1←f(yn,tn) k2←f(yn+k1/2,tn+Dt/2) k3←f(yn+k2/2,tn+Dt/2) k4←f(yn+k3,tn+Dt) yn+1←yn+(Dt/6)*(k1+2*k2+2*k3+k4) Imprimir (tn, yn, k1, k2, k3, k4, yn+1) yn←yNI tn←tn+Dt

Programación en Matlab Definimos la función rk_1.m que resuelve la ecuación diferencial de primer orden, cuando le pasamos la función f(t,x), la condición inicial de que en el instante t0 el valor inicial es x0, el instante final tf y el número de pasos de integración n comprendidos entre el instante inicial t0 y el final tf. Nos devolverá un vector t y su correspondiente vector x. function [t,x] =rk_1(f,t0,tf,x0,n) h=(tf-t0)/n; t=t0:h:tf; x=zeros(n+1,1); %reserva memoria para n elementos del vector x x(1)=x0; for i=1:n k1=h*f(t(i),x(i)); k2=h*f(t(i)+h/2,x(i)+k1/2); k3=h*f(t(i)+h/2,x(i)+k2/2); k4=h*f(t(i)+h,x(i)+k3); x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6; end end

Ejemplo en Matlab Considérese el circuito en serie de la figura. Inicialmente el condensador está descargado. Si se cierra el interruptor I la carga empieza a fluir produciendo corriente en el circuito, el condensador se empieza a cargar. Una vez que el condensador adquiere la carga máxima, la corriente cesa en el circuito.

Escribimos el script carga.m para que establezca mediante comandos input: la resistencia R del circuito, la capacidad C del condensador, el tiempo final tf, el número de pasos n. También fijará las condiciones iniciales, en las que en el instante inicial t=0, el condensador está descargado x=0. Definirá la función f(t,x), Llamará al procedimiento numérico rk_1 y mediante el comando plot realizará una representación gráfica de la solución numérica. Después realizará una representación gráfica de la solución exacta.

Ejemplo: R=2.0, C=0.8, y tf=10. V0=10; R=input('Resistencia R: '); C=input('Capacidad C: '); tf=input('tiempo final, tf: '); n=input('número de pasos, n: '); f=@(t,x) V0/R-x/(R*C); %condiciones iniciales t0=0; x0=0; [t,x]=rk_1(f,t0,tf,x0,n); hold on plot(t,x,'b') y=C*V0*(1-exp(-t/(R*C))); plot(t,y,'r') xlabel('t') ylabel('q'); legend('aproximada','exacta','Location','Southeast') title('Carga del condensador') hold off

En la ventana de comandos corremos el script carga.m >> carga Resistencia R: 2 Capacidad C: 0.8 tiempo final, tf: 10 número de pasos, n: 50 >>

No se aprecia la diferencia entre la solución exacta y la numérica, aplicando el procedimiento de Runge-Kutta debido a su alta exactitud.

APROXIMACIÓN FUNCIONAL

Una función p(x) definida en un intervalo [a,b] aproxima a una función f(x) en el mismo intervalo si: p(x) ≈ f(x) Para describir un método de aproximación funcional, son necesarios dos ingredientes fundamentales. En primer lugar, hay que definir el espacio de funciones donde se elige el aproximante p(x). Es decir, hay que elegir el tipo de aproximación. En segundo lugar, es necesario formalizar matemáticamente el significado de “es aproximadamente igual a” ≈. Es decir, hay que definir el criterio de aproximación. La aproximación funcional se utiliza en dos situaciones típicas: con la función dato f(x) conocida a partir de su expresión analítica, o con f(x) conocida de forma discreta, únicamente en unos puntos base xi (i=0,…,n). Si se dispone de la expresión de f(x), el objetivo es aproximarla por una función p(x) más manejable y computable, que habitualmente nos permitirá además aproximar integrales o derivadas de la función original. Si sólo se dispone de valores discretos f(xi), además del objetivo reseñado, nos interesa obtener p(x) para poder evaluarlo en puntos distintos a los puntos base. Las funciones discretas son comunes en muchas áreas, desde el diseño asistido por ordenador (CAD) al tratamiento de señales, pasando por las medidas experimentales.

Mínimos cuadrados Mínimos cuadrados es una técnica de análisis numérico enmarcada dentro de la optimización matemática, en la que, dados un conjunto de pares ordenados: variable independiente, variable dependiente, y una familia de funciones, se intenta encontrar la función continua, dentro de dicha familia, que mejor se aproxime a los datos (un "mejor ajuste"), de acuerdo con el criterio de mínimo error cuadrático. En su forma más simple, intenta minimizar la suma de cuadrados de las diferencias en las ordenadas (llamadas residuos) entre los puntos generados por la función elegida y los correspondientes valores en los datos. Específicamente, se llama mínimos cuadrados promedio (LMS) cuando el número de datos medidos es 1 y se usa el método de descenso por gradiente para minimizar el residuo cuadrado. Se puede demostrar que LMS minimiza el residuo cuadrado esperado, con el mínimo de operaciones (por iteración), pero requiere un gran número de iteraciones para converger. Desde un punto de vista estadístico, un requisito implícito para que funcione el método de mínimos cuadrados es que los errores de cada medida estén distribuidos de forma aleatoria. También es importante que los datos a procesar estén bien escogidos, para que permitan visibilidad en las variables que han de ser resueltas. La técnica de mínimos cuadrados se usa comúnmente en el ajuste de curvas. Muchos otros problemas de optimización pueden expresarse también en forma de mínimos cuadrados. Definición del método La aproximación mínimo cuadrática consiste en minimizar el error cuadrático mencionado más arriba, y tiene solución general cuando se trata de un problema de aproximación lineal (lineal en sus coeficientes ) cualesquiera que sean las funciones base: antes mencionadas. Por lineal se entiende que la aproximación buscada se expresa como una combinación lineal de dichas funciones base. Sea

un conjunto de n pares con abscisas distintas, y sea

un conjunto de m funciones linealmente independientes (en un espacio vectorial de funciones), que se llamarán funciones base. Se desea encontrar una función de dicho espacio, o sea, combinación lineal de las funciones base, tomando por ello la forma:

.

Ello equivale por tanto a hallar los m coeficientes:

. En concreto, se

desea que tal función sea la mejor aproximación a los n pares empleando, como criterio de "mejor", el criterio del mínimo error cuadrático medio de la función

con respecto a los puntos

.

El error cuadrático medio será para tal caso:

Minimizar el error cuadrático medio es equivalente a minimizar el error cuadrático, definido como el radicando del error cuadrático medio, esto es:

Así, los que minimizan también minimizan derivando e igualando a cero este último:

, y podrán ser calculados

Siendo i=1,2, . . .,m Se obtiene un sistema de m ecuaciones con m incógnitas.

, para i=1, 2,..., m

, para i=1, 2,..., m Si se desarrolla la suma, se visualiza la ecuación "i-ésima" del sistema de m ecuaciones normales:

Para cada i=1, 2,..., m Lo cual, en forma matricial, se expresa como:

Siendo y g(x) como:

el producto escalar discreto, definido para dos funciones dadas h(x)

, Y para una función h(x) y vector cualquiera u, como:

La resolución de dicho sistema permite obtener, para cualquier base de funciones derivables localmente, la función f(x) que sea mejor aproximación mínimo cuadrática al conjunto de puntos antes mencionado. La solución es óptima –esto es, proporciona la mejor aproximación siguiendo el criterio de mínimo error cuadrático–, puesto que se obtiene al optimizar el problema. Si se tratara de hallar el conjunto de coeficientes exactamente por todos los pares a

tal que

, esto es, tales que

pase interpole

, entonces tendría que cumplirse que:

Que en forma matricial se expresa como:

Esto establece un sistema de n ecuaciones y m incógnitas, y como en general n>m, quedaría sobre determinado: no tendría siempre una solución general. Por tanto, la aproximación tratará en realidad de hallar el vector c que mejor aproxime .

Descripción gráfica

Ejemplo: función lineal Programación en Matlab function correl xi=input('ingrese los valores de xi:\n'); yi=input('ingreso los valores de yi:\n'); xm=input('ingreso un valor de x para analisar la pendiente:\n'); n=length(xi); p=xi.'; q=yi.'; r=p.*q; s=p.*p; a1=sum(xi); a2=sum(yi); a3=sum(r); a4=sum(s); m=a2/n; j=sum((yi-m).^2); A=[a4 a1;a1 n]; B=[a3 a2]; X=inv(A)*(B'); l=X(1);r=X(2); fprintf('\nLos coeficientes de la ecuacion son:\n') fprintf(' a b\n') fprintf('%8.8f%10.8f\n',l,r) fun=input('\nLa funcion lineal es y=','s'); y=l*xi+r; t=sum((yi-y).^2); r2=1-(t/j); z=(r2)^0.5; m=diff(fun); x=xm; n=eval(m); fprintf('\n El coeficiente de determinacion r^2='); fprintf('%5.5f\n',r2); if(n>0) fprintf('\n El coeficiente de correlacion r='); fprintf('%5.5f\n',z); else fprintf('\n El coeficiente de correlacion r='); fprintf('%5.5f\n',-z); end hold on ezplot(fun),grid on, plot(xi,yi,'m.'), hold off

Corremos el programa con un ejemplo de valores >> correl ingrese los valores de xi: [4 5 2 5 6 7 1 8 3 7] ingreso los valores de yi: [5 6 4 5 7 10 3 11 4 9] ingreso un valor de x para analisar la pendiente: 2 Los coeficientes de la ecuacion son: a b 1.130252100.97478992 La funcion lineal es y=1.13025210*x+0.97478992 El coeficiente de determinacion r^2=0.88900 El coeficiente de correlacion r=0.94287 >> Gráfica