s1Eva_IIIT2007_T2_AN Función Cobb-Douglas

usando los algoritmos para encontrar los polinomios de lagrange y búsqueda de raíz por la bisección, se obtiene:

 **** literal a:
Polinomio de Lagrange
[0, 6.16333333333368e-5*x**3 - 0.00527800000000012*x**2 + 0.254746666666668*x + 11.0, 0.00010365*x**3 - 0.008876*x**2 + 0.428425*x + 18.4997, 0.000140483333333337*x**3 - 0.0120304999999998*x**2 + 0.580686666666667*x + 25.0746]
Puntos f[25]:
[0, 15.0329375000000, 25.2823562500000, 34.2677562500002]
Polinoio con f[25]:
0.000586583333333373*x**3 - 0.0415150937500013*x**2 + 1.85978635416668*x
Estimado de [la,ka]:
29.7130898437501
 **** Literal b
La raiz se encuentra en:  25.3295898438

Revisar el algoritmo por partes: literal a y literal b.

# Interpolacion de Lagrange
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

def lagrange(xi,fi):
    n = len(xi)
    x = sp.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if j!=i:
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*fi[i]
    # Expande el polinomio
    px = polinomio.expand()
    return(px)

# INGRESO , Datos de prueba
M = np.array([[0, 0, 0, 0],
               [11.0000, 13.0813, 14.4768, 15.5563],
               [18.4997, 22.0000, 24.3470, 26.1626],
               [25.0746, 29.8189, 33.0000, 35.4608]])
li = np.array([0, 10, 20, 30])
kj = np.array([10, 20, 30, 40])

la = 25
ka = 25

# PROCEDIMIENTO
tamano = np.shape(M)
n = tamano[0]
m = tamano[1]
x = sp.Symbol('x')
poli = []
for i in range(0,n,1):
    xi = li
    fi = M[i,:]
    p = lagrange(xi,fi)
    poli.append(p)
# literal a, evalua en la
f_la =[]
for i in range(0,n,1):
    puntola = poli[i].subs(x,la)
    f_la.append(puntola)
poli_la = lagrange(li,f_la)
# literal a, evalua en ka
puntolaka = poli_la.subs(x,ka)


# Salida
print(' **** literal a:')
print('Polinomio de Lagrange')
print(poli)
print('Puntos f['+str(la)+']:')
print(f_la)
print('Polinoio con f['+str(ka)+']:')
print(poli_la)
print('Estimado de [la,ka]:')
print(puntolaka)

# literal b, usa f_ka con resultado 30

def biseccion(funcionx,a,b,tolera):
    """
    rango de búsqueda [a,b]
    busca hasta que (a-b)=tolera and cambia<0):
        c = (a+b)/2
        fc = funcionx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if (cambia<0):
            b=c
        else:
            a=c
        fa = funcionx(a)
        fb = funcionx(b)
        cambia = np.sign(fa)*np.sign(fb)
        tramo = abs(b-a)
    respuesta = c

    return(respuesta)
# Ingreso
a = kj[0]
b = kj[m-1]
tolera = 0.01

# PROCEDIMIENTO
# Busca raiz
fx = sp.lambdify(x,poli_la - 30)
fxi = fx(kj)

raiz = biseccion(fx,a,b,tolera)

# Salida
print(' **** Literal b')
print('La raiz se encuentra en: ',raiz)
plt.plot(kj,fxi)
plt.axhline(0)
plt.show()