file mbe10code.txt has Maxima code for Maxima by Example, Chapter 10, Fourier Series, Fourier and Laplace Transforms. Edwin L Woollett, 2009 woollett@charter.net http://www.csulb.edu/~woollett 10.1 Introduction 10.2 Fourier Series Expansion of a Function 10.2.1 Expansion of a Function over (-pi, pi) 10.2.2 Fourier Expansion of f(x) = x over (-pi, pi) (%i1) (declare(n,integer), assume(n > 0 ), facts() ); (%i2) define(b(n),integrate(x*sin(n*x),x,-%pi,%pi)/%pi ); (%i3) map('b,makelist(i,i,1,7)); (%i4) fs(nmax) := sum( b(m)*sin(m*x),m,1,nmax)$ (%i5) map('fs,[1,2,3,4] ); (%i6) (load(draw), load(qdraw) )$ (%i7) qdraw( xr(-5.6,5.6),yr(-4,4), ex([x,fs(1),fs(2)],x,-%pi,%pi),key(bottom) )$ (%i8) qdraw( xr(-5.6,5.6),yr(-4,4), ex([x,fs(1),fs(2),fs(3),fs(4)],x,-%pi,%pi),key(bottom) )$ 10.2.3 The calculus/fourie.mac Package: fourier, foursimp, fourexpand (%i1) facts(); (%i2) f(m) := block([n],declare(n,integer),assume( n > 0 ), if m < 2 then n :2 else n:3,(2*n*m) )$ (%i3) f(1); (%i4) facts(); (%i5) is(n>0); restart (%i1) facts(); (%i2) (load(fourie), facts() ); (%i3) (declare(n,integer), facts() ); (%i4) clist : fourier(x,x,%pi); (%i7) facts(); (%i8) fs(nmax) := fourexpand(clist,x,%pi, nmax )$ (%i9) map( 'fs, [1,2,3,4] ); (%i10) b(n); (%i11) define(b(n),rhs(%t6) ); (%i12) map( 'b, makelist(i,i,1,7) ); 10.2.4 Fourier Series Expansion of a Function over (-p, p) 10.2.5 Fourier Series Expansion of the Function abs(x) (%i13) integrate(abs(x)*cos(n*%pi*x/2),x,-2,2)/2; restart (%i1) (declare(n,integer),assume(n > 0), facts() ); (%i2) a0 :integrate(-x,x,-2,0)/2 + integrate(x,x,0,2)/2; (%i3) an : integrate((-x)*cos(n*%pi*x/2),x,-2,0)/2 + integrate(x*cos(n*%pi*x/2),x,0,2)/2; (%i4) an : (ratsimp(an), factor(%%) ); (%i5) define(a(n),an); (%i6) map( 'a, [1,2,3,4,5] ); (%i7) fs(nmax) := a0/2 + sum(a(m)*cos(m*%pi*x/2),m,1,nmax)$ (%i8) map('fs, [1, 3, 5] ); (%i9) (load(draw),load(qdraw))$ (%i10) qdraw( ex([abs(x),fs(1)],x,-2,2),key(bottom) )$ (%i11) qdraw( ex([abs(x),fs(1),fs(3)],x,-2,2),key(bottom) )$ (%i12) qdraw( ex([abs(x),fs(5) ],x,-2,2),key(bottom) )$ restart (%i1) ( load(fourie),facts() ); (%i2) (declare(n,integer),facts()); (%i3) fourier(abs(x),x,2); (%i6) clist : foursimp(%); (%i9) facts(); (%i10) fs(nmax) := fourexpand(clist,x,2,nmax )$ (%i11) map( 'fs, [1,3,5] ); 10.2.6 Fourier Series Expansion of a Rectangular Pulse (%i1) f(x):= if x >= -1 and x <= 1 then 3/2 else 0$ (%i2) map('f, [-3/2,-1,0,1,3/2] ); (%i3) (load(draw),load(qdraw) )$ (%i4) qdraw( yr(-0.5,2), ex1(f(x),x,-2,2,lw(5),lc(blue) ) )$ (%i5) integrate(f(x),x,-2,2); (%i6) (declare(n,integer), assume(n>0),facts() ); (%i7) a0 : (1/2)*integrate( (3/2),x,-1,1 ); (%i8) define(a(n),(1/2)*integrate((3/2)*cos(n*%pi*x/2),x,-1,1)); (%i9) map( 'a, makelist(i,i,1,7) ); (%i10) fs(nmax) := a0/2 + sum( a(m)*cos(m*%pi*x/2),m,1,nmax )$ (%i11) map( 'fs, [1,3] ); (%i12) qdraw( yr(-0.5,2),ex([f(x),fs(1),fs(3) ],x,-2,2) )$ (%i13) qdraw( yr(-0.5,2),ex([f(x),fs(11) ],x,-2,2) )$ 10.2.7 Fourier Series Expansion of a Two Element Pulse (%i1) f(x):= if x >= -5 and x < 0 then -5 elseif x >= 0 and x <= 5 then 5 else 0$ (%i2) map('f,[-6,-5,-1,0,1,5,6]); (%i3) ( load(draw),load(qdraw) )$ (%i4) qdraw( yr(-8,8), ex1(f(x),x,-10,10,lw(5),lc(blue) ) )$ (%i5) a0 : (1/10)*( integrate( -5, x, -5, 0 ) + integrate( 5, x, 0, 5 ) ); (%i6) an : (1/10)*(integrate( -5*cos( n*%pi*x/10 ), x, -5, 0 ) + integrate(5*cos(n*%pi*x/10), x, 0, 5 ) ); (%i7) bn : ( (1/10)*(integrate( -5*sin(n*%pi*x/10), x, -5, 0 ) + integrate( 5*sin(n*%pi*x/10), x, 0, 5 ) ), ratsimp(%%) ); (%i8) define( b(n), bn ); (%i9) map('b,makelist(i,i,1,7)); (%i10) fs(nmax) := sum( b(m)*sin(m*%pi*x/10), m, 1, nmax )$ (%i11) map('fs,[1,2,3]); (%i12) qdraw( xr(-15, 15), yr(-10, 10), ex( [f(x), fs(1), fs(2) ], x, -10, 10 ) )$ (%i13) qdraw( xr(-15, 15), yr(-10, 10), ex( [ f(x), fs(11) ], x, -10, 10 ) )$ (%i14) qdraw(yr(-10, 10),ex( fs(11), x, -10 , 30 ) )$ 10.2.8 Exponential Form of a Fourier Series 10.3 Fourier Integral Transform Pairs 10.3.1 Fourier Cosine Integrals and fourintcos(..) (%i1) f:sin(x)*exp(-x)$ (%i2) integrate(f*cos(w*x),x,0,inf); (%i3) fcw : (2/%pi)*ratsimp(%); (%i4) integrate(fcw*cos(w*x),w,0,inf); (%i5) load(fourie); (%i6) fourintcos(f,x); (%i7) az : ratsimp(rhs(%t6)); (%i8) (2/%pi)*ratsimp(%pi*az/2); 10.3.2 Fourier Sine Integrals and fourintsin(..) (%i1) f:cos(x)*exp(-x)$ (%i2) assume(w>0)$ (%i3) integrate(f*sin(w*x),x,0,inf); (%i4) fsw : (2/%pi)*%; (%i5) integrate(fsw*sin(w*x),w,0,inf); (%i6) load(fourie); (%i7) facts(); (%i8) (forget(w>0),facts()); (%i9) fourintsin(f,x); (%i10) bz : rhs(%t9); 10.3.3 Exponential Fourier Integrals and fourint(..) 10.3.4 Example 1: Even Function (%i1) assume(w>0)$ (%i2) i1:integrate(exp(%i*w*x)*cos(x)*exp(x),x,minf,0); (%i3) i2:integrate(exp(%i*w*x)*cos(x)*exp(-x),x,0,inf); (%i4) i12:ratsimp(i1+i2); (%i5) 2*ratsimp(i12/2); (%i6) iexp:%/(2*%pi); (%i7) integrate(exp(-%i*w*x)*iexp,w,minf,inf); (%i8) integrate(exp(-%i*w*x)*iexp,w,minf,inf); (%i9) i3:ratsimp(integrate(cos(x)*exp(-x)*cos(w*x),x,0,inf)); (%i10) i3:(2/%pi)*i3; restart (%i1) load(fourie); (%i2) fourint(cos(x)*exp(-abs(x)),x); (%i4) ratsimp(rhs(%t2)); (%i5) (2/%pi)*ratsimp(%pi*%/2); 10.3.5 Example 2: Odd Function (%i1) ( assume(w>0), facts()); (%i2) i1:ratsimp(integrate(exp(%i*w*x)*sin(x)*exp(x),x,minf,0)); (%i3) i2:ratsimp(integrate(exp(%i*w*x)*sin(x)*exp(-x),x,0,inf)); (%i4) iexp:ratsimp(i1+i2)/(2*%pi); (%i5) facts(); (%i6) integrate(exp(-%i*w*x)*iexp,w,minf,inf); (%i7) integrate(exp(-%i*w*x)*iexp,w,minf,inf); (%i8) facts(); (%i9) ratsimp(integrate(sin(x)*exp(-x)*sin(w*x),x,0,inf)); (%i10) (2/%pi)*%; (%i11) load(fourie); (%i12) fourint(sin(x)*exp(-abs(x)),x); 10.3.6 Example 3: A Function Which is Neither Even nor Odd (%i1) ( assume(w>0), facts()); (%i2) i1:ratsimp(integrate(exp(%i*w*x)*cos(x-1)^2*exp(x),x,minf,0)); (%i3) i2:ratsimp(integrate(exp(%i*w*x)*cos(x-1)^2*exp(-x),x,0,inf)); (%i4) i12 : rectform(ratsimp(i1 + i2)); (%i5) i12 : map('ratsimp,i12); (%i6) i12 : realpart(i12)/(2*%pi) + %i*imagpart(i12)/(2*%pi); (%i7) integrate(exp(-%i*w*x)*i12,w,minf,inf); (%i8) cos(x-1)^2; (%i9) trigreduce(%); (%i10) ratsimp(%); (%i11) trigexpand(%); (%i12) trigreduce(%); (%i13) factor(%); (%i14) trigexpand(%); (%i15) trigsimp(%); (%i16) load(fourie); (%i17) fourint(cos(x-1)^2*exp(-abs(x)),x); (%i19) az : rhs(%t17); (%i20) bz : rhs(%t18); (%i21) iexp_f : az/2 + %i*bz/2; (%i22) subst(z=w,iexp_f) - i12; 10.3.7 Dirac Delta Function delta(x) 10.3.8 Laplace Transform of the Delta Function Using a Limit Method (%i1) assume(s>0,e>0,t0>0)$ (%i2) i1: integrate(exp(-s*t),t,t0-e,t0+e)/(2*e); (%i3) limit(i1,e,0,plus); 10.4 Laplace Transform and Inverse Transform Integrals 10.4.1 Laplace Transform Integrals: laplace(..), specint(..) Comparison of laplace and specint (%i1) assume( s>0, a > 0, b > 0, s > a, s > b )$ (%i2) laplace(1,t,s); (%i3) specint(exp(-s*t),t); (%i4) ilt(%,s,t); (%i5) laplace(t,t,s); (%i6) specint(exp(-s*t)*t,t); (%i7) ilt(%,s,t); (%i8) laplace(exp(a*t),t,s); (%i9) specint(exp(-s*t)*exp(a*t),t); (%i10) laplace(sin(a*t)/a,t,s); (%i11) (specint(exp(-s*t)*sin(a*t)/a,t),ratsimp(%%) ); (%i12) laplace(cos(a*t),t,s); (%i13) (specint(exp(-s*t)*cos(a*t),t ), ratsimp(%%) ); (%i14) laplace(sin(a*t)*t/(2*a),t,s); (%i15) (specint(exp(-s*t)*sin(a*t)*t/(2*a),t ),ratsimp(%%) ); (%i16) map('factorsum,%); (%i17) laplace(exp(a*t)*cos(b*t),t,s); (%i18) map('factorsum,%); (%i19) (specint(exp(-s*t)*exp(a*t)*cos(b*t),t ),ratsimp(%%) ); (%i20) map('factorsum,%); (%i21) expr : t^(1/2) * bessel_j(1, 2 * a^(1/2) * t^(1/2)); (%o21) bessel_j(1, 2 sqrt(a) sqrt(t)) sqrt(t) (%i22) laplace(expr,t,s); (%i23) specint(exp(-s*t)*expr,t); (%i24) ilt(%,s,t); (%i25) assume(s>0)$ (%i26) laplace(erf(sqrt(t)),t,s); (%i27) specint(exp(-s*t)*erf(sqrt(t)),t); (%i28) radcan(%); (%i29) ilt(%,s,t); (%i30) laplace(erf(t),t,s); (%i31) ilt(%,s,t); (%i32) specint(exp(-s*t)*erf(t),t); 10.4.2 Use of the Dirac Delta Function (Unit Impulse Function) delta with laplace(..) (%i1) laplace(delta(t),t,s); (%i2) ilt(1,s,t); (%i3) laplace(delta(t-a),t,s); (%i4) laplace(delta(t-a)*sin(b*t),t,s); 10.4.3 The Inverse Laplace Transform: ilt(..), residue(..) Inverse Laplace Transform: ilt(...) residue(...) (%i1) residue (s/(s^2+a^2), s, a*%i); (%i2) residue (sin(a*x)/x^4, x, 0); (%i3) fs : -s/(s^3 +4*s^2 + 5*s +2)$ (%i4) ft : ( ilt(fs,s,t),collectterms(%%,exp(-t),exp(-2*t) ) ); (%i5) (laplace(ft,t,s), ratsimp(%%) ); (%i6) fs - %; (%i7) fst : exp(s*t)*fs; (%i8) partfrac(fst,s); (%i9) assume( t > 0 )$ (%i10) r1 : residue(fst,s,-2); (%i11) r2 : residue(fst,s,-1); (%i12) r1 + r2 - ft;