> #sfb feb 21,2002
> # Part 0 -- Debugging
> # Is it me or is it Maple?
> # Try something simple
> # 2+2;
> # plot(x^2,x=0..1);
> # Look at the intermediate steps.
> # Change a:=foo: to a:=foo; or a:=plot(foo): to a:=plot(foo):a;
> # 95% of the errors are typos or something simple like something isn't defined.
> # Send email, or come by.
> #
> # Part 1 -- Pretty Integrals
> # Maple uses `int(f,x)' for indefinite integrals and `int(f,x=a..b)' for definite ones
> int(x^2,x);   int(x^2,x=0..1);
> # With a little cut 'n paste one can always make an integral pretty
> Int(x^2,x)=int(x^2,x);   Int(x^2,x=0..1)=int(x^2,x=0..1);
> 'int(x^2,x)'=int(x^2,x);   'int(x^2,x=0..1)'=int(x^2,x=0..1);
> # Works for double and triple integrals too
> Int(Int(x^2+y^2,y=0..x^2),x=-1..1)=int(int(x^2+y^2,y=0..x^2),x=-1..1);
> Int(Int(Int(x*y*z^2,z=0..x+y),y=0..x),x=0..1)=int(int(int(x*y*z^2,z=0..x+y),y=0..x),x=0..1);
> #What do you need to remember about integration? Well everything, but for this class
> #Standard functions, Substitution and the Half-angle formula's do 99% of the problems.
> (sin(x))^2=(1-cos(2*x))/2;     (cos(x))^2=(1+cos(2*x))/2;
> #
> # Part 2 -- Least Squares Fit
> # Data [[x1,y1],[x2,y2], ... [xn,yn]]
> # Problem 15, Section 14.2, Page 735
> postage:=[[1920,2],[1932,3],[1958,4],[1963,5],[1968,6],[1971,8],[1974,10],[1975,13],[1978,15],[1981,20],[1985,22],[1988,25],[1991,29],[1995,32]];
> plot(postage,title="Cost of a first class stamp");
> plot(postage,title="Cost of a first class stamp",style=point);
> dataplot:=plot(postage,title="Cost of a first class stamp",style=point,symbol=box,symbolsize=24):dataplot;
> ?leastsquare
> ?fit
> postx:=[seq(postage[i][1],i=1..nops(postage))];
> posty:=[seq(postage[i][2],i=1..nops(postage))];
> fit[leastsquare[[x,y]]],(postx,posty);
> with(stats);
> eqn:=fit[leastsquare[[x,y]]]([postx,posty]);
> eqn[1];rhs(eqn);
> lineplot:=plot(rhs(eqn),x=1920..1995):lineplot;
> with(plots):
> display(dataplot,lineplot);
> 
> lnPost:=[seq([postage[i][1],ln(postage[i][2])],i=1..nops(postage))];
> dataplot:=plot(lnPost,title="Log cost of a first class stamp",style=point,symbol=box,symbolsize=24):dataplot;
> lnPosty:=[seq(lnPost[i][2],i=1..nops(postage))];
> eqn:=fit[leastsquare[[x,y]]]([postx,lnPosty]);
> lineplot:=plot(rhs(eqn),x=1920..1995):lineplot;
> display(dataplot,lineplot);
> #
> # Part 3 Lagrange Pictures.
> display(contourplot(x^2+y^2,x=-3..3,y=-3..3,contours=[4],scaling=constrained,color=blue),contourplot(x^2+2*y^2,x=-3..3,y=-3..3,contours=[2,4,6,8,10]));
> display(contourplot(x^2+y^2,x=-3..3,y=-3..3,contours=[4],scaling=constrained,color=blue),contourplot(x+2*y,x=-3..3,y=-3..3));
> display(contourplot(x+2*y,x=-3..3,y=-3..3,contours=[2],scaling=constrained,color=blue),contourplot(x^2+y^2,x=-3..3,y=-3..3,contours=30));
> #
> # Part 4 -- Classification of Critical Points
> # From one variable theory
> # f(x) ~ f(a) + f'(a)(x-a) + f''(a)(x-a)^2 + higher order terms
> # Problem is the same if we let a = 0 and f(a) = 0 (translations)
> # At a critical point, f'(a) = 0, so f(x) ~ f''(0) x^2
> # So it is enough to know what a*x^2 does
> # a>0 smiling local min
> # a<0 frowning local max
> # a=0 test fails (you can figure it out from the next non-zero term) but it is `flat'
> #
> # Similarly for two variables, we just need to use the quadratic part of the
> # Taylor series   (f_xx(a,b) (x-a)^2 + 2 f_xy (x-a)(y-b) + f_yy (y-b)^2 )/2
> # Let a = b = 0; We are looking at
> f:=a*x^2+b*x*y+c*y^2;
> Diff(f,x,x)=diff(f,x,x);Diff(f,x,y)=diff(f,x,y);Diff(f,y,y)=diff(f,y,y);
> (diff(f,x,x)*x^2+2*diff(f,x,y)*x*y+diff(f,y,y)*y^2)/2;
> diff(f,x,x)*diff(f,y,y)-(diff(f,x,y))^2;
> with(linalg);
> A:=matrix([[a,b/2],[b/2,c]]);
> 4*det(A);
> #
> contourplot(x^2+2*y^2,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> contourplot(x^2+4*y^2,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> contourplot(x^2-2*y^2,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> contourplot(x^2-4*y^2,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> contourplot(x^2,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> #
> # How does the x*y term change things
> contourplot(x^2+2*y^2-x*y,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> contourplot(x^2+2*y^2+x*y,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> contourplot(x^2-2*y^2-x*y,x=-3..3,y=-3..3,scaling=constrained,title="Where are my contours");
> #
> # somehow x*y rotates things
> #
> # Let's try rotating our co-ordinates
> theta:=Pi/3;theta:='theta';
> u:=vector([cos(theta),sin(theta)]);v:=vector([-sin(theta),cos(theta)]);dp:=dotprod(u,v,orthogonal);
> i:=vector([1,0]);j:=vector([0,1]);null:=vector([0,0]);origin:=[0,0];
> with(plottools);
> display(plot([cos(t),sin(t),t=0..2*Pi]),arrow(null,i,.1,.2,.2,color=red),arrow(null,j,.1,.2,.2,color=red),arrow(null,u,.1,.2,.2,color=blue),arrow(null,v,.1,.2,.2,color=green),title="Rotation",scaling=constrained);
> #f depends on x and y which depend on s and t.
> rot:=matadd(scalarmul(u,s),scalarmul(v,t));
> fst:=subs({x=rot[1],y=rot[2]},f);
> Diff('f(s,t)',s,t)=diff(fst,s,t);
> b*(cos(theta)^2-sin(theta)^2)=(a-c)*(2*sin(theta)*cos(theta));
> expand([cos(2*theta)+(1-cos(theta)^2-sin(theta)^2),sin(2*theta)]);
> b*cos(2*theta)=(a-c)*sin(2*theta);
> tan(2*theta)=b/(a-c);
> # Example
> g:=x^2+2*y^2-x*y;
> solve(tan(2*theta)=-1/(1-2));
> theta:=%;
> u:=vector([cos(theta),sin(theta)]);v:=vector([-sin(theta),cos(theta)]);rot:=matadd(scalarmul(u,s),scalarmul(v,t));
> gst:=subs({x=rot[1],y=rot[2]},g);
> coeff(coeff(expand(gst),s),t);
> cos(Pi/8):=solve(x^2=(1+cos(Pi/4))/2,x)[1];
> sin(Pi/8):=solve(x^2=(1-cos(Pi/4))/2,x)[1];
> coeff(coeff(expand(gst),s),t);
> radnormal(coeff(coeff(expand(gst),s),t));
> #
> # Part 5 -- What does Diff(f,x,y) measure
> #
> Diff(F,x)=Limit((F(x+h)-F(x))/h,h=0);Diff(F,x)=Limit((F(x+h)-F(x-h))/(2*h),h=0);
> Diff(F,x)=(F(x+h)-F(x-h))/(2*h);# for small h
> Diff(F,x,x)=simplify(( ( F(x+h+h)-F(x+h-h) )/(2*h) )-(F(x-h+h)-F(x-h-h))/(2*h))/(2*h);
> Diff(F,x,y)=(Diff(F,x)(x,y+h)-Diff(F,x)(x,y-h))/(2*h);
> Diff(F,x,y)=simplify( (F(x+h,y+h)-F(x-h,y+h))/(2*h) - (F(x+h,y-h)-F(x-h,y-h))/(2*h) )/(2*h);
> # Look at the lattice pattern
> # Look at the solutions to the PDE Diff(F,x,y) = 0, namely f(x) + g(y)
>