# Práctica del tema 2
Blanca Cano Camarero


```python
# import necesarios para la práctica
import sympy as sp
import numpy as np

## ficheros propios
#### De vuelve el polinomio de interpolación de newton 
import polinomioNewton as pn # polinomioNewton([a-h,a,a+h])

##### Algunas fórmulas de integración 
import formulasIntegracion as fs

#### cuadratura gaussina
import nodosGaussiana as ng
```

# Ejercicio 1
1.- Obtenga mediante interpolación en el espacio $\mathbb{P}_2$ una fórmula para aproximar $f''(a)$ del tipo combinación
de $f(a-h)$, $f(a)$ y $f(a+h)$.



```python
a,h,z = sp.symbols('a,h,z')
x = [a-h, a, a+h]

#llamamos a nuestra funcioncilla que nos calcula polinomios de interpolación de newton
## (para ver su código mirar en polinomioNewton.py)
p = pn.polinomioNewton(x)

## Finalemente, derivamos pa obtener nuestra fórmula buscada
formula=sp.diff(p,z,2).subs({z:a}).simplify() # esta será la fórmula
# obtenida esta vez para aproximar la derivada segunda f''(a)
print(formula)
formula
```

    (-2*f(a) + f(a - h) + f(a + h))/h**2





$\displaystyle \frac{- 2 f{\left(a \right)} + f{\left(a - h \right)} + f{\left(a + h \right)}}{h^{2}}$



## Ejercicio 2
2.- Con la fórmula obtenida en el ejercicio 1, halle una tabla de aproximaciones y errores de $f_1''(2.5)$, siendo $f_1(x)=x^x$, para $h=10^{-i},\; i=1,\ldots,5.$



```python
## derivada aproximada 
def d2(f):
    return (-2*f(a) + f(a - h) + f(a + h))/h**2

def g(x):
    return x**x



### derivada exacta
g2 = sp.lambdify(z,sp.diff(z**z, z, 2), 'numpy')

#####
def g2_aprox(formula, g,mia, eh):
    """Devuelve la derivada aproximada de g con la fórmula "formula" en el punto mia, con h valiendo eh
    Argumentos:
        formula: aplicación lineal, con la que este caso aproximamos derivadas
        g la función de la que queremos calcular la derivada
        mia: punto en el que geremos calcular a
        eh: variación de h
    """
    return formula(g).evalf(subs={a:mia,h:10**(-eh)})

### hacemos una función pa representar la tablilla

def tabla(formula,faprox, fexac, nodos, vh):
    print(f' h          | a \t  | f(a) aprox \t | error'.expandtabs(25))
    print('-'*65)
    for x in nodos:
        valor_exacto = fexac(x)
        for mih in vh:
            v = g2_aprox(formula,faprox,x,mih)
            print(f'10^(-{mih})    | {x} \t | {v} \t|{abs(v-valor_exacto)}'.expandtabs(25))


tabla(d2,g,g2,[2.5],[*range(1,6)])    

```

    h          | a            | f(a) aprox            | error
    -----------------------------------------------------------------
    10^(-1)    | 2.5          | 40.4205682979584      |0.178903441570895
    10^(-2)    | 2.5          | 40.2434502309567      |0.00178537456920225
    10^(-3)    | 2.5          | 40.2416827097676      |0.0000178533801502567
    10^(-4)    | 2.5          | 40.2416650349212      |1.78533753114607E-7
    10^(-5)    | 2.5          | 40.2416648581728      |1.78533099415290E-9


## Ejercicio 3
3.- Sea $f_2(x)=\frac{x^2+40}{x+\sqrt{5x}+7}$. Calcule una tabla que recoja las derivadas de $f_2$ en $x_i=1,2,\ldots,10$, utilizando alguna de las fórmulas de derivación numérica de primer orden obtenidas al inicio de la práctica, con $h=10^{-3}$, y muestre al mismo tiempo el error cometido en cada punto. Repita el ejercicio con la fórmula centrada obtenida para la derivada primera y, finalmente, para la obtenida en el ejercicio 1 (con respecto a la segunda derivada).


```python
def f2(x):
    return (x**2 +40) / (x+(5*x)**(1/2)+7)


### Vamos a comenzar con la segunda derivada
ddf=sp.diff( (z**2 +40) / (z+(5*z)**(1/2)+7), z, 2)
ddf=sp.lambdify(z,ddf, 'numpy')
x = [*range(1,11)]

### tabla para la segunda derivada 
print('Las aproximaciones para la segunda derivada son ')
tabla(d2,f2,ddf,x,[3])

## ________Para la fórmula centrada ____
print('\nDerivamos una vez el polinomio de newton, para obtener la primera derivada')
formula=sp.diff(p,z,1).subs({z:a}).simplify() # esta será la fórmula
# obtenida esta vez para aproximar la derivada segunda f''(a)
print(formula)

#copiamso en el return la fórmula anterior
def df_aprox(f):
    return (-f(a - h) + f(a + h))/(2*h)

#calculamos primera derivada exacta
df=sp.diff( (z**2 +40) / (z+(5*z)**(1/2)+7), z, 1)
df=sp.lambdify(z,df, 'numpy')

print('\n\n Las proximaciones para la primera derivada')
tabla(df_aprox,f2,df,x,[3] )

```

    Las aproximaciones para la segunda derivada son 
     h          | a            | f(a) aprox            | error
    -----------------------------------------------------------------
    10^(-3)    | 1            | 0.676265261235160     |1.62949784665578E-7
    10^(-3)    | 2            | 0.283220380779308     |1.66032020909590E-8
    10^(-3)    | 3            | 0.168340324507543     |4.57942220255525E-9
    10^(-3)    | 4            | 0.114907314792673     |1.89762004720873E-9
    10^(-3)    | 5            | 0.0846224312804510    |9.93457316411650E-10
    10^(-3)    | 6            | 0.0654364319768901    |6.12947123612706E-10
    10^(-3)    | 7            | 0.0523721747995802    |4.30544377927333E-10
    10^(-3)    | 8            | 0.0430109452392692    |3.36394023747744E-10
    10^(-3)    | 9            | 0.0360420060101230    |2.86374354785401E-10
    10^(-3)    | 10           | 0.0306970069224294    |2.60408326846484E-10
    
    Derivamos una vez el polinomio de newton, para obtener la primera derivada
    (-f(a - h) + f(a + h))/(2*h)
    
    
     Las proximaciones para la primera derivada
     h          | a            | f(a) aprox            | error
    -----------------------------------------------------------------
    10^(-3)    | 1            | -0.633413983453862    |1.41948961474014E-7
    10^(-3)    | 2            | -0.203730021215820    |2.98523985808874E-8
    10^(-3)    | 3            | 0.0135536643822118    |1.22135472690327E-8
    10^(-3)    | 4            | 0.152356375979550     |6.46680242688547E-9
    10^(-3)    | 5            | 0.250865047976407     |3.92670773674553E-9
    10^(-3)    | 6            | 0.325234483749073     |2.59700044802358E-9
    10^(-3)    | 7            | 0.383753087445717     |1.82151538297148E-9
    10^(-3)    | 8            | 0.431201819322740     |1.33390964940361E-9
    10^(-3)    | 9            | 0.470566738047861     |1.00977470829378E-9
    10^(-3)    | 10           | 0.503824069630687     |7.84850628932077E-10



```python

```


```python

# Algunas de las definidas en formulasSimples.py
# declaramos las fórmulas
def rectangulo_izquierdo(f,a,b):
    return f(a)*(b-a)

def simpson(f, a, b):
    return 1/6*(f(a)+4*f((a+b)/2)+f(b))*(b-a)

def trapecio(f,a,b):
    return (f(a)+f(b))/2*(b-a)

def formula_compuesta (f,formula,a,b,n):
    h = (b-a)/n
    return sum([formula(f,a+i*h,a+(i+1)*h) for i in range(n)])



```

# Ejericicios 4 y 5

4.- Divida el intervalo $[1,2]$ en 100 partes iguales y aplique las fórmulas del rectángulo, Simpson y trapecio compuestas para aproximar la integral en dicho intervalo de $f_1$. Compare dichos resultados.

5.- Repita el ejercicio 4 para $f_2$. 


```python
def f1(x):
    return x**x
print("___Integral para f1___")
print(f'Integra por rectángulo izquierdo {formula_compuesta(f1,rectangulo_izquierdo,1,2,100)}')
print(f'Integra por trapecio {formula_compuesta(f1,trapecio,1,2,100)}')
print(f'Integra por simpson {formula_compuesta(f1,simpson,1,2,100)}')

def f2(x):
    return (x*x+40)/(x+sqrt(5*x)+7)
print("___Integral para f2___")
print(f'Integra por rectángulo izquierdo {formula_compuesta(f2,rectangulo_izquierdo,1,2,100)}')
print(f'Integra por trapecio {formula_compuesta(f2,trapecio,1,2,100)}')
print(f'Integra por simpson {formula_compuesta(f2,simpson,1,2,100)}')

```

    ___Integral para f1___
    Integra por rectángulo izquierdo 2.0354943390855578
    Integra por trapecio 2.0504943390855574
    Integra por simpson 2.0504462346235286
    ___Integral para f2___
    Integra por rectángulo izquierdo 3.77852320278209
    Integra por trapecio 3.77658469845732
    Integra por simpson 3.77658111777024


Comentario: 
Como podemos observar dan un resultado similar, siend trapecio y simpson muy parecidas entre sí.  

## Ejercicio 6
6.- Sea $f_3(x)=x^{15} e^x$ en $[0,2]$. Vamos a dividir el intervalo en $10\times 2^n$ subintervalos, es decir, $10,\,20,\,40,\, 80,\ldots $ y a aplicar la fórmula de Simpson compuesta hasta que la diferencia entre dos aproximaciones consecutivas (por ejemplo, podrían
ser con $20$ y $40$ subintervalos) sea menor que $10^{-2}$, dando en tal caso por buena la última aproximación obtenida. Programe
y calcule dicha aproximación. Compare ambas aproximaciones con
el valor exacto.


```python
def f3(x):
    return x**(15) * np.exp(x)

## cálculo de valor exacto por regla de barrow
x = sp.symbols("x")
f3s = x**15 * sp.exp(x) 
If3s = sp.integrate(f3s,x)
exacto = If3s.subs(x,2)-If3s.subs(x,0)
print(f'Aplicando la regla de barrow el valor es {exacto.evalf()}')
diferencia = 9
y = fs.formula_compuesta(f3,fs.simpson,0,2,10)
n= 1
while diferencia > 10**(-2):
    ys = fs.formula_compuesta(f3,fs.simpson,0,2,10*2**(n))
    n += 1
    diferencia = abs(ys-y)
    y = ys
print (f'La a integral aproximada (en la iteración {n})  es {y}')
print (f'El error es de {abs(exacto-y).evalf()}')
```

    Aplicando la regla de barrow el valor es 27062.7024138996
    La a integral aproximada (en la iteración 6)  es 27062.702480891196
    El error es de 0.0000213214760707824


## Ejercicio 7

7.- Calcule las fórmulas gaussianas con $2$ y $3$ nodos,en el intervalo $[-1,1]$, siendo la función peso el valor absoluto de la variable. Aplíquelas para aproximar la función $x\; e^x$ en $[-1,1]$ y compare los resultados con el valor exacto (organizando los cálculos de forma adecuada).



```python
#calcula nodos
def nodos_gaussiana (a, b, n,w):
    """ Devuelve los nodos de la fórmula de la gaussiana
    a,b: intervalo [a,b] de definición
    n: número de nodos 
    w: función peso
    """
    x = sp.Symbol("x")
    c = list(sp.symbols('c0:'+ str(n)))
    pi = np.prod([ (x - c[i]) for i in range(n)])
    #print('Voy a computar la integral (esto puede tardar)...')
    I = [sp.integrate(pi*w(x)*x**i,(x, a, b)) for i in range(n)]
    #print(f'El valor de la integrales es {I}')
    #print('la solución buscada es: ')
    s = sp.solve(I,c)
    return list(s[0])

```


```python
a = -1
b= 1
n = 3
def w(x):
    return abs(x)#sp.sqrt(1/(1-x**2))
z=sp.Symbol('z')
for i in range(2,n+1):
    q =pn.polinomioNewton(nodos_gaussiana(a,b,i,w))
    print(f'fóruma para {i} nodos')
    print(sp.integrate(q,(z,a,b)))

f_exacta = x*sp.exp(x)
print(f'\nEl valor exacto es {sp.N(sp.integrate(f_exacta,(x,a,b)))}')

## vemos ahora a por el valro aproximado
def f_aprox(x):
    return x*np.exp(x)

def dosNodos(f):
    return f(-np.sqrt(2)/2) + f(np.sqrt(2)/2)
def tresNodos(f):
    return f(0) + f(-np.sqrt(6)/3)/2 + f(np.sqrt(6)/3)/2

print(f'Para la aproximación de dos nodos el valor es {dosNodos(f_aprox)}')
print(f'Para la aproximación de tres nodos el valor es {tresNodos(f_aprox)}')
```

    fóruma para 2 nodos
    f(-sqrt(2)/2) + f(sqrt(2)/2)
    fóruma para 3 nodos
    f(0) + f(-sqrt(6)/3)/2 + f(sqrt(6)/3)/2
    
    El valor exacto es 0.735758882342885
    Para dos nodos el valor es 1.0854416412726071
    Para tres nodos el valor es 0.7432494342785245


## Ejercicio 8

8.- Programar las técnicas de integración de Romberg y adaptativa, para después aplicarlas a la aproximación de la siguiente integral $$\int_a^b p(x)\, dx$$
siendo  $\;a=\displaystyle\min_{0\leq i\leq 7}{d_i}$, $\;b=\displaystyle\max_{0\leq i\leq 7}{d_i}$ y 
$$p(x)=d_0 + d_1 x + d_2 x^2 + d_3 x^3+ d_4 x^4 + d_5 x^5 + d_6 x^6 + d_7 x^7 $$
(siendo $d_0, d_1, \ldots, d_7$ los dígitos).


```python
## El polinomio pedido es: 
coef = [2,5,3,7,8,2,9,1] 
a = min(coef)
b = max(coef)

def p(x):
    return sum([coef[i]*x*i for i in range(len(coef))])

# fórmula en fichero formulasIntegracion.py
### ADAPTATIVA #######

print(f'El resultado de la adaptativa es {fs.adaptativa(fs.simpson , p, a , b)}')

### ROMBER 
def T2n(tn,f,h,a,n):
    """Fórmula recursiva trapecio
    tn: valor de T_n
    f: función a evaluar
    h: diferencia partición h_{n+1}=h_{n}/2, por eso los nodos por calcular son x_i+h_{n+1}     = (a + m*h_n)+h_{n+1} = (a + m*2h_{n+1})+h_{n+1}
    a: intervalo inferior partición
    n: proviene de h{n+1}=(b-a)/2^{n+1}
    """
    return 1/2*(tn)+h*sum(map((lambda e: f(a+(e*2+1)*h )),[*range(n)]))

def Romberg(f,a,b,tol=(10**(-5)),length=20): 
    """ Integrar calculada por el método de Romberg, 
    devuelve el número ancho de columan exploradas 
    f: función sobre la que integrar
    a: intervalo inferior sobre el que integrar
    b: intervalo superior integración
    length: Ancho de la tabla maximo
    tol: diferencia que debe haber entre R(n-1,n-1) y R(n,n)
    """
    #matriz que contentendrá datos integración Romberg
    l = [ [None]*(i+1) for i in range(length)] # ojo que viene dada por l[fila][columan]

    #Calulamos R(N,0), es decir la primera columa de la tabla
    l[0][0]=fs.formula_compuesta(f,fs.trapecio, a, b,1)
    h = (b-a)
    n=1
    for i in range(1,length):
        h=h/2
        n = n*2
        l[i][0] = T2n(l[i-1][0],f,h,a,n)
             
        

    ## camos a proceder ahora con los elementos de la forma R(j,k)
    columna = 1
    diferencia = 99999 
    #for columna in range(1,length):
    while( columna < length and diferencia > tol):
        for fila in range(columna, length):
            l[fila][columna] = (4**columna*l[fila][columna -1] - l[fila - 1][columna -1] ) / (4**columna -1) 
        
        diferencia = abs(l[columna][columna]-l[columna-1][columna-1])
        columna += 1
    columna -= 1
    return (l[columna][columna],columna)

r = Romberg(p,a,b)
print(f'El resultado de mi fórmula de romber es {r[0]} tras {r[1]} iteraciones')
```

    El resultado de la adaptativa es 5400.0
    El resultado de mi fórmula de romber es 19439.983687473403 tras 19 iteraciones



```python

```

c79665ebdf0d16d3595e10571bdbe9e7740eafad
ipynb
Jupyter Notebook
practica2-integracia_numeria/Entregable.ipynb
BlancaCC/metodosNumericosII
73fd6d8bc202a34af0d95d9ac17a3b671895314d
MIT
practica2-integracia_numeria/Entregable.ipynb
BlancaCC/metodosNumericosII
73fd6d8bc202a34af0d95d9ac17a3b671895314d
MIT
practica2-integracia_numeria/Entregable.ipynb
BlancaCC/metodosNumericosII
73fd6d8bc202a34af0d95d9ac17a3b671895314d
MIT
Qwen/Qwen-72B
1. YES
2. YES