###################################################################### ##Math421: Save this file as Math421. To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read Math421 : # ##(in Windows (DOS) you need to type: # ##read `C:\\direct1\\direct2\\Math421`: # #assuming that Math421 resides in C:\direct1\direct2 ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg@math.rutgers.edu. # ####################################################################### #Created: Nov. 25, 2002 #This version: Dec. 5, 2002 (More Heat Eq. stuff and ODE) #Previous version: Dec. 2, 2002 (Thanks to Philip Engel!) #Math421: A Maple package to accopmany Math421 #Please report bugs to zeilberg@math.rutgers.edu print(`To save your work in dumb terminals use the command`): print(`writeto; Use the Maple help to find out about it`): print(``): print(`Created: Nov. 25, 2002`): print(`Version of Dec. 5, 2002`): print(`Previous version: Dec. 2, 2002 `): print(`(Thanks to Philip Engel, for correcting)`): print(` numerous bugs and erros!`): print(`This version of Dec. 5 also has several kinds of Heat Eqs`): print(`And Stability and Nature of Critical Points (sections 10.3, 10.4)`): print(`The exercises refer to the textbook, Advanced Engineering Mathematics`): print(`by Peter V. O'Neil, fifth edition `): lprint(``): print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`): lprint(``): print(`Please report bugs to zeilberg@math.rutgers.edu`): lprint(``): print(`The most current version of this package `): print(` is available from`): print(`http://www.math.rutgers.edu/~zeilberg`): print(`For a list of the procedures type Help(), for help with`): print(`a specific procedure, type Help(procedure_name)`): print(`Math421, a Maple package to accompany Math421 Dr. Z.'s class`): print(`This is for MANDATORY self-study computer-labs`): print(`To replace Thu. Nov. 19`): print(`MANDATORY HOMEWORK: For each of the procedures listed`): print(`by typing Help(); find TWO problems`): print(`from ALL the homework assigned in this class`): print(`that can be answered by that procedure.`): Help:=proc() if args=NULL then print(`Contains the following procedures: Critical,`): print(` dAlembert, FourierSeries `): print(` FourierSineSeries, FourierCosineSeries , `): print(` HeatEqDifferentTemp, HeatEqInsulatedEnds, `): print(`InHomogWaveEq, WaveEq, WaveEqZeroDisp, WaveEqZeroInitVel `): fi: if nops([args])=1 and op(1,[args])=Critical then print(`Critical(a,b,c,d): Given a 2 by 2 matrix A=[[a,b],[c,d]]`): print(`finds the nature and stability and asympt. stability of`): print(`the system x'(t)=a*x(t)+b*y(t), y'(t)=c*x(t)+d*y(t)`): print(`For example, to do ex. 10.3.6 (and 10.4.6) do`): print(`Critical(2,-7,5,-10):`): fi: if nops([args])=1 and op(1,[args])=dAlembert then print(`dAlembert(f,g,x,t,c): d'Alembert's solution of the Wave Eq`): print(`u_tt=c^2 u_xx for -infinity0`): print(`with u(x,0)=f, u_t(x,0)=g(x) -infinity0) with u(0,t)=T1, u(L,t)=T2 (t>=0), `): print(` u(x,0)=f(x) for 0<=x<=L. For example to do `): print( ` Ex. 17.2.17, type : HeatEqDifferentTemp(x^2,16,x,t,1,2,5); `): fi: if nops([args])=1 and op(1,[args])=HeatEqInsulatedEnds then print(`HeatEqInsulatedEnds(f,k,x,t,L): The solution of the Heat equation`): print(`u_t=k u_xx (00) with u_x(0,t)=u_x(L,t)=0 (t>=0), `): print(`u(x,0)=f(x) for 0<=x<=L. For example to do`): print(`Ex. 17.2.4, type : HeatEqInsulatedEnds(sin(x),1,x,t,Pi);`): fi: if nops([args])=1 and op(1,[args])=HeatEqZeroTemp then print(`HeatEqZeroTemp(f,k,x,t,L): The solution of the Heat equation`): print(`u_t=k u_xx (00) with u(0,t)=u(L,t)=0 (t>=0), `): print(`u(x,0)=f(x) for 0<=x<=L. For example to do`): print(`Ex. 17.2.1, type : HeatEqZeroTemp(x*(L-x),k,x,t,L);`): fi: if nops([args])=1 and op(1,[args])=InHomogWaveEq then print(`InHomogWave(F,f,g,x,t,c): Solution of`): print(`u_tt(x,t)=c^2 u_xx (x,t)+ F(x,t) for -infinity0`): print(`with u(x,0)=f, u_t(x,0)=g(x) -infinity0 then gu:=[op(gu),[sol[i],limit(A,n=sol[i])]]: fi: od: gu: end: #FourierSeries(f,x,L): The Fourier series of a function #of x given by the expression f, over the interval [-L,L] #For example try FourierSeries(x^2,x,Pi); FourierSeries:=proc(f,x,L) local a0,an,bn,n,F1,g,f1,F2: F1:=ExtractSines(f,x,L): F2:=ExtractCosines(f,x,L): g:=F1[1]+F2[1]: f1:=expand(f-g): assume(n,integer): a0:=(1/L)*simplify(int(f1,x=-L..L)): an:=simplify((1/L)*simplify(int(f1*cos(n*\Pi*x/L),x=-L..L))): bn:=simplify((1/L)*simplify(int(f1*sin(n*\Pi*x/L),x=-L..L))): if an=0 and bn=0 then RETURN(a0/2+g): elif an<>0 and bn=0 then RETURN(a0/2+g+Sum(an*cos(n*Pi*x/L),n=1..infinity)): elif an=0 and bn<>0 then RETURN(g+1/2*a0+Sum(bn*sin(n*Pi*x/L),n=1..infinity)): else RETURN(g+1/2*a0+ Sum(an*cos(n*Pi*x/L),n=1..infinity)+Sum(bn*sin(n*Pi*x/L),n=1..infinity)): fi: end: #FourierSineSeries(f,x,L): The Fourier Sine series of a function #of x given by the expression f, over the interval [0,L] #For example try FourierSineSeries(x^2,x,Pi); FourierSineSeries:=proc(f,x,L) local bn,n,F,g,f1: F:=ExtractSines(f,x,L): g:=F[1]: f1:=F[2]: assume(n,integer): bn:=2/L*int(f1*sin(n*\Pi*x/L),x=0..L): bn:=simplify(bn): if bn=0 then RETURN(g): else RETURN(g+Sum(bn*sin(n*Pi*x/L),n=1..infinity)): fi: end: #FourierCosineSeries(f,x,L): The Fourier Cosine series of a function #of x given by the expression f, over the interval [0,L] #For example try FourierCosineSeries(x^2,x,Pi); FourierCosineSeries:=proc(f,x,L) local a0, an,n,F, g,f1: F:=ExtractCosines(f,x,L): g:=F[1]: f1:=F[2]: a0:=2/L*int(f1,x=0..L): assume(n,integer): an:=2/L*int(f1*cos(n*\Pi*x/L),x=0..L): an:=simplify(an): if an<>0 then RETURN(a0/2+g+Sum(an*cos(n*Pi*x/L),n=1..infinity)): else RETURN(a0/2+g): fi: end: #ExtractSines(f,x,L): Given an expression with sin(n*x*Pi/L) in it #extracts those, followed by the rest ExtractSines:=proc(f,x,L) local i,f1,f11,f111,g,c,F: f1:=f: F:=f1: g:=0: if type(f1,`+`) then for i from 1 to nops(f1) do f11:=op(i,f1): if type(f11,`*`) then c:=op(1,f11): f111:=op(2,f11): if type(c,numeric) and type(f111,function) then if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: if type(f11,function) then c:=1: f111:=f11: if type(c,numeric) and type(f111,function) then if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: od: fi: if not type(f1,`+`) then f11:=f1: if type(f11,`*`) then c:=op(1,f11): f111:=op(2,f11): if type(c,numeric) and type(f111,function) then if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: if type(f11,function) then c:=1: f111:=f11: if type(c,numeric) and type(f111,function) then if op(0,f111)=sin and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: fi: g,F: end: #ExtractCosines(f,x,L): Given an expression with cos(n*x*Pi/L) in it #extracts those, followed by the rest ExtractCosines:=proc(f,x,L) local i,f1,f11,f111,g,c,F: f1:=f: F:=f1: g:=0: if type(f1,`+`) then for i from 1 to nops(f1) do f11:=op(i,f1): if type(f11,`*`) then c:=op(1,f11): f111:=op(2,f11): if type(c,numeric) and type(f111,function) then if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: if type(f11,function) then c:=1: f111:=f11: if type(c,numeric) and type(f111,function) then if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: od: fi: if not type(f1,`+`) then f11:=f1: if type(f11,`*`) then c:=op(1,f11): f111:=op(2,f11): if type(c,numeric) and type(f111,function) then if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: if type(f11,function) then c:=1: f111:=f11: if type(c,numeric) and type(f111,function) then if op(0,f111)=cos and type(op(1,f111)/(Pi*x/L),integer) then g:=g+f11: F:=F-f11: fi: fi: fi: fi: g,F: end: #ExtractSinesWithCosines(f,x,t,c,L): #Given an expression with sin(n*x*Pi/L) in it #extracts those with cosines attached, followed by the rest #(without the cosines) ExtractSinesWithCosines:=proc(f,x,t,c,L) local i,f1,f11,f111,g,c1,F,mumu: f1:=f: F:=f1: g:=0: if type(f1,`+`) then for i from 1 to nops(f1) do f11:=op(i,f1): if type(f11,`*`) then c1:=op(1,f11): f111:=op(2,f11): if type(c1,numeric) and type(f111,function) then mumu:=op(1,f111)/(Pi*x/L) : if op(0,f111)=sin and type(mumu,integer) then g:=g+f11*cos(mumu*Pi*c*t/L): F:=F-f11: fi: fi: fi: if type(f11,function) then c1:=1: f111:=f11: mumu:=op(1,f111)/(Pi*x/L) : if type(c1,numeric) and type(f111,function) then if op(0,f111)=sin and type(mumu,integer) then g:=g+f11*cos(mumu*Pi*c*t/L): F:=F-f11: fi: fi: fi: od: fi: if not type(f1,`+`) then f11:=f1: if type(f11,`*`) then c1:=op(1,f11): f111:=op(2,f11): if type(c1,numeric) and type(f111,function) then mumu:=op(1,f111)/(Pi*x/L): if op(0,f111)=sin and type(mumu,integer) then g:=g+f11*cos(mumu*Pi*c*t/L): F:=F-f11: fi: fi: fi: if type(f11,function) then c1:=1: f111:=f11: if type(c1,numeric) and type(f111,function) then mumu:=op(1,f111)/(Pi*x/L): if op(0,f111)=sin and type(mumu,integer) then g:=g+f11*cos(mumu*Pi*c*t/L): F:=F-f11: fi: fi: fi: fi: g,F: end: #WaveEqZeroInitVel(f,c,x,t,L): The solution of the wave equation #y_tt=c^2 y_xx in [0,L] with y(0,t)=y(L,t)=0, y(x,0)=f(x) and #y_t(x,0)=0 (00 #with u(x,0)=f, u_t(x,0)=g(x) -infinity0 #with u(x,0)=f, u_t(x,0)=g(x) -infinity0) with u(0,t)=u(L,t)=0 (t>=0), #u(x,0)=f(x) for 0<=x<=L. For example to do #Ex. 17.2.1, type : HeatEqZeroTemp(x*(L-x),k,x,t,L); HeatEqZeroTemp:=proc(f,k,x,t,L) local cn,n,F,g1,f1: F:=ExtractSinesWithExp(f,x,t,k,L): g1:=F[1]: f1:=F[2]: assume(n,integer): cn:=(2/L)*int(f1*sin(n*\Pi*x/L),x=0..L): cn:=simplify(cn): if cn=0 then RETURN(g1): else RETURN(g1+Sum(cn*sin(n*Pi*x/L)*exp(-(n*Pi/L)^2*k*t),n=1..infinity)): fi: end: #ExtractCosinesWithExp(f,x,t,k,L): #Given an expression with cos(n*x*Pi/L) in it #extracts those with exp(-k*(n*Pi/L)^2*t) attached, followed by the rest #(without the sines) ExtractCosinesWithExp:=proc(f,x,t,k,L) local i,f1,f11,f111,g,c1,F,mumu: f1:=f: F:=f1: g:=0: if type(f1,`+`) then for i from 1 to nops(f1) do f11:=op(i,f1): if type(f11,`*`) then c1:=op(1,f11): f111:=op(2,f11): if type(c1,numeric) and type(f111,function) then mumu:=op(1,f111)/(Pi*x/L) : if op(0,f111)=cos and type(mumu,integer) then g:=g+f11*exp(-(mumu*Pi/L)^2*k*t): F:=F-f11: fi: fi: fi: if type(f11,function) then c1:=1: f111:=f11: mumu:=op(1,f111)/(Pi*x/L) : if type(c1,numeric) and type(f111,function) then if op(0,f111)=cos and type(mumu,integer) then g:=g+f11*exp(-(mumu*Pi/L)^2*k*t): F:=F-f11: fi: fi: fi: od: fi: if not type(f1,`+`) then f11:=f1: if type(f11,`*`) then c1:=op(1,f11): f111:=op(2,f11): if type(c1,numeric) and type(f111,function) then mumu:=op(1,f111)/(Pi*x/L): if op(0,f111)=cos and type(mumu,integer) then g:=g+f11*exp(-(mumu*Pi/L)^2*k*t): F:=F-f11: fi: fi: fi: if type(f11,function) then c1:=1: f111:=f11: if type(c1,numeric) and type(f111,function) then mumu:=op(1,f111)/(Pi*x/L): if op(0,f111)=cos and type(mumu,integer) then g:=g+f11*exp(-(mumu*Pi/L)^2*k*t): F:=F-f11: fi: fi: fi: fi: g,F: end: #HeatEqInsulatedEnds(f,k,x,t,L): The solution of the Heat equation #u_t=k u_xx (00) with u_x(0,t)=u_x(L,t)=0 (t>=0), #u(x,0)=f(x) for 0<=x<=L. For example to do #Ex. 17.2.4, type : HeatEqInsulatedEnds(sin(x),1,x,t,Pi); HeatEqInsulatedEnds:=proc(f,k,x,t,L) local cn,n,F,g1,f1,c0: F:=ExtractCosinesWithExp(f,x,t,k,L): g1:=F[1]: f1:=F[2]: c0:=(2/L)*int(f1,x=0..L): assume(n,integer): cn:=(2/L)*int(f1*cos(n*\Pi*x/L),x=0..L): cn:=simplify(cn): if cn=0 then RETURN(g1): else RETURN(g1+1/2*c0+Sum(cn*cos(n*Pi*x/L)*exp(-(n*Pi/L)^2*k*t),n=1..infinity)): fi: end: #HeatEqDifferentTemp(f,k,x,t,L,T1,T2): The solution of the Heat equation #u_t=k u_xx (00) with u(0,t)=T1, u(L,t)=T2 (t>=0), #u(x,0)=f(x) for 0<=x<=L. For example to do #Ex. 17.2.17, type : HeatEqDifferentTemp(x^2,16,x,t,1,2,5); HeatEqDifferentTemp:=proc(f,k,x,t,L,T1,T2) local psi,f1: psi:=(T2-T1)/L*x+T1: f1:=f-psi: HeatEqZeroTemp(f1,k,x,t,L)+psi: end: #Critical(a,b,c,d): Given a 2 by 2 matrix A=[[a,b],[c,d]] #finds the nature and stability and asympt. stability of #the system x'(t)=a*x(t)+b*y(t), y'(t)=c*x(t)+d*y(t) #For example, to do ex. 10.3.6 (and 10.4.6) do #Critical(2,-7,5,-10): Critical:=proc(a,b,c,d) local lam,eq,gu: eq:=expand((a-lam)*(d-lam)-b*c): gu:=evalf(evalc({solve(eq,lam)})): if nops(gu)=1 then lam:=evalf(gu[1]): if a=lam and d=lam and b=0 and c=0 then if lam<0 then RETURN([proper_node,sink, stable,asymptotically_stable]): elif lam>0 then RETURN([proper_node,source, unstable]): fi: else if lam<0 then RETURN([improper_node,sink, stable,asymptotically_stable]): elif lam>0 then RETURN([improper_node,source, unstable]): fi: fi: fi: if nops(gu)=2 then if coeff(gu[1],I,1)=0 and evalf(gu[1]*gu[2])>0 then if gu[1]<0 then RETURN([nodal_sink,stable, asymptotically_stable]): else RETURN([nodal_source,unstable]): fi: elif coeff(gu[1],I,1)=0 and evalf(gu[1]*gu[2])<0 then RETURN([saddle_point,unstable]): elif coeff(gu[1],I,1)<>0 and evalf(gu[1]-I*coeff(gu[1],I,1))>0 then RETURN([spiral_source,unstable]): elif coeff(gu[1],I,1)<>0 and evalf(gu[1]-I*coeff(gu[1],I,1))<0 then RETURN([spiral_sink,stable,asymptotically_stable]): elif coeff(gu[1],I,1)<>0 and evalf(gu[1]-I*coeff(gu[1],I,1))=0 then RETURN([center,stable,not_asymptotically_stable]): else ERROR(`Something is wrong`): fi: fi: end: