#Let's solve for even higher degree
a = 2
b = 3
def makePoly(p):
g(x) = sum(p[j]*x**(j+1) for j in range(len(p)))
return g(x)-g(2)+b
def objFunc(p):
f = makePoly(p)
x1 = 99
try:
x1 = find_root(f, a, 99)
except:
pass
return numerical_integral(sqrt(1+f.derivative()**2), 0, x1)[0]
def solveRand(n):
startP = [10*random()-5 for _ in range(n)]
sol = minimize_constrained(objFunc, [(None, None)]*(len(startP)-1) +[(None, 0)], startP)
return (sol, objFunc(sol))
n = 4
simuN=10
bestSol = None
bestVal = 999999999
for j in range(simuN):
sol, val = solveRand(n)
if sol==None or val<bestVal:
bestSol = sol
bestVal = val
print (bestSol)
print (objFunc(bestSol))
g = Graphics()
solPoly = makePoly(bestSol)
g += plot(solPoly, 0, find_root(solPoly, a, 10))
g += points([(a,b)], color='red')
g += line([(0, b), (a,b)], color='red')
g += line([(a, 0), (a,b)], color='red')
g.show()
