/* wxMaxima commands */ /* NOTE: after cutting-and-pasting into wxMaxima, use CTRL-RET to execute the command. */ /* (These three lines are each enclosed by Maxima comment brackets.) */ /* Inspired by Chapter 1 of Hugh L. Montgomery's "Early Fourier Analysis" */ /* M.V.Wickerhauser Math 410 lecture notes 2017-08-30 2017-09-01 2017-09-06 */ /**** Basic commands to perform complex-number arithmetic */ /* Square root of -1 is "%i". When squared it gives -1: */ %i; %i^2; /* Define a complex number z with real part x, imaginary part y */ declare(z, complex); declare(x,real,y,real); z:x+y*%i; /* colon ":" is the assignment operator in Maxima */ realpart(z); imagpart(z); cabs(z); /* complex absolute value */ carg(z); /* complex argument */ polarform(z); assume(x>0, y>0); carg(z); /* simplifying assumption */ polarform(z); conjugate(z); /* complex conjugate */ z*conjugate(z); expand(z*conjugate(z)); /* Define another complex number */ declare(w,complex,u,real,v,real); w:u+v*%i; realpart(z*w); /* Real and imaginary parts */ imagpart(z*w); c: 17 * %e ^ (%i * %pi /3); /* some complex constant in polar form */ polarform(c); /* force display in polar form */ rectform(polarform(c)); /* revert to rectangular form... */ realpart(c); /* ...with real part */ imagpart(c); /* ......and imaginary part */ /**** Matrix form of complex numbers */ I: matrix([1,0],[0,1]); J: matrix([0,1],[-1,0]); Z:x*I+y*J; /* Use capital "Z" to be distinct from z=x+%i*y */ invert(Z); /* matrix form of 1/z */ determinant(Z); /* equals |z|^2 = x^2+y^2 */ adjoint(Z); /* gives matrix form of conjugate(z) */ W:u*I+v*J; Z.W; /* matrix multiplication symbol is "." */ /* (Unfortunately, "*" in Maxima gives componentwise multiplication.) */ /**** Prove that the diagonals of a rhombus (0,a,b,a+b) are perpendicular */ declare(r,real, theta,real); assume(r>0); assume(theta>0, theta<%pi/2); /* insures acute angle at 0 */ declare(a,complex, b,complex); declare(diag1, complex, diag2, complex); /* Vertices at 0, a on real axis, b and a+b in first quadrant */ a:r; b:r*exp(%i*theta); diag1:a+b; diag2:b-a; /* It suffices to show that arg(diag1/diag2) is in {pi/2, 3pi/2} */ carg(diag1/diag2); /* not useful */ /* ... or equivalently, that Re(diag1/diag2)=0 */ realpart(diag1/diag2); expand(realpart(diag1/diag2)); radcan(expand(realpart(diag1/diag2))); /* "RADical CANonical" simplification */ /* Note that the numerator is zero since sin(theta)^2+cos(theta)^2=1. */ /**** Prove Napoleon's Theorem */ /* Start with a triangle abc placed canonically with base ab along the real axis, vertex a at the origin, and an acute angle cab so that c lies in the first quadrant: */ declare(p,real, q,real); assume(p>0,q>0); /* side lengths p=|b-a|, q=|c-a| */ declare(a,complex, b,complex, c,complex); /* vertices in the complex plane */ declare(theta,real); assume(theta>0,theta<%pi/2); a:0; b:p; c:q*exp(%i*theta); /* canonical placement */ /* Note: if p,q,theta are specified numerically, then the triangle can be plotted for visualization: */ eg_triangle: subst( [p=3, q=2, theta=1.0], /* specific numerical values */ [ [realpart(a),imagpart(a)], [realpart(b),imagpart(b)], [realpart(c),imagpart(c)], [realpart(a),imagpart(a)] ] /* list of [x,y] pairs... */ ); /* ...is now a list of numerical [x,y] pairs */ /* (Vertex "a" appears first and last to close the triangle.) */ plot2d([ discrete, eg_triangle] ); /* Find A,B,C */ rot60: exp(%i*%pi/3); /* rotation by +60 degrees */ A: c+(b-c)*rot60; B: a+(c-a)*rot60; C: b+(a-b)*rot60; /* Find alpha,beta,gamma, the barycenters of the equilateral triangles A,B,C */ alpha: (A+b+c)/3; beta: (a+B+c)/3; gamma: (a+b+C)/3; /* Check the lengths of the sides, which should be equal */ ab: beta-alpha; bc: gamma-beta; ca: alpha-gamma; expand(cabs(ab)); /* expand() actually simplifies! */ expand(cabs(bc)); expand(cabs(ca)); /* ...and these three are evidently equal */ /* (Try checking the lengths without expand() and be disappointed.) */ /* Checking the angles at alpha, beta, and gamma is also disappointing */ carg(ca/ab); /* should be pi/3 */ carg(ab/bc); /* should be pi/3 */ carg(bc/ca); /* should be pi/3 */ /**** Polynomials */ P:p3*x^3 + p2*x^2 + p1*x + p0; /* Cubic polynomial, unspecified coefficients */ Q:q2*x^2 + q1*x + q0; /* Quadratic polynomial, unspecified coefficients */ divide(P,Q); /* returns [quotient, remainder] form of P/Q */ PP: subst([p3=1,p2=%pi,p1=%e,p0=0],P); /* specific coefficients for P */ QQ: subst([q2=4,q1=-1,q0=0],Q); /* specific coefficients for Q */ divide(PP,QQ); /* specific example [quotient,remainder] */ gcd(PP,QQ); /* greatest common divisor polynomial */ /**** Power series */ /* Taylor polynomials of some functions */ taylor(sin(x),x,0,10); /* sin(x), in x, about x=0, degree 10 or less */ taylor(cos(x),x,0,10); /* cos(x), in x, about x=0, degree 10 or less */ taylor(exp(%i*x),x,0,10); /* exp(ix), in x, about x=0, degree 10 or less */ /* Check Euler's formula to degree 10 */ taylor(cos(x),x,0,10) + %i*taylor(sin(x),x,0,10) - taylor(exp(%i*x),x,0,10); /* ...same result to degree 1000 */ taylor(cos(x),x,0,1000) + %i*taylor(sin(x),x,0,1000) - taylor(exp(%i*x),x,0,1000); /* Check the double angle formula for sine to degree 100 */ taylor(sin(2*x),x,0,100) - 2*taylor(sin(x),x,0,100)*taylor(cos(x),x,0,100); /* Define a function by its power series (Taylor series about x=0) */ deftaylor(F(x), x^2 + sum(x^i/(2^i*i!^2), i, 4, inf) ); powerseries(F(x),x,0); /* Find the Taylor polynomials (about x=0) of expressions in F: */ F(x)*taylor(exp(x),x,0,5); /**** Analytic functions and domains of analyticity */ taylor(1/(1-x),x,0,5); /* converges for |x|<1, diverges if |x|>1 */ taylor(1/(1-x),x,-1/2,5); /* R=3/2; converges if |x+1/2|<3/2 */ taylor(1/(1-x),x,%i,5); /* R=sqrt(2); converges if |x-i|