Trazador C\u00fabico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n\n\n\n librer\u00eda<\/a><\/p>\n<\/div>\n\n\n\n Referencia<\/strong><\/em>: Burden 3.5 p105, Chapra 18.6.3 p532, Rodr\u00edguez 6.11.1 p244<\/p>\n\n\n\n Tiene como objetivo incorporar condiciones adicionales para la funci\u00f3n en los extremos del intervalo donde se encuentran los puntos o \"nodos\".<\/p>\n\n\n\n Por ejemplo, si los datos son los de una trayectoria en un experimento de f\u00edsica, se podr\u00eda disponer de la aceleraci\u00f3n el punto inicial de las muestras (izquierda) y de salida (derecha) del intervalo de las muestras.<\/p>\n\n\n\n El m\u00e9todo busca obtener un polinomio de tercer grado para cada sub-intervalo o tramo entre \"nodos\" consecutivos (j, j+1), de la forma:<\/p>\n\n\n\n S_j(x_j) = a_j + b_j (x-x_j) + c_j(x-xj)^2 + d_j(x-x_j)^3 <\/span>\n\n\n\n para cada j = 0, 1, ..., n-1<\/p>\n\n\n\n Para n<\/strong> datos existen n-1<\/strong> tramos, cuatro inc\u00f3gnitas (coeficientes) que evaluar por cada tramo, como resultado 4(n-1)<\/strong> inc\u00f3gnitas para todo el intervalo .<\/p>\n\n\n\n Para los t\u00e9rminos (x<\/strong>j+1<\/sub>- x<\/strong>j<\/sub>) de los tramos que se usan varias veces en el desarrollo, se simplifican como hj<\/sub>:<\/p>\n\n\n\n h_j = x_{j+1} - x_{j}<\/span>\n\n\n\n S_j(x_j) = a_j + b_j h_j + c_j h_j^2 + d_jh_j^3 <\/span>\n\n\n\n Para generar el sistema de ecuaciones, se siguen los siguientes criterios:<\/p>\n\n\n\n 1. Los valores de la funci\u00f3n deben ser iguales en los nodos interiores<\/p>\n\n\n\n S_j(x_{j+1}) = f(x_{j+1}) = S_{j+1}(x_{j+1}) <\/span>\n\n\n\n 2. Las primeras derivadas en los nodos interiores deben ser iguales<\/p>\n\n\n\n S'_j(x_{j+1}) = S'_{j+1}(x_{j+1}) <\/span>\n\n\n\n 3. Las segundas derivadas en los nodos interiores deben ser iguales<\/p>\n\n\n\n S''_j(x_{j+1}) = S''_{j+1}(x_{j+1}) <\/span>\n\n\n\n 4. El primer y \u00faltimo polinomio deben pasar a trav\u00e9s de los puntos extremos<\/p>\n\n\n\n f(x_0) = S_0(x_0) = a_0 <\/span>\n\n\n\n f(x_n) = S_n(x_n) = a_n <\/span>\n\n\n\n 5. Una de las condiciones de frontera se satisface:<\/p>\n\n\n\n 5a. frontera libre o natural<\/em><\/strong>: Las segundas derivadas en los nodos extremos son cero<\/p>\n\n\n\n S''(x_0) = S''(x_n) = 0 <\/span>\n\n\n\n 5b. frontera sujeta<\/em><\/strong>: las primeras derivadas en los nodos extremos son conocidas<\/p>\n\n\n\n S'(x_0) = f'(x_0)<\/span>\n\n\n\n S'(x_n) = f'(x_n) <\/span>\n\n\n\n Tarea<\/strong>: Revisar en los libros, textos el planteamiento de las ecuaciones para resolver el sistema que se genera al plantear el polinomio.<\/p>\n\n\n\n Ubicar en el texto las ecuaciones resultantes, que son las que se aplicar\u00e1n en el algoritmo.<\/p>\n\n\n\n Trazador C\u00fabico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n\n\n\n librer\u00eda<\/a><\/a><\/p>\n<\/div>\n\n\n\n El algoritmo parte de lo desarrollado para \"trazadores lineales\", donde se presentaron los bloques para:<\/p>\n\n\n\n Del algoritmo b\u00e1sico se modifica entonces el bloque del c\u00e1lculo de los polinomios de acuerdo al planteamiento de formulas y procedimientos para trazadores c\u00fabicos naturales (enviado a revisar como tarea).<\/p>\n\n\n con lo que se obtiene el resultado por cada tramo<\/p>\n\n\n\n Trazador C\u00fabico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n\n\n\n librer\u00eda<\/a><\/a><\/p>\n<\/div>\n\n\n\n Para el trazado de la gr\u00e1fica se a\u00f1ade al final del algoritmo:<\/p>\n\n\n Trazador C\u00fabico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n\n\n\n librer\u00eda<\/a><\/a><\/p>\n<\/div>\n\n\n\n Referencia<\/strong>: https:\/\/docs.scipy.org\/doc\/scipy\/reference\/generated\/scipy.interpolate.CubicSpline.html#scipy.interpolate.CubicSpline<\/a><\/p>\n\n\n\n La librer\u00eda Scipy dispone de una funci\u00f3n para calcular los coeficientes de los trazadores c\u00fabicos. Para el ejemplo del ejercicio se realiza el c\u00e1lculo y los polinomios.<\/p>\n\n\n\n Las instrucciones en Python con la librer\u00eda Scipy para trazadores c\u00fabicos natural son:<\/p>\n\n\n La gr\u00e1fica se puede generar con el bloque anterior<\/p>\n\n\n\n Trazador C\u00fabico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n\n\n\n librer\u00eda<\/a><\/a><\/p>\n<\/div>\n\n\n\n Si los polinomios entre cada tramo no se igualan sus velocidades y aceleraciones<\/strong>, se presentan puntos de inflexi\u00f3n en cada nodo.<\/p>\n\n\n\n El resultado para en la trayectoria de un veh\u00edculo tripulado podr\u00eda generar una experiencia semejante a la mostrada en los siguientes videos.<\/p>\n\n\n\n Referencia<\/strong>: [1] Car sales woman scares customers. maxman.tv. 4 enero 2016<\/p>\n\n\n\n
\n\n\n\n1. Trazador C\u00fabico - Planteamiento<\/h2>\n\n\n\n
![\"spline](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/07\/spline03_curvigrafo.jpg\")
<\/figure>\n\n\n\n![\"spline](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/07\/spline04naturallibre.png\")
<\/figure>\n\n\n\n
\n\n\n\n
\n\n\n\n2. Algoritmo en Python<\/h2>\n\n\n\n
\n
\n# Trazador cubico natural\nimport numpy as np\nimport sympy as sym\n\n# INGRESO\nxi = [0.1 , 0.2, 0.3, 0.4]\nfi = [1.45, 1.8, 1.7, 2.0]\n\ntitulo = 'Trazador c\u00fabico natural o libre'\n# natural, 2da derivadas en los nodos extremos son cero\n\n# Algoritmos como funcion\n\ndef traza3natural(xi,yi, vertabla=False):\n ''' trazador cubico natural,\n 2da derivadas en los nodos extremos son cero\n '''\n # Vectores como arreglo, numeros reales\n xi = np.array(xi,dtype=float)\n yi = np.array(yi,dtype=float)\n n = len(xi)\n\n h = np.diff(xi,n=1) # h tamano de paso\n\n # Sistema de ecuaciones\n A = np.zeros(shape=(n-2,n-2), dtype=float)\n B = np.zeros(n-2, dtype=float)\n S = np.zeros(n, dtype=float)\n A[0,0] = 2*(h[0]+h[1])\n A[0,1] = h[1]\n B[0] = 6*((yi[2]-yi[1])\/h[1] - (yi[1]-yi[0])\/h[0])\n for i in range(1,n-3,1):\n A[i,i-1] = h[i]\n A[i,i] = 2*(h[i]+h[i+1])\n A[i,i+1] = h[i+1]\n B[i] = 6*((yi[i+2]-yi[i+1])\/h[i+1] - (yi[i+1]-yi[i])\/h[i])\n A[n-3,n-4] = h[n-3]\n A[n-3,n-3] = 2*(h[n-3]+h[n-2])\n B[n-3] = 6*((yi[n-1]-yi[n-2])\/h[n-2] - (yi[n-2]-yi[n-3])\/h[n-3])\n\n # Resolver sistema de ecuaciones, S\n r = np.linalg.solve(A,B)\n #S[0] = 0; S[n-1] = 0 # Definidos en cero\n S[1:n-1] = r[0:n-2]\n\n # Coeficientes\n a = np.zeros(n-1, dtype = float)\n b = np.zeros(n-1, dtype = float)\n c = np.zeros(n-1, dtype = float)\n d = np.zeros(n-1, dtype = float)\n for j in range(0,n-1,1):\n a[j] = (S[j+1]-S[j])\/(6*h[j])\n b[j] = S[j]\/2\n c[j] = (yi[j+1]-yi[j])\/h[j] - (2*h[j]*S[j]+h[j]*S[j+1])\/6\n d[j] = yi[j]\n\n # Polinomio trazador\n x = sym.Symbol('x')\n px_tabla = []\n for j in range(0,n-1,1):\n pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2 + c[j]*(x-xi[j])+ d[j]\n px_tabla.append(pxtramo)\n \n if vertabla==True:\n print('h:',h)\n print('A:') ; print(A)\n print('B:',B) ; print('S:',S)\n print('r:',r)\n print('coeficientes[ ,tramo]:')\n print('a:',a) ; print('b:',b)\n print('c:',c) ; print('d:',d)\n\n return(px_tabla)\n\n# PROGRAMA ---------------\n# PROCEDIMIENTO\nn = len(xi)\n# Obtiene los polinomios por tramos\npx_tabla = traza3natural(xi,fi,vertabla=True)\n\n# SALIDA\nprint(titulo)\nprint('Polinomios por tramos: ')\nfor i in range(0,n-1,1):\n print('x = [',xi[i],',',xi[i+1],']')\n print('expresion:',px_tabla[i])\n print('simplificado:')\n sym.pprint(px_tabla[i].expand())\n<\/pre><\/div>\n\n\nh: [0.1 0.1 0.1]\nA:\n[[0.4 0.1]\n [0.1 0.4]]\nB: [-27. 24.]\nr: [-88. 82.]\nS: [ 0. -88. 82. 0.]\ncoeficientes[ ,tramo]:\na: [-146.66666667 283.33333333 -136.66666667]\nb: [ 0. -44. 41.]\nc: [4.96666667 0.56666667 0.26666667]\nd: [1.45 1.8 1.7 ]\nTrazador c\u00fabico natural o libre\nPolinomios por tramos: \nx = [ 0.1 , 0.2 ]\nexpresion: 4.96666666666667*x - 146.666666666667*(x - 0.1)**3 + 0.953333333333333\nsimplificado:\n 3 2 \n- 146.666666666667\u22c5x + 44.0\u22c5x + 0.566666666666666\u22c5x + 1.1\nx = [ 0.2 , 0.3 ]\nexpresion: 0.566666666666666*x + 283.333333333333*(x - 0.2)**3 - 44.0*(x - 0.2)**2 + 1.68666666666667\nsimplificado:\n 3 2 \n283.333333333333\u22c5x - 214.0\u22c5x + 52.1666666666667\u22c5x - 2.34\nx = [ 0.3 , 0.4 ]\nexpresion: 0.266666666666665*x - 136.666666666667*(x - 0.3)**3 + 41.0*(x - 0.3)**2 + 1.62\nsimplificado:\n 3 2 \n- 136.666666666667\u22c5x + 164.0\u22c5x - 61.2333333333333\u22c5x + 9.0\n>>><\/code><\/pre>\n\n\n\n
\n\n\n\n
\n\n\n\n3. Gr\u00e1fica con Python<\/h2>\n\n\n\n
\n# GRAFICA --------------\nimport matplotlib.pyplot as plt\n\nmuestras = 10 # entre cada par de puntos\n\n# Puntos para grafica en cada tramo\nxtraza = np.array([],dtype=float)\nytraza = np.array([],dtype=float)\ni = 0\nwhile i<n-1: # cada tramo\n xtramo = np.linspace(xi[i],xi[i+1],muestras)\n # evalua polinomio del tramo\n pxk = sym.lambdify('x',px_tabla[i])\n ytramo = pxk(xtramo)\n # vectores de trazador en x,y\n xtraza = np.concatenate((xtraza,xtramo))\n ytraza = np.concatenate((ytraza,ytramo))\n i = i + 1\n\n# Graficar\nplt.plot(xtraza,ytraza, label='trazador')\nplt.plot(xi,fi,'o', label='puntos')\nplt.title(titulo)\nplt.xlabel('xi')\nplt.ylabel('px(xi)')\nplt.legend()\nplt.show()\n<\/pre><\/div>\n\n\n
\n\n\n\n
\n\n\n\n4. Trazador c\u00fabico natural en librer\u00eda Scipy<\/h2>\n\n\n\n
Trazador c\u00fabico natural o libre\ncoeficientes:\na: [-146.66666667 283.33333333 -136.66666667]\nb: [ 8.8817842e-15 -4.4000000e+01 4.1000000e+01]\nc: [4.96666667 0.56666667 0.26666667]\nd: [1.45 1.8 1.7 ]\nPolinomios por tramos: \nx = [ 0.1 , 0.2 ]\nexpresi\u00f3n: 4.96666666666667*x - 146.666666666667*(x - 0.1)**3 + 8.88178419700125e-15*(x - 0.1)**2 + 0.953333333333333\nsimplifica:\n 3 2 \n- 146.666666666667\u22c5x + 44.0\u22c5x + 0.566666666666662\u22c5x + 1.1\nx = [ 0.2 , 0.3 ]\nexpresi\u00f3n: 0.566666666666667*x + 283.333333333333*(x - 0.2)**3 - 44.0*(x - 0.2)**2 + 1.68666666666667\nsimplifica:\n 3 2 \n283.333333333333\u22c5x - 214.0\u22c5x + 52.1666666666667\u22c5x - 2.34\nx = [ 0.3 , 0.4 ]\nexpresi\u00f3n: 0.266666666666665*x - 136.666666666667*(x - 0.3)**3 + 41.0*(x - 0.3)**2 + 1.62\nsimplifica:\n 3 2 \n- 136.666666666667\u22c5x + 164.0\u22c5x - 61.2333333333333\u22c5x + 9.0<\/code><\/pre>\n\n\n\n\n# Trazador c\u00fabico natural o libre con Scipy\nimport numpy as np\nimport scipy as sp\nimport sympy as sym\n\n# INGRESO\nxi = [0.1 , 0.2, 0.3, 0.4]\nfi = [1.45, 1.8, 1.7, 2.0]\n\ntitulo = 'Trazador c\u00fabico Natural o libre'\n# natural, 2da derivadas en los nodos extremos son cero\ncs_tipo = ((2, 0.0), (2, 0.0)) # natural\n\n# PROCEDIMIENTO\n# Vectores como arreglo, numeros reales\nxi = np.array(xi,dtype=float)\nfi = np.array(fi,dtype=float)\nn = len(xi)\n\n# coeficientes de trazadores c\u00fabicos\ntraza_cub = sp.interpolate.CubicSpline(xi,fi,bc_type=cs_tipo )\ntraza_coef = traza_cub.c # coeficientes\na = traza_coef[0]\nb = traza_coef[1]\nc = traza_coef[2]\nd = traza_coef[3]\n\n# Polinomio trazador\nx = sym.Symbol('x')\npx_tabla = []\nfor j in range(0,n-1,1):\n pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2 + c[j]*(x-xi[j])+ d[j]\n px_tabla.append(pxtramo)\n\n# SALIDA\nprint(titulo)\nprint('coeficientes:')\nprint('a:',a)\nprint('b:',b)\nprint('c:',c)\nprint('d:',d)\nprint('Polinomios por tramos: ')\nfor i in range(0,n-1,1):\n print('x = [',xi[i],',',xi[i+1],']')\n print('expresion:',px_tabla[i])\n print('simplifica:')\n sym.pprint(px_tabla[i].expand())\n<\/pre><\/div>\n\n\n
\n\n\n\n
\n\n\n\n![\"spline](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/07\/spline02_dimensiones.png\")
<\/figure>\n\n\n\n