222 - Topic 11 - Lagrange Multipliers

270 days ago by Professor222

Solving a 3D problem maximization/minimization problem with one constraint using the method of Lagrange Multipliers.

John Travis

Mississippi College


%auto var('x,y') # To start things off, enter the function f that you want to use. global func @interact def _(func = [1,2,3,4,5,6]): global f if func==1: f = x^2-y^2 elif func==2: f = x^2/4+10*x+3*y^2/20+12*y elif func==3: f = x^3*y^2 elif func==3: f = x^2-y^2+x*y elif func==4: f = y^3-3*y*x^2-3*y^2-3*x^2+1 elif func==5: f = x+y # This one is interesting with constraint 5 since constraint no continuous # This one also works with the more complicated constraint 6. elif func==6: f = 3*x+2*y+4 #func=8 #f = (x+y)/(x^2+y^2+1) #func=9 #f = x^2+3*x*y+y^2 #func=10 #f = 8*x^2+9*y^2 show(f) 
       
func 

Click to the left again to hide and once more to show the dynamic interactive window

       
-3*x^2*y + y^3 - 3*x^2 - 3*y^2 + 1
-3*x^2*y + y^3 - 3*x^2 - 3*y^2 + 1
%auto # Next, we need a constraint equation of the form g(x,y)= 0. (Note, "C" is inside g.) # Also, the constraint curve needs to be entered as a vector valued function r(t) = <u(t),v(t)> var('t') func=6 if func==1: cons=1 g(x,y) = x^2+y^2-1 u(t)=cos(t) v(t)=sin(t) tstart=0 tend=2*pi elif func==2: cons=2 g(x,y)=x^2/9+y^2/4-900 u(t) = 90*cos(t) v(t) = 60*sin(t) tstart=0 tend=2*pi elif func==3: cons=3 g(x,y)=x+y-100 u(t)=t v(t)=100-t tstart=120 tend=180 elif func==4: cons=4 g(x,y)=x^2+y^2-2 u(t) = sqrt(2)*cos(t) v(t) = sqrt(2)*sin(t) tstart = 0 tend = 2*pi elif func==5: cons=5 g(x,y)=x*y-185 u(t) = t v(t) = 185/t tstart=3 tend=40 elif func==6: cons=6 # A figure 8 g(x,y)=4*x^2-4*x^4-y^2 u(t) = cos(t) v(t) = sin(2*t) tstart=0 tend=2*pi else: print "Need to evaluate function cell first to pick a function to use." print "Using constraint",cons print "Using g = ",g print "Using u = ",u print "Using v = ",v print "Using interval ",tstart,"->",tend 
       
Using constraint 6
Using g =  (x, y) |--> -4*x^4 + 4*x^2 - y^2
Using u =  t |--> cos(t)
Using v =  t |--> sin(2*t)
Using interval  0 -> 2*pi
Using constraint 6
Using g =  (x, y) |--> -4*x^4 + 4*x^2 - y^2
Using u =  t |--> cos(t)
Using v =  t |--> sin(2*t)
Using interval  0 -> 2*pi

(Be certain to evaluate each cell--in order--from the top.  The results from each active cell should print below the cell.)

%auto fx = diff(f,x) fy = diff(f,y) gradf = [fx,fy] gx = diff(g,x) gy = diff(g,y) gradg = [gx,gy] var('lamb') solns_boundary = solve([gx==lamb*fx, gy==lamb*fy, g==0],(x,y,lamb),solution_dict=True) print 'Extrema on the boundary may occur at the lagrange points' show((sol[x],sol[y]) for sol in solns_boundary) for sol in solns_boundary: print "(",sol[x],",",sol[y],")" 
       
Extrema on the boundary may occur at the lagrange points
( 0 , 0 )
( 1/8*sqrt(-3/2*sqrt(137) + 55/2) ,
-1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137)
+ 55)) )
( -1/8*sqrt(-3/2*sqrt(137) + 55/2) ,
1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) +
55)) )
( 1/8*sqrt(3/2*sqrt(137) + 55/2) , 1/128*sqrt(2)*sqrt(3*sqrt(137) +
55)*(sqrt(137) - 3) )
( -1/8*sqrt(3/2*sqrt(137) + 55/2) , -1/128*sqrt(2)*sqrt(3*sqrt(137) +
55)*(sqrt(137) - 3) )
Extrema on the boundary may occur at the lagrange points
( 0 , 0 )
( 1/8*sqrt(-3/2*sqrt(137) + 55/2) , -1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)) )
( -1/8*sqrt(-3/2*sqrt(137) + 55/2) , 1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)) )
( 1/8*sqrt(3/2*sqrt(137) + 55/2) , 1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) )
( -1/8*sqrt(3/2*sqrt(137) + 55/2) , -1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) )
lagrange=[] xs=[] ys=[] for soln in solns_boundary: # prune off all complex roots if soln[x] not in RR: print 'Not using complex critical value at (',soln[x],',',soln[y],')' else: lagrange.append((soln[x],soln[y])) xs.append(soln[x]) ys.append(soln[y]) print print "Using" print "x = ",xs print "y = ",ys print "Lag = ",lagrange for sol in lagrange: print sol[0],sol[1],f(x=sol[0],y=sol[1]) 
       
Using
x =  [0, 1/8*sqrt(-3/2*sqrt(137) + 55/2), -1/8*sqrt(-3/2*sqrt(137) +
55/2), 1/8*sqrt(3/2*sqrt(137) + 55/2), -1/8*sqrt(3/2*sqrt(137) + 55/2)]
y =  [0, -1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) +
3*sqrt(-3*sqrt(137) + 55)), 1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) +
55) + 3*sqrt(-3*sqrt(137) + 55)), 1/128*sqrt(2)*sqrt(3*sqrt(137) +
55)*(sqrt(137) - 3), -1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) -
3)]
Lag =  [(0, 0), (1/8*sqrt(-3/2*sqrt(137) + 55/2),
-1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137)
+ 55))), (-1/8*sqrt(-3/2*sqrt(137) + 55/2),
1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) +
55))), (1/8*sqrt(3/2*sqrt(137) + 55/2), 1/128*sqrt(2)*sqrt(3*sqrt(137) +
55)*(sqrt(137) - 3)), (-1/8*sqrt(3/2*sqrt(137) + 55/2),
-1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3))]
0 0 4
1/8*sqrt(-3/2*sqrt(137) + 55/2)
-1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137)
+ 55)) -1/64*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) +
3*sqrt(-3*sqrt(137) + 55)) + 3/8*sqrt(-3/2*sqrt(137) + 55/2) + 4
-1/8*sqrt(-3/2*sqrt(137) + 55/2)
1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) +
55)) 1/64*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) +
3*sqrt(-3*sqrt(137) + 55)) - 3/8*sqrt(-3/2*sqrt(137) + 55/2) + 4
1/8*sqrt(3/2*sqrt(137) + 55/2) 1/128*sqrt(2)*sqrt(3*sqrt(137) +
55)*(sqrt(137) - 3) 1/64*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3)
+ 3/8*sqrt(3/2*sqrt(137) + 55/2) + 4
-1/8*sqrt(3/2*sqrt(137) + 55/2) -1/128*sqrt(2)*sqrt(3*sqrt(137) +
55)*(sqrt(137) - 3) -1/64*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3)
- 3/8*sqrt(3/2*sqrt(137) + 55/2) + 4
Using
x =  [0, 1/8*sqrt(-3/2*sqrt(137) + 55/2), -1/8*sqrt(-3/2*sqrt(137) + 55/2), 1/8*sqrt(3/2*sqrt(137) + 55/2), -1/8*sqrt(3/2*sqrt(137) + 55/2)]
y =  [0, -1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)), 1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)), 1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3), -1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3)]
Lag =  [(0, 0), (1/8*sqrt(-3/2*sqrt(137) + 55/2), -1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55))), (-1/8*sqrt(-3/2*sqrt(137) + 55/2), 1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55))), (1/8*sqrt(3/2*sqrt(137) + 55/2), 1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3)), (-1/8*sqrt(3/2*sqrt(137) + 55/2), -1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3))]
0 0 4
1/8*sqrt(-3/2*sqrt(137) + 55/2) -1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)) -1/64*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)) + 3/8*sqrt(-3/2*sqrt(137) + 55/2) + 4
-1/8*sqrt(-3/2*sqrt(137) + 55/2) 1/128*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)) 1/64*sqrt(2)*(sqrt(137)*sqrt(-3*sqrt(137) + 55) + 3*sqrt(-3*sqrt(137) + 55)) - 3/8*sqrt(-3/2*sqrt(137) + 55/2) + 4
1/8*sqrt(3/2*sqrt(137) + 55/2) 1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) 1/64*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) + 3/8*sqrt(3/2*sqrt(137) + 55/2) + 4
-1/8*sqrt(3/2*sqrt(137) + 55/2) -1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) -1/64*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) - 3/8*sqrt(3/2*sqrt(137) + 55/2) + 4
# Need to make these dependent upon the ranges of soln_boundary delx = (max(xs)-min(xs))/2 dely = (max(ys)-min(ys))/2 print delx print dely if delx==0: delx = abs(min(xs)/2) if dely==0: dely = abs(min(ys)/2) xmin = min(xs)-delx xmax = max(xs)+delx ymin = min(ys)-dely ymax = max(ys)+dely pretty_print(html('Function = %s'%str(func)+' and Constraint = %s'%str(cons))) F = plot3d(f,(x,xmin,xmax),(y,ymin,ymax),color='yellow') Curve = parametric_plot3d( (u(t),v(t),f(x=u(t),y=v(t))),(t,tstart,tend),color='red') Points = sum(point3d((sol[0],sol[1],f(x=sol[0],y=sol[1])),size=20) for sol in lagrange) # Points += point3d((-20,-40,f(x=-20,y=-40)),size=50) @interact def _(In3D = checkbox(default=True),Spin=checkbox(default=True)): show(F+Curve+Points,spin=Spin) if In3D: print 'You will want your 3D glasses...' show(F+Curve+Points,stereo='redcyan',spin=Spin,title='A plot') 
       
1/8*sqrt(3/2*sqrt(137) + 55/2) 1/128*sqrt(2)*sqrt(3*sqrt(137) + 55)*(sqrt(137) - 3) Function = 6 and Constraint = 6
In3D 
Spin 

Click to the left again to hide and once more to show the dynamic interactive window

 
       
H = contour_plot(f,(x,xmin,xmax),(y,ymin,ymax),contours=50,cmap='Blues') H += parametric_plot( (u(t),v(t)),(t,tstart,tend),color='red') H.show() print 'Numerical Approximations to the exact Lagrange solutions:' for k in range(len(solns_boundary)): x0 = N(solns_boundary[k][x]) y0 = N(solns_boundary[k][y]) show((x0,y0,f(x=x0,y=y0))) 
       
Numerical Approximations to the exact Lagrange solutions:

                                
                            
Numerical Approximations to the exact Lagrange solutions: