# Construct function out of extrema, roots...

Is there a way in sage to construct and plot a smooth function which has a given set of data:

• the minima $(miA_1,miA_2)$, $(miB_1,miB_2)$,...
• the maxima $(maA_1,maA_2)$, $(maB_1,maB_2)$,...
• the roots $(nA1,nA2)$, $(nA1,nA2)$...
• a given behaviour for $x \to \pm \infty$.
edit retag close merge delete

I don't think so (this would be a strange feature to have), but you can do it analytically through calculus and other methods. Do you need help doing that?

( 2012-06-08 19:26:21 +0200 )edit

If you have a good idea how to do it, it would be nice if you could explain it.

( 2012-06-09 02:57:13 +0200 )edit

Sort by » oldest newest most voted

OK, the third and last time ;-)

L_roots = [-1,4]
L_max = [(-3,1)]
L_min = [(1,1)]
L_other = []

distinguish_min_max = False  # set to True or False

N = len(L_roots)+2*len(L_min)+2*len(L_max)+len(L_other)-1

param_list = [var("a%d" % k) for k in [0..N]]

f(x) = sum([ param_list[k]*x^k for k in [0..N]])
df = diff(f,x)
ddf = diff(df,x)
show(f)

eqns = [ f(x0)==0 for x0 in L_roots ]
eqns += [ f(P[0])==P[1] for P in L_min ]
eqns += [ df(P[0])==0 for P in L_min]
eqns += [ f(P[0])==P[1] for P in L_max ]
eqns += [ df(P[0])==0 for P in L_max ]
eqns += [ f(P[0])==P[1] for P in L_other ]

if distinguish_min_max:
eqns += [ ddf(P[0])> 0 for P in L_min]
eqns += [ ddf(P[0])<0 for P in L_max]

#for eq in eqns:
#    print eq

sol = solve(eqns,param_list,solution_dict=True)
show(sol)

try:
g = f.substitute(sol[0])
show(g)
G = Graphics()
G += g.plot(xmin=-10,xmax=10,ymin=-10,ymax=10)
G += points(L_min,color='red')
G += points(L_max,color='red')
G += points(L_other,color='black')
G += points([(x0,0) for x0 in L_roots],color='green')
G.show()
except:
print 'no solution'

more

Great, thanks very much!

( 2012-06-19 10:13:31 +0200 )edit

Here an example where you don't have to rewrite the code. I do not distinguish between minima and maxima. Controlling behavior for x to infinity: set some additional points.

L_roots = [-1,0]
L_min_max = [(1,1),(5,1)]
L_other = [(10,10)]
N = len(L_roots)+2*len(L_min_max)+len(L_other) -1

param_str = reduce(lambda x,y:x+y,['a%s,'%(k) for k in [0..N]])
param_list = var(param_str[:-1])

f(x) = sum([ param_list[k]*x^k for k in [0..N]])
df = diff(f,x)
show(f)

eqns = [ f(x0)==0 for x0 in L_roots ]
eqns += [ f(P[0])==P[1] for P in L_min_max ]
eqns += [ df(P[0])==0 for P in L_min_max ]
eqns += [ f(P[0])==P[1] for P in L_other ]
sol=solve(eqns,param_list,solution_dict=True)
show(f.substitute(sol[0]))

more

I think param_list = [var("a%d" % k) for k in [0..N]] is simpler.

( 2012-06-10 13:39:43 +0200 )edit

Thanks, that's great, almost what I want. Is there a way to implement into your code that it can distinguish between minima and maxima?

( 2012-06-10 15:50:34 +0200 )edit

DSM, you are right, thanks, we don't need param_str

( 2012-06-10 16:44:42 +0200 )edit

What kind of function do you want? Here a basic example for beginners. (http://www.hr.shuttle.de:9000/home/pu...)

# Minima / Maxima
A = (-1 , 1)
B = ( 3 , 3)

# Root
C = (1 , 0)

var("a0 a1 a2 a3 a4 a5 a6")

# f(x) = a5*x^5+ a4*x^4+a3*x^3 + a2*x^2+ a1*x + a0
f(x) = a4*x^4+a3*x^3 + a2*x^2+ a1*x + a0
# f(x) = a3*x^3 + a2*x^2+ a1*x + a0

df = diff(f,x);
ddf = diff(df,x);

xa,ya = A; xb,yb = B; xc,yc = C
# equations
eq1 = ( f(xa) == ya )
eq2 = ( df(xa) == 0 )
eq3 = ( f(xb) == yb )
eq4 = ( df(xb) == 0 )
eq5 = ( f(xc) == yc )
sol = solve([eq1,eq2,eq3,eq4,eq5],a0,a1,a2,a3,a4, solution_dict=True)
show(sol)
g = f.substitute(sol[0])
show(g)
plot(g,-4,4,ymin=-5,ymax=5)

more

Thanks, that's nice. Polynomial functions are fine. The problem is that if I want to add two roots or a minima to your example, I have to rewrite essentially the whole thing. How do you control the behaviour for $x\to \pm \infinity$? How do you distinguish between minima and maxima in your example?

( 2012-06-09 04:12:29 +0200 )edit

Minima will be at points where the second derivative is positive and maxima where it is negative. See http://en.wikipedia.org/wiki/Second_derivative_test.

( 2012-06-09 04:48:59 +0200 )edit

Yes that's clear. My question is about how to use this in the sage program

( 2012-06-09 07:26:21 +0200 )edit