> # Sep 26, 2002 sfb
> # Extrema. Review of one variable
> fun:=x^3-6*x^2+11*x-6;
> #find critical points: derivative, set = 0 or undefined, solve
> fun1:=diff(fun,x);
> eqn:=fun1=0;
> ans:=solve(eqn);
> #look at the first critical point
> ans[1];
> #second derivative test: sign of the second derivative
> fun2:=diff(fun1,x);
> sign2:=eval(fun2,x=ans[1]);
> #positive, smile, local min
> #second critical point
> sign2:=eval(fun2,x=ans[2]);
> #negative, frown, local max
> #The test can fail if sign2 is zero for example x^3
> # In calculus 2 we did series which shows why this works
> #taylor=f(a)+f'(a)(x-a)+f''(a)(x-a)^2/2!+f'''(a)(x-a)^3/3!+....
> mapleseries:=series(exp(x),x=0,4);#note 4 is off by one.
> taylorpoly:=convert(mapleseries,'polynom');
> series(fun,x=ans[1],3);
> tpoly1:=convert(series(fun,x=ans[1],3),'polynom');
> check:=coeff(tpoly1,x-ans[1]);
> simplify(check);
> tpoly2:=convert(series(fun,x=ans[2],3),'polynom');
> check:=coeff(tpoly2,x-ans[2]);
> simplify(check);
> plot([fun,tpoly1,tpoly2],x=0.5..3.5,y=-2..2,color=[red,blue,green],title="cubic and two quadratic approximations"); #note,spell,yrange,list
> ## functions of two variables
> #need taylor series for two variables
> #f(a,b)+f_x(a,b)(x-a)+f_y(a,b)(y-b)+[f_xx(a,b)(x-a)^2+2f_xy(a,b)(x-a)(y-b)+f_yy(a,b)(y-b)^2]/2!+...
> mtaylor(exp(x+y),[x=2,y=1],3);
> #f+fx(x-a)+fy(y-b)+(fxx(x-a)^2+2fxy(x-a)(y-b)+fyy(y-b)^2)/2
> #need critical points, tangent plane horizontal or vertical
> #both fx and fy zero or undefined
> fun:=x^2*y^2-x^2-y^2+1;
> funx:=diff(fun,x);funy:=diff(fun,y);
> eqn:={funx=0,funy=0};
> ans:=solve(eqn);
> ans[1];
> tpoly1:=mtaylor(fun,ans[1],3);
> plot3d({fun,tpoly1,1},x=-1..1,y=-1..1,style=wireframe,axes=boxed,title="fun and taylor at (0,0)");
> ans[2];
> tpoly2:=mtaylor(fun,ans[2],3);
> plot3d({fun,tpoly2,0},x=0..2,y=0..2,style=wireframe,axes=boxed,title="fun and taylor at (1,1)");
> a:=plot3d(fun,x=0..2,y=0..2,style=wireframe,color=red):
> b:=plot3d(tpoly2,x=0..2,y=0..2,style=wireframe,color=blue):
> c:=plot3d(0,x=0..2,y=0..2,style=wireframe,color=green):
> with(plots):
> display(a,b,c,axes=boxed,title="fun and taylor at (1,1)");
> #the function has a saddle at the critical point (1,1)
> #
> #there is a second derivative test for functions of 2 variables
> #but it is more complicated than the one variable case.
> #First find 2nd partials
> funxx:=diff(funx,x);diff(fun,x,x);diff(fun,x$2);
> funyy:=diff(funy,y);
> funxy:=diff(funx,y);diff(fun,x,y);diff(fun,y,x);
> discriminate:=funxx*funyy-(funxy)^2;
> #determinate of the array
> #[ f_xx  f_xy]
> #[ f_yx  f_yy]
> # 2nd derivative test for 2 variables: @critical point
> # case 1 discriminate > 0 f_xx > 0  => local min
> # case 2 discriminate > 0 f_xx < 0  => local max
> # case 3 discriminate < 0           => saddle
> # case 4 discriminate = 0 The Test FAILS, it could be anything or nothing
> #critical point 1
> eval(discriminate,ans[1]);
> eval(funxx,ans[1]);
> #so it is a local max
> eval([discriminate,funxx],ans[2]);
> #so it is a saddle
> for i from 1 to nops([ans]) do eval([discriminate,funxx,ans[i]],ans[i]);end do;
> #we see the other critical points are saddles (discriminate < 0)
> # The prototypes
> f:=x^2+y^2;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> f:=-x^2-y^2;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> f:=x^2-y^2;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> f:=x*y;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> f:=x^3-3*x*y^2;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> eval(bigD,{x=0,y=0});
> f:=A*x^2+B*x*y+C*y^2;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> f:=A*x^2+2*B*x*y+C*y^2;f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> # Analytic geometry using rotation to remove the x*y term
> # you can do same here. There are perpendicular directions u, v
> # st the quadratic expressed in these terms has no u*v term.
> f:=x^2-x*y+y^2;
> #warning the following is magic
> with(linalg):
> eigenvectors(matrix([[1,-1/2],[-1/2,1]]));
> #the vectors are perpendicular
> dotprod(vector([1,1]),vector([-1,1]));
> ans:=solve({u=x-y,v=x+y},{x,y});
> eval(f,ans);simplify(%);
> eval(f,{x=u,y=u});
> eval(f,{x=-v,y=v});
> f_xx:=diff(f,x,x);bigD:=f_xx*diff(f,y,y)-diff(f,x,y)^2;
> #CPEM feb14 worksheet from last spring's calculus 3