/* Fourier analysis Maxima package Authors: Jose A. Vallejo Faculty of Sciences Universidad Autonoma de San Luis Potosi (Mexico) http://galia.fc.uaslp.mx/~jvallejo Emmanuel Roque Faculty of Sciences Universidad Autonoma de San Luis Potosi (Mexico) */ load(fourie)$ load(simplify_sum)$ load(draw)$ /*The following rules are intented to work with intervals in "a canonical way" leq,geq,lss,grtr are for bounded intervals Note: Instead of using constantp inside of matchdeclare, numberp but this approach won't work if %pi is used, so probably is not an option since this is a Fourier package and %pi will appear most of the time. Another possible approach is to use freeof(var) instead of constantp */ matchdeclare(constn, constantp)$ defrule(leq,constn>=xx,xx<=constn)$ defrule(geq,xx>=constn,constn<=xx)$ defrule(lss,constn>xx,xx<constn)$ defrule(grtr,xx>constn,constn<xx)$ /* fun_parts(expr,var) Input: piecewise defined function written in one of the following ways if x>=a_0 and x<=a_1 then expr1 elseif x>a_1 and x<= a_2 then expr2 ... elseif x>a_n and x<=a_{n+1} then expr_{n+1} if x<=a_0 then expr0 elseif x>a0 and x<=a1 then expr1 ... elseif x>=a_n then exprn Output: A list [[a0,a1,expr1],[a1,a2,expr2,]...] Important notes: -In all cases it is expected that a_i<a_{i+1} Writing (x<=7 and x>=6) instead of (x>=6 and x<=7) does not work. -If in the function definition some of the intervals are empty (e.g x>=3 and x<=-3) returns error. -Use of else is currently unsupported, the expr after else is ignored! Use elseif instead. -functionality with logical operator "or" is currently unsupported, will be treated same as "and". -fun_parts does not check if the intervals are disjoint! */ fun_parts(expr,var):=block( [subint,subval,tmp,ll,tmp1,tmp2,xx,ans,tmp3,tmpb,leftU,rightU], if not piecewisep(expr) then error("The function is not piecewise defined") else ( xx:var, ll:(length(expr)-2)/2, tmp:makelist(inpart(expr,i),i,makelist(2*k-1,k,1,ll)), subval:makelist(inpart(expr,i),i,makelist(2*k,k,1,ll)), tmp1:makelist(operatorp(tmp[j],["<",">","<=",">="]),j,1,ll), tmp2:sublist_indices(tmp1,lambda([x],x=true)), /*if tmp2 is an empty list then all the intervals in the domain of the function are bounded */ if emptyp(tmp2) then( /*Get rid of logical operators*/ tmp:makelist(makelist(inpart(tmp[k],i),i,1,2),k,1,ll), /*Use canonical ordering*/ tmp:apply1(tmp,leq,geq,lss,grtr), /*Extract subintervals*/ subint:makelist(makelist(inpart(tmp[i][j],j),j,1,2),i,1,ll), /*Check for ill defined subintervals*/ tmp3:map(lambda([L],lfreeof(L,var)),subint), if not emptyp(sublist_indices(tmp3,lambda([x],x=false))) then error("Check function definition, some interval(s) appear to be empty or ill-defined. Read documentation for further details.") else (ans:makelist([subint[k],subval[k]],k,1,ll), return(ans) ) ) /*Function is not bounded*/ else( if not is(tmp2=[1,ll]) then error("Check function definition, some interval(s) appear to be empty or ill-defined. Read documentation for further details.") else (if is(ll=2) then( tmp:apply1(tmp,leq,geq,lss,grtr), leftU:[[minf,inpart(tmp[1],2)],subval[1]], rightU:[[inpart(tmp[2],1),inf],subval[2]], ans:[leftU,rightU], return(ans) ) elseif is(ll>2) then( leftU:apply1(tmp[1],leq,geq,lss,grtr), rightU:apply1(tmp[ll],leq,geq,lss,grtr), leftU:[[minf,inpart(leftU,2)],subval[1]], rightU:[[inpart(rightU,1),inf],subval[ll]], /*Get rid of logical operators in the bounded intervals*/ tmpb:makelist(makelist(inpart(tmp[k],i),i,1,2),k,2,ll-1), /*Use canonical ordering*/ tmpb:apply1(tmpb,leq,geq,lss,grtr), /*Extract bounded subintervals*/ subint:makelist(makelist(inpart(tmpb[i][j],j),j,1,2),i,1,ll-2), /*Check for ill-defined subintervals*/ tmp3:map(lambda([L],lfreeof(L,var)),subint), if not emptyp(sublist_indices(tmp3,lambda([x],x=false))) then error("Check function definition, some interval(s) appear to be empty or ill-defined. Read documentation for further details.") else ( ans:append([leftU],makelist([subint[k],subval[k+1]],k,1,ll-2),[rightU]), return(ans)) ) ) ) ))$ /*bint_comp(L1,L2,var) Bounded intervals comparison Input: Two lists, L1 and L2, each one having the following format [[a_i,b_i],expr_i] Output: A flag used by parityL */ bint_comp(L1,L2,var):=block( if is(L1[1][1]=-L2[1][2]) and is(L1[1][2]=-L2[1][1]) then( if is(L1[2]=0) and is(L2[2]=0) then return('zero) elseif equalp(L1[2],ratsubst(-var,var,L2[2])) then return('even) elseif equalp(L1[2],-ratsubst(-var,var,L2[2])) then return('odd) else return('none) ) else return('none) )$ /* uint_comp(L1,L2,var) Unbounded intervals comparison Input: Two list, L1 and L2, corresponding to unbounded intervals in the following format [[minf,b_1],expr1] or [[a_2,inf],expr2] Output: A flag used by parityL It is assumed that parity check has already checked if a_1,b_2 are equal to minf,inf respectively. */ uint_comp(L1,L2,var):=block( if is(L1[1][2]=-L2[1][1]) then( if is(L1[2]=0) and is(L2[2]=0) then return('zero) elseif equalp(L1[2],ratsubst(-var,var,L2[2])) then return('even) elseif equalp(L1[2],-ratsubst(-var,var,L2[2])) then return('odd) else return('none) ) else return('none) )$ /*parityb(L,var) Input: A list returned by fun_parts of a bounded function Output: A list used by parityL */ parityb(L,var):=block( [ll,icentral,aux,ans,llaux,Laux,paux], ll:length(L), /* Trivial case, just check if the only interval has an even expression, an odd expression or neither of them.*/ if is(ll=1) and is(L[1][1][1]=-L[1][1][2]) then( if evenfunp(L[1][2],var) then return('even) elseif oddfunp(L[1][2],var) then return('odd) else return('none) ) elseif is(ll>1) then( icentral:0, /*Search if there is a central interval*/ for i:1 thru ll do (if is(L[i][1][1]*L[i][1][2]<0) then icentral:i), /*Case 1: icentral=0 */ if is(icentral=0) then( if not evenp(ll) then return('none) else( aux:makelist(bint_comp(L[i],L[ll+1-i],var),i,1,ll/2), if is(length(sublist_indices(aux,lambda([x],x=even or x=zero)))=ll/2) then return('even) elseif is(length(sublist_indices(aux,lambda([x],x=odd or x=zero)))=ll/2) then return('odd) else return('none) ) ) /*icentral>0*/ elseif is(icentral>0) and oddp(ll) then( /*Check parity of central element*/ if is(L[icentral][1][1]=-L[icentral][1][2]) then ( if evenfunp(L[icentral][2],var) then paux:'even elseif oddfunp(L[icentral][2],var) then paux:'odd else return('none), Laux:delete(L[icentral],L), aux:makelist(bint_comp(Laux[i],Laux[ll-i],var),i,1,(ll-1)/2), if is(length(sublist_indices(aux,lambda([x],x=even or x=zero)))=(ll-1)/2) and is(paux=even) then return('even) elseif is(length(sublist_indices(aux,lambda([x],x=odd or x=zero)))=(ll-1)/2) and is(paux=odd) then return('odd) else return('none) ) else return('none) ) ))$ /*parityL(L,var) Input: A list returned by fun_parts Output: the parity of the function Important notes: - */ parityL(L,var):=block( [aux1,aux2,Laux,ll], ll:length(L), if (not boundedp(L)) and is(ll=2) then( if is(uint_comp(L[1],L[2],var)=zero) or is(uint_comp(L[1],L[2],var)=even) then return('even) elseif is(uint_comp(L[1],L[2],var)=odd) then return('odd) else return('none) ) /*Check if there are unbounded intervals*/ elseif (not boundedp(L)) and is(ll>2) then( aux1:uint_comp(L[1],L[ll],var), Laux:delete(L[1],L), Laux:delete(L[ll],Laux), aux2:parityb(Laux,var), if (is(aux1=zero) and is(aux2=even)) or (is(aux1=even) and is(aux2=even)) then return('even) elseif (is(aux1=zero) and is(aux2=odd)) or (is(aux1=odd) and is(aux2=odd)) then return('odd) else return('none) ) else return(parityb(L,var)) )$ boundedp(L):=block( [ll], ll:length(L), if is(L[1][1][1]=minf) and is (L[ll][1][2]=inf) then return(false) else return(true) )$ piecewisep(expr):=if atom(expr) then false else is(inpart(expr,0)="if")$ paritycheck(expr,var):=block( if piecewisep(expr) then return(parityL(fun_parts(expr,var),var)) elseif listp(expr) then return(parityL(expr,var)) else ( if equalp(ratsubst(-var,var,expr),expr) then return('even) elseif equalp(ratsubst(-var,var,expr),-expr) then return('odd) else return('none) ) )$ /*Check if the domain is bounded and symmetric [-L,L] symbintp(L) intput: a list L returned by fun_parts output: true or false */ symbintp(L):=if boundedp(L) and is(L[1][1][1]=-L[length(L)][1][2]) then L[length(L)][1][2] else false$ /*integratepw(L,var) Integrate the list form of a piecewise defined function */ integratepw(L,var):=block( [ll], ll:length(L), return(sum(integrate(L[i][2],var,L[i][1][1],L[i][1][2]),i,1,ll)) )$ /* sum2list(expr) Auxiliary function Input: An expression Output: A list whose elements are the terms of expression Important notes: sum2list does not check by itself if op(expr)="+" or not*/ sum2list(expr):= if is(expr=0) then [0] elseif is(nterms(expr)=1) then [expr] else args(expr)$ /*secsum2bl(L) Input: A list as in fun_parts output format Output: */ secsum2bl(L):=block( [ll,laux], ll:length(L), laux:makelist(i,i,1,ll), create_list([L[i][1],y],i,laux,y,sum2list(L[i][2])) )$ /*Match rules for special cases*/ matchdeclare(nexp,lambda([e],e#0 and nonnegintegerp(e)))$ matchdeclare(nfreq,lambda([e],e#0 and nonnegintegerp(e)))$ matchdeclare(const,lambda([e],e#0 and freeof(xargument,e)))$ /*matchdeclare(argtrig,lambda([e],e#0 and constantp(e)))$ defmatch(powx,const*xargument^nexp,xargument)$*/ defmatch(powxsin,const*xargument^nexp*sin(nfreq*%pi*xargument/pargument),xargument,pargument)$ defmatch(powxcos,const*xargument^nexp*cos(nfreq*%pi*xargument/pargument),xargument,pargument)$ defmatch(multsin,const*sin(nfreq*%pi*xargument/pargument),xargument,pargument)$ defmatch(multcos,const*cos(nfreq*%pi*xargument/pargument),xargument,pargument)$ /*powxcos_int Computes \int_a^b (x^r*cos(n*%pi*x/L)dx for r>=0 */ powxcos_int(pow,freq,p,a,b):=block( if is(freq=0) then (b^(pow+1)-a^(pow+1))/(pow+1) elseif is(pow=0) then (p/(%pi*freq))*(sin(freq*%pi*b/p)-sin(freq*%pi*a/p)) elseif is(pow>0) then (p/(%pi*freq))*(b^pow*sin(freq*%pi*b/p)-a^pow*sin(freq*%pi*a/p))-(p*pow/(freq*%pi)*powxsin_int(pow-1,freq,p,a,b)) )$ /*powxsin_int Computes \int_a^b x^r*sin(n*%pi*x/L) dx for r>=0 */ powxsin_int(pow,freq,p,a,b):=block( if is(freq=0) then 0 elseif is(pow=0) then (p/(%pi*freq))*(cos(freq*%pi*a/p)-cos(freq*%pi*b/p)) elseif is(pow>0) then (p/(%pi*freq))*(a^pow*cos(freq*%pi*a/p)-b^pow*cos(freq*%pi*b/p))+(p*pow/(freq*%pi)*powxcos_int(pow-1,freq,p,a,b)) )$ /*Heuristic a_n */ heuristic_an(expr,var,p,a,b):=block( [n,ans], declare(n,integer), if listp(powxsin(expr,var,p)) then ( if is(a=-b) and evenp(nexp) then 0 else( ans:(const/(2*p))*powxsin_int(nexp,nfreq+n,p,a,b)+ev((const/(2*p))*powxsin_int(nexp,nfreq-n,p,a,b),noeval,n), remove(n,integer), ans ) ) elseif listp(powxcos(expr,var,p)) then ( if is(a=-b) and oddp(nexp) then 0 else( ans:(const/(2*p))*powxcos_int(nexp,nfreq+n,p,a,b)+ev((const/(2*p))*powxcos_int(nexp,nfreq-n,p,a,b),noeval,n), remove(n,integer), ans) ) elseif listp(multsin(expr,var,p)) then( if is(a=-b) then 0 else( ans:(const/(2*p))*powxsin_int(0,nfreq+n,p,a,b)+ev((const/(2*p))*powxsin_int(0,nfreq-n,p,a,b),noeval,n), remove(n,integer), ans ) ) elseif listp(multcos(expr,var,p)) then ( ans:(const/(2*p))*powxcos_int(0,nfreq+n,p,a,b)+ev((const/(2*p))*powxcos_int(0,nfreq-n,p,a,b),noeval,n), remove(n,integer), ans ) else( ans:adefint(expr*cos(n*%pi*var/p),var,a,b)/p, remove(n,integer), ans ) )$ /*Heuristic b_n */ heuristic_bn(expr,var,p,a,b):=block( [n,ans], declare(n,integer), if listp(powxsin(expr,var,p)) then ( if is(a=-b) and oddp(nexp) then 0 else( ans:(-const/(2*p))*powxcos_int(nexp,nfreq+n,p,a,b)+ev((const/(2*p))*powxcos_int(nexp,nfreq-n,p,a,b),noeval,n), remove(n,integer), ans ) ) elseif listp(powxcos(expr,var,p)) then ( if is(a=-b) and evenp(nexp) then 0 else( ans:(const/(2*p))*powxsin_int(nexp,nfreq+n,p,a,b)+ev((const/(2*p))*powxsin_int(nexp,n-nfreq,p,a,b),noeval,n), remove(n,integer), ans ) ) elseif listp(multsin(expr,var,p)) then( ans:(-const/(2*p))*powxcos_int(0,nfreq+n,p,a,b)+ev((const/(2*p))*powxcos_int(0,nfreq-n,p,a,b),noeval,n), remove(n,integer), ans ) elseif listp(multcos(expr,var,p)) then ( if is(a=-b) then 0 else( ans:(const/(2*p))*powxsin_int(0,nfreq+n,p,a,b)+ev((const/(2*p))*powxsin_int(0,n-nfreq,p,a,b),noeval,n), remove(n,integer), ans ) ) else( ans:adefint(expr*sin(n*%pi*var/p),var,a,b)/p, remove(n,integer), ans ) )$ /*foucoeffpw(L,var) Fourier coefficients of a piecewise defined function Currently it only works with functions defined on an interval of the form [-p,p]. The output is a list of the form [a0,a1,b1,an,bn] */ foucoeffpw(L,var):=block( [llaux,ll,a0,an,bn,answ,lm], ll:length(L), lm:(L[ll][1][2]-L[1][1][1])/2, /*For more general intervals*/ Laux:secsum2bl(L), llaux:length(Laux), if is(parityL(L,var)=odd) then ( a0:0, an:0, bn:sum(heuristic_bn(Laux[i][2],var,lm,Laux[i][1][1],Laux[i][1][2]),i,1,llaux), answ:[a0,an,simplify_sum(bn)], return(ratsimp(answ))) elseif is(parityL(L,var)=even) then ( a0:(1/(2*lm))*integratepw(L,var), an:sum(heuristic_an(Laux[i][2],var,lm,Laux[i][1][1],Laux[i][1][2]),i,1,llaux), bn:0, answ:[simplify_sum(a0),simplify_sum(an),bn], return(ratsimp(answ))) else( a0:(1/(2*lm))*integratepw(L,var), an:sum(heuristic_an(Laux[i][2],var,lm,Laux[i][1][1],Laux[i][1][2]),i,1,llaux), bn:sum(heuristic_bn(Laux[i][2],var,lm,Laux[i][1][1],Laux[i][1][2]),i,1,llaux), answ:[a0,an,bn], answ:map(simplify_sum,answ), return(ratsimp(answ) )) )$ foucoeffterm(term,var,p):=block( [a0,bn,an,ans], if listp(term) then( /* if is(symbintp(term,var)=p) then */ return(foucoeffpw(term,var)) /* else error("Domain of the piecewise function is not valid. Read the documentation for further details")*/ ) elseif piecewisep(term) then( /*Check if the domain is valid*/ /*if is(symbintp(fun_parts(term,var))=p) then*/ if is(boundedp(fun_parts(term,var))) then return(foucoeffpw(map(lambda([e],[e[1],expand(trigrat(e[2]))]),fun_parts(term,var)),var)) else error("Domain of the piecewise function is not valid. Read the documentation for further details") ) else( /*Check parity of term*/ if is(paritycheck(term,var)=even) then( a0:ratsimp((1/(2*p))*adefint(term,var,-p,p)), an:heuristic_an(term,var,p,-p,p), bn:0, ans:ratsimp([a0,an,bn]), return(ans) ) elseif is(paritycheck(term,var)=odd) then( a0:0, an:0, bn:heuristic_bn(term,var,p,-p,p), ans:ratsimp([a0,an,bn]), return(ans) ) else( a0:ratsimp((1/(2*p))*adefint(term,var,-p,p)), an:heuristic_an(term,var,p,-p,p), bn:heuristic_bn(term,var,p,-p,p), ans:ratsimp([a0,an,bn]), return(ans) ) ) )$ list2bl(L):=block( [ll,laux], ll:length(L), laux:makelist(i,i,1,ll), create_list(y,i,laux,y,sum2list(L[i])) )$ trigpattern(term,var,p):=block( if listp(powxsin(term,var,p)) then return(nfreq) elseif listp(powxcos(term,var,p)) then return(nfreq) elseif listp(multcos(term,var,p)) then return(nfreq) elseif listp(multsin(term,var,p)) then return(nfreq) else return(false) )$ searchpoles_term(term,var,p):=block( [L,Laux,indx_aux,indx], if listp(term) then( L:secsum2bl(term), indx_aux:map(lambda([e],trigpattern(e[2],var,p)),L), indx:unique(sublist(indx_aux,integerp)), return(indx) ) elseif piecewisep(term) then( Laux:map(lambda([e],[e[1],expand(trigrat(e[2]))]),fun_parts(term,var)), L:secsum2bl(Laux), indx_aux:map(lambda([e],trigpattern(e[2],var,p)),L), indx:unique(sublist(indx_aux,integerp)), return(indx) ) else( indx_aux:trigpattern(term,var,p), if integerp(indx_aux) then return([indx_aux]) else return([]) ) )$ searchpoles(expr,var,p):=block( [L,Laux,indx], Laux:sum2list(expand(expr)), Laux:map(lambda([e],if not piecewisep(e) then expand(trigrat(e)) else e),Laux), L:list2bl(Laux), indx:apply(append,map(lambda([e],searchpoles_term(e,var,p)),L)), return(unique(indx)) )$ fouriercoeff(expr,var,p):=block( [ans,Llist,Laux,coeffpoles,poles,n], Laux:sum2list(expand(expr)), Llist:map(lambda([e],if not piecewisep(e) then expand(trigrat(e)) else e),Laux), Llist:list2bl(Llist), ans:lsum(y,y,map(lambda([e],foucoeffterm(e,var,p)),Llist)), ans:ratsimp(ans), poles:searchpoles(expr,var,p), coeffpoles:makelist(at([i,ans[2],ans[3]],n=i),i,poles), declare(n,integer), coeffpoles:ev(coeffpoles), ans:ev(ans), remove(n,integer), return(ratsimp([ans,coeffpoles])) )$ fouriercoeff_expand(fcoeff,var,p,NN):=block( [ans,a0,an,bn,poles,indx_poles,polessum,sl_indx,laux], [[a0,an,bn],poles]:fcoeff, an:ev(an), bn:ev(bn), if emptyp(poles) then( ans:a0+sum(at(an*cos(n*%pi*var/p)+bn*sin(n*%pi*var/p),n=n),n,1,NN), return(ans) ) else( indx_poles:makelist(poles[i][1],i,1,length(poles)), if is(NN=inf) then( polessum:sum(poles[i][2]*cos(poles[i][1]*%pi*var/p)+poles[i][3]*sin(poles[i][1]*%pi*var/p),i,1,length(poles)), ans:a0+polessum+sum(at(an*cos(n*%pi*var/p)+bn*sin(n*%pi*var/p),n=n),n,1,inf), if is(an=0) and is(bn=0) then return(ans), print("The sum is over \\N -",setify(indx_poles)), return(ans) ) else( sl_indx:sublist_indices(indx_poles,lambda([e],is(e<=NN))), polessum:lsum(poles[i][2]*cos(poles[i][1]*%pi*var/p)+poles[i][3]*sin(poles[i][1]*%pi*var/p),i,sl_indx), laux:listify(setdifference(setify(makelist(i,i,1,NN)),setify(indx_poles))), ans:a0+polessum+lsum(at(an*cos(n*%pi*var/p)+bn*sin(n*%pi*var/p),n=i),i,laux), return(ans) ) ) )$ fourier_series(expr,var,p,NN):=block( [fcoeff,ans], fcoeff:fouriercoeff(expr,var,p), ans:fouriercoeff_expand(fcoeff,var,p,NN), ans )$ fourier_amplitudes(expr,var,p,N):=block( [a0,an,bn,ans,poles,indx_poles,polelist,sl_indx,laux], [[a0,an,bn],poles]:fouriercoeff(expr,var,p), an:ev(an), bn:ev(bn), if emptyp(poles) then( ans:makelist(at([n,sqrt(an^2+bn^2)],n=i),i,1,N), /*ans:float(ans),*/ return(ans)) else( indx_poles:makelist(poles[i][1],i,1,length(poles)), sl_indx:sublist_indices(indx_poles,lambda([e],is(e<N))), polelist:makelist([poles[i][1],sqrt(poles[i][2]^2+poles[i][3]^2)],i,sl_indx), laux:listify(setdifference(setify(makelist(i,i,1,N)),setify(indx_poles))), ans:append(polelist,makelist(at([n,sqrt(an^2+bn^2)],n=i),i,laux)), /*ans:float(ans),*/ ans:sort(ans,lambda([a,b],a[1]<b[1])), return(ans) ) )$ fourier_freq_list(expr,var,p,N):=block( [modules], modules:float(fourier_amplitudes(expr,var,p,N)) )$ fourier_freq(expr,var,p,N):=block( [modules], modules:float(fourier_amplitudes(expr,var,p,N)), draw2d(points_joined=impulses,line_width=4,color="blue",points(modules), xlabel="w(n)=nw_0",ylabel="|c_n|",axis_top=false,axis_right=false,xtics=1,user_preamble="set grid ytics") )$ wxfourier_freq(expr,var,p,N):=block( [modules], modules:float(fourier_amplitudes(expr,var,p,N)), wxdraw2d(points_joined=impulses,line_width=4,color="blue",points(modules), xlabel="w(n)=nw_0",ylabel="|c_n|",axis_top=false,axis_right=false,xtics=1,user_preamble="set grid ytics") )$ atan_fourier(an,bn):=if is(an=0) and is(bn=0) then 0 else atan2(an,bn)$ fourier_phshift(expr,var,p,N):=block( [a0,an,bn,ans,poles,indx_poles,polelist,sl_indx,laux], [[a0,an,bn],poles]:fouriercoeff(expr,var,p), an:ev(an), bn:ev(bn), if emptyp(poles) then( ans:makelist([i,atan_fourier(at(bn,n=i),at(an,n=i))],i,1,N), /*ans:float(ans),*/ return(ans)) else( indx_poles:makelist(poles[i][1],i,1,length(poles)), sl_indx:sublist_indices(indx_poles,lambda([e],is(e<N))), polelist:makelist([poles[i][1],atan_fourier(poles[i][3],poles[i][2])],i,sl_indx), laux:listify(setdifference(setify(makelist(i,i,1,N)),setify(indx_poles))), ans:append(polelist,makelist([i,atan_fourier(at(bn,n=i),at(an,n=i))],i,laux)), /*ans:float(ans),*/ ans:sort(ans,lambda([a,b],a[1]<b[1])), return(ans) ) )$ fourier_harm(expr,var,p,N):=block( [amplitudes,phshift,ans], amplitudes:fourier_amplitudes(expr,var,p,N), phshift:fourier_phshift(expr,var,p,N), ans:makelist(amplitudes[i][2]*cos(i*%pi*var/p-phshift[i][2]),i,1,N), return(ans) )$ chop(expr,[N]):=if emptyp(N) then scanmap(lambda([x],if is(numberp(x)) then (if is(abs(x)<1.0*10^(-12)) then 0.0 else x) else x),expr) else scanmap(lambda([x],if is(numberp(x)) then (if is(abs(x)<1.0*10^(-N[1])) then 0.0 else x) else x),expr)$ oddextension_pwterm(L,var):=block( [Laux], Laux:[[-L[1][2],-L[1][1]],-ratsubst(-x,x,L[2])], Laux )$ oddextensionpw(L,var):=block( [Lodd], Lodd:append(reverse(map(lambda([e],oddextension_pwterm(e,var)),L)),L), Lodd )$ evenextension_pwterm(L,var):=block( [Laux], Laux:[[-L[1][2],-L[1][1]],ratsubst(-x,x,L[2])], Laux )$ evenextensionpw(L,var):=block( [Lodd], Lodd:append(reverse(map(lambda([e],evenextension_pwterm(e,var)),L)),L), Lodd )$ odd_extension(expr,var,p):=block( [Laux], if piecewisep(expr) then( Laux:fun_parts(expr,var), if is(Laux[1][1][1]=0) and is(Laux[length(Laux)][1][2]=p) then return(oddextensionpw(Laux,var)) else error("The domain must be of the form [0,p]") ) else( Laux:[[[0,p],expr]], return(oddextensionpw(Laux,var)) ) )$ even_extension(expr,var,p):=block( [Laux], if piecewisep(expr) then( Laux:fun_parts(expr,var), if is(Laux[1][1][1]=0) and is(Laux[length(Laux)][1][2]=p) then return(evenextensionpw(Laux,var)) else error("The domain must be of the form [0,p]") ) else( Laux:[[[0,p],expr]], return(evenextensionpw(Laux,var)) ) )$ /*foucoscoeff(expr,var,p):=block( [Laux,fcoeff,a0,an,bn,poles], Laux:even_extension(expr,var,p), [[a0,an,bn],poles]:fouriercoeff(Laux,var,p), fcoeff:[[a0,an],map(lambda([e],[e[1],e[2]]),poles)], return(fcoeff) )$ fousincoeff(expr,var,p):=block( [Laux,fcoeff,a0,an,bn,poles], Laux:odd_extension(expr,var,p), [[a0,an,bn],poles]:fouriercoeff(Laux,var,p), fcoeff:[[bn],map(lambda([e],[e[1],e[3]]),poles)], return(fcoeff) )$ */ /*Here we have the new implementation of fouriersincoeff and fouriercoscoeff*/ fouriersincoeffpw(Llist,var):=block( [llaux,ll,bn,answ,lm], ll:length(Llist), lm:(Llist[ll][1][2]-Llist[1][1][1]), /*No need to divide by 2*/ Laux:secsum2bl(Llist), llaux:length(Laux), bn:2*sum(heuristic_bn(Laux[i][2],var,lm,Laux[i][1][1],Laux[i][1][2]),i,1,llaux), answ:[bn], answ:map(simplify_sum,answ), return(ratsimp(answ) ) )$ fouriercoscoeffpw(Llist,var):=block( [llaux,ll,a0,an,answ,lm], ll:length(Llist), lm:(Llist[ll][1][2]-Llist[1][1][1]), /*No need to divide by 2*/ Laux:secsum2bl(Llist), llaux:length(Laux), a0:(1/(lm))*integratepw(Llist,var), an:2*sum(heuristic_an(Laux[i][2],var,lm,Laux[i][1][1],Laux[i][1][2]),i,1,llaux), answ:[a0,an], answ:map(simplify_sum,answ), return(ratsimp(answ) ) )$ fouriersincoeffterm(term,var,p):=block( [bn,ans], if listp(term) then( /* if is(symbintp(term,var)=p) then */ return(fouriersincoeffpw(term,var)) /* else error("Domain of the piecewise function is not valid. Read the documentation for further details")*/ ) elseif piecewisep(term) then( /*Check if the domain is valid*/ /*if is(symbintp(fun_parts(term,var))=p) then*/ if is(boundedp(fun_parts(term,var))) then return(fouriersincoeffpw(map(lambda([e],[e[1],expand(trigrat(e[2]))]),fun_parts(term,var)),var)) else error("Domain of the piecewise function is not valid. Read the documentation for further details") ) else( bn:2*heuristic_bn(term,var,p,0,p), ans:ratsimp([bn]), return(ans) ) )$ fouriercoscoeffterm(term,var,p):=block( [a0,an,ans], if listp(term) then( /* if is(symbintp(term,var)=p) then */ return(fouriercoscoeffpw(term,var)) /* else error("Domain of the piecewise function is not valid. Read the documentation for further details")*/ ) elseif piecewisep(term) then( /*Check if the domain is valid*/ /*if is(symbintp(fun_parts(term,var))=p) then*/ if is(boundedp(fun_parts(term,var))) then return(fouriercoscoeffpw(map(lambda([e],[e[1],expand(trigrat(e[2]))]),fun_parts(term,var)),var)) else error("Domain of the piecewise function is not valid. Read the documentation for further details") ) else( a0:ratsimp((1/p)*adefint(term,var,0,p)), an:2*heuristic_an(term,var,p,0,p), ans:ratsimp([a0,an]), return(ans) ) )$ fouriersincoeff(expr,var,p):=block( [ans,Llist,Laux,coeffpoles,poles,n], Laux:sum2list(expand(expr)), Llist:map(lambda([e],if not piecewisep(e) then expand(trigrat(e)) else e),Laux), Llist:list2bl(Llist), ans:lsum(y,y,map(lambda([e],fouriersincoeffterm(e,var,p)),Llist)), ans:ratsimp(ans), poles:searchpoles(expr,var,p), /*Now there is only a list of the form [bn]*/ coeffpoles:makelist(at([i,ans[1]],n=i),i,poles), declare(n,integer), coeffpoles:ev(coeffpoles), ans:ev(ans), remove(n,integer), return(ratsimp([ans,coeffpoles])) )$ fouriercoscoeff(expr,var,p):=block( [ans,Llist,Laux,coeffpoles,poles,n], Laux:sum2list(expand(expr)), Llist:map(lambda([e],if not piecewisep(e) then expand(trigrat(e)) else e),Laux), Llist:list2bl(Llist), ans:lsum(y,y,map(lambda([e],fouriercoscoeffterm(e,var,p)),Llist)), ans:ratsimp(ans), poles:searchpoles(expr,var,p), /*Now there is only a list of the form [a0,an]*/ coeffpoles:makelist(at([i,ans[2]],n=i),i,poles), declare(n,integer), coeffpoles:ev(coeffpoles), ans:ev(ans), remove(n,integer), return(ratsimp([ans,coeffpoles])) )$ fouriersincoeff_expand(sincoeff,var,p,NN):=block( [ans,bn,poles,indx_poles,polessum,sl_indx,laux], [[bn],poles]:sincoeff, bn:ev(bn), if emptyp(poles) then( ans:sum(at(bn*sin(n*%pi*var/p),n=n),n,1,NN), return(ans) ) else( indx_poles:makelist(poles[i][1],i,1,length(poles)), if is(NN=inf) then( polessum:sum(poles[i][2]*sin(poles[i][1]*%pi*var/p),i,1,length(poles)), ans:polessum+sum(at(bn*sin(n*%pi*var/p),n=n),n,1,inf), if is(bn=0) then return(ans), print("The sum is over \\N-",setify(indx_poles)), return(ans) ) else( sl_indx:sublist_indices(indx_poles,lambda([e],is(e<=NN))), polessum:lsum(poles[i][2]*sin(poles[i][1]*%pi*var/p),i,sl_indx), laux:listify(setdifference(setify(makelist(i,i,1,NN)),setify(indx_poles))), ans:polessum+lsum(at(bn*sin(n*%pi*var/p),n=i),i,laux), return(ans) ) ) )$ fouriercoscoeff_expand(coscoeff,var,p,NN):=block( [ans,a0,an,poles,indx_poles,polessum,sl_indx,laux], [[a0,an],poles]:coscoeff, an:ev(an), if emptyp(poles) then( ans:a0+sum(at(an*cos(n*%pi*var/p),n=n),n,1,NN), return(ans) ) else( indx_poles:makelist(poles[i][1],i,1,length(poles)), if is(NN=inf) then( polessum:sum(poles[i][2]*cos(poles[i][1]*%pi*var/p),i,1,length(poles)), ans:a0+polessum+sum(at(an*cos(n*%pi*var/p),n=n),n,1,inf), if is(an=0) then return(ans), print("The sum is over \\N -",setify(indx_poles)), return(ans) ) else( sl_indx:sublist_indices(indx_poles,lambda([e],is(e<=NN))), polessum:lsum(poles[i][2]*cos(poles[i][1]*%pi*var/p),i,sl_indx), laux:listify(setdifference(setify(makelist(i,i,1,NN)),setify(indx_poles))), ans:a0+polessum+lsum(at(an*cos(n*%pi*var/p),n=i),i,laux), return(ans) ) ) )$ fouriersin_series(expr,var,p,NN):=block( [ans,sincoeff], sincoeff:fouriersincoeff(expr,var,p), ans:fouriersincoeff_expand(sincoeff,var,p,NN), ans )$ fouriercos_series(expr,var,p,NN):=block( [ans,coscoeff], coscoeff:fouriercoscoeff(expr,var,p), ans:fouriercoscoeff_expand(coscoeff,var,p,NN), ans )$ real2complex_fcoeff(fcoeff):=block( [c0,cn,cpoles,rpoles,a0,an,bn], [[a0,an,bn],rpoles]:fcoeff, c0:a0, cn:(an-%i*bn)/2, cpoles:map(lambda([e],[e[1],(e[2]-%i*e[3])/2]),rpoles), ratsimp([[c0,cn],cpoles]) )$ cfouriercoeff(expr,var,p):=real2complex_fcoeff(fouriercoeff(expr,var,p))$ cfouriercoeff_expand(cfcoeff,var,p,NN):=block( [ans,c0,cn,cpoles,indx_cpoles,cpolessumpos,cpolessumneg,sl_indx,laux], [[c0,cn],cpoles]:cfcoeff, an:ev(cn), if emptyp(cpoles) then( if is(NN=inf) then( ans:c0+sum(at(cn*exp(%i*n*%pi*var/p),n=n),n,minf,inf), print("The sum is over \\N-\{0\}"), return(ans) ) else( ans:c0+sum(at(conjugate(cn)*exp(%i*-n*%pi*var/p),n=n),n,1,NN) +sum(at(cn*exp(%i*n*%pi*var/p),n=n),n,1,NN), return(ratsimp(demoivre(ans)))) ) else( indx_cpoles:map(first,cpoles), if is(NN=inf) then( cpolessumpos:sum(cpoles[i][2]*exp(%i*cpoles[i][1]*%pi*var/p),i,1,length(cpoles)), cpolessumneg:sum(conjugate(cpoles[i][2])*exp(-%i*cpoles[i][1]*%pi*var/p),i,1,length(cpoles)), ans:c0+cpolessumpos+cpolessumneg+sum(at(cn*exp(%i*n*%pi*var/p),n=n),n,minf,inf), if is(cn=0) then return(ans), print("The sum is over \\N -",setify(append(indx_cpoles,[0]))), return(ans) ) else( sl_indx:sublist_indices(indx_cpoles,lambda([e],is(e<=NN))), cpolessumpos:lsum(cpoles[i][2]*exp(%i*cpoles[i][1]*%pi*var/p),i,sl_indx), cpolessumneg:lsum(conjugate(cpoles[i][2])*exp(-%i*cpoles[i][1]*%pi*var/p),i,sl_indx), laux:listify(setdifference(setify(makelist(i,i,1,NN)),setify(indx_cpoles))), ans:c0+cpolessumpos+cpolessumneg+lsum(at(cn*exp(%i*n*%pi*var/p),n=i),i,laux) +lsum(at(conjugate(cn)*exp(-%i*n*%pi*var/p),n=i),i,laux), return(ratsimp(demoivre(ans))) ) ) )$ cfourier_series(expr,var,p,NN):=cfouriercoeff_expand(cfouriercoeff(expr,var,p),var,p,NN)$ evalpole_scn(bfn,polesf,indx_pf,indx):=block([n], if member(indx,indx_pf) then return(assoc(indx,polesf)) else return(at(bfn,n=indx)) )$ fourier_heatcoeff(Qexpr,fexpr,var,p,kap):=block( [bnf,bnQ,polesf,polesQ,indxpf,indxpQ,indxp,n,un,polesU,ans], /* [[bnf],polesf]:fousincoeff(fexpr,var,p), [[bnQ],polesQ]:fousincoeff(Qexpr,var,p), */ [[bnf],polesf]:fouriersincoeff(fexpr,var,p), [[bnQ],polesQ]:fouriersincoeff(Qexpr,var,p), declare(n,integer), un:bnf*exp(-(n*%pi/p)^2*kap*t)+integrate(at(bnQ,t=s)*exp(-(n*%pi/p)^2*kap*(t-s)),s,0,t), if emptyp(polesf) and emptyp(polesQ) then (ans:[[un],[]], remove(n,integer), return(ratsimp(ans))) /*now let's deal with the poles*/ else( indxpf:map(first,polesf), indxpQ:map(first,polesQ), indxp:unique(append(indxpf,indxpQ)), polesU:map(lambda([e],[e,evalpole_scn(bnf,polesf,indxpf,e)*exp(-(e*%pi/p)^2*kap*t) +integrate(at(evalpole_scn(bnQ,polesQ,indxpQ,e),t=s)*exp(-(e*%pi/p)^2*kap*(t-s)),s,0,t)]),indxp), ans:[[un],polesU], remove(n,integer), return(ratsimp(ans)) ) )$ fourier_heat(Qexpr,fexpr,var,p,kap,NN):=fouriersincoeff_expand(fourier_heatcoeff(Qexpr,fexpr,var,p,kap),var,p,NN)$ fourier_wavecoeff(Qexpr,fexpr,gexpr,var,p,A):=block( [indxp,indxpf,indxpQ,indxpg,n,polesf,polesg,polesQ,bnf,bnQ,bng,Ln,un,polesU], Ln:(n*%pi*A)/p, /* [[bnf],polesf]:fousincoeff(fexpr,var,p), [[bng],polesg]:fousincoeff(gexpr,var,p), [[bnQ],polesQ]:fousincoeff(Qexpr,var,p), */ [[bnf],polesf]:fouriersincoeff(fexpr,var,p), [[bng],polesg]:fouriersincoeff(gexpr,var,p), [[bnQ],polesQ]:fouriersincoeff(Qexpr,var,p), bng:bng/Ln, polesg:map(lambda([e],[e[1],e[2]*at(1/Ln,n=e[1])]),polesg), declare(n,integer), un:bnf*cos(Ln*t)+bng*sin(Ln*t)+(1/Ln)*integrate(at(bnQ,t=s)*sin(Ln*(t-s)),s,0,t), if emptyp(polesf) and emptyp(polesQ) and emptyp(polesg) then (ans:[[un],[]], remove(n,integer), return(ratsimp(ans))) /*now let's deal with the poles*/ else( indxpf:map(first,polesf), indxpQ:map(first,polesQ), indxpg:map(first,polesg), indxp:unique(append(indxpf,indxpQ,indxpg)), polesU:map(lambda([e],[e,evalpole_scn(bnf,polesf,indxpf,e)*cos(e*%pi*A*t/p)+ evalpole_scn(bng,polesg,indxpg,e)*sin(e*%pi*A*t/p)+ (p/(e*%pi*A))*integrate(at(evalpole_scn(bnQ,polesQ,indxpQ,e),t=s)*sin(e*%pi*A*(t-s)/p),s,0,t)]),indxp), ans:[[un],polesU], remove(n,integer), return(ratsimp(ans)) ) )$ fourier_wave(Qexpr,fexpr,gexpr,var,p,A,NN):=fouriersincoeff_expand(fourier_wavecoeff(Qexpr,fexpr,gexpr,var,p,A),var,p,NN)$ dirichlet_heat(Qexpr,fexpr,h1expr,h2expr,var,p,kap,NN):=block( [Qaux,faux,vaux], vaux:(h1expr-h2expr)*var/p-h1expr, Qaux:Qexpr+diff(vaux,t), faux:fexpr+at(vaux,t=0), fourier_heat(Qaux,faux,var,p,kap,NN)-vaux )$ dirichlet_wave(Qexpr,fexpr,gexpr,h1expr,h2expr,var,p,A,NN):=block( [Qaux,faux,gaux,vaux], vaux:(h1expr-h2expr)*var/p-h1expr, Qaux:Qexpr+diff(vaux,t,2), faux:fexpr+at(vaux,t=0), gaux:gexpr+at(diff(vaux,t),t=0), fourier_wave(Qaux,faux,gaux,var,p,A,NN)-vaux )$ dirichlet_laplace_disk(a,fexpr,var,NN):=block( [cn,c0,cpoles,n], [[c0,cn],cpoles]:cfouriercoeff(fexpr,var,%pi), cn:(r/a)^abs(n)*cn, cpoles:map(lambda([e],[e[1],(r/a)^e[1]*e[2]]),cpoles), cfouriercoeff_expand([[c0,cn],cpoles],var,%pi,NN) )$ neumann_laplace_disk(a,fexpr,var,NN):=block([a0,an,bn,poles,n], [[a0,an,bn],poles]:fouriercoeff(fexpr,var,%pi), if not is(a0=0) then return("The solution does not exist.") else( an:(a/n)*an*(r/a)^n, bn:(a/n)*bn*(r/a)^n, poles:map(lambda([e],[e[1],e[2]*(a/e[1])*(r/a)^e[1],e[3]*(a/e[1])*(r/a)^e[1]]),poles), return(fouriercoeff_expand([[a0,an,bn],poles],var,%pi,NN)) ) )$ dirichlet_laplace_annulus(a,b,fexpr,gexpr,var,NN):=block( [n,cnf,c0f,cng,c0g,cpolesf,cpolesg,c0aux,cnaux,cpolesaux, indxpf,indxpg,indxp,ans], [[c0f,cnf],cpolesf]:cfouriercoeff(fexpr,var,%pi), [[c0g,cng],cpolesg]:cfouriercoeff(gexpr,var,%pi), c0aux:(c0g*log(b)-c0f*log(a))*log(b/a)+(c0f-c0g)*log(r)/log(b/a), cnaux:(b^abs(n)*r^abs(n)-a^abs(2*n)*b^abs(n)*r^(-abs(n)))*cnf/(b^abs(2*n)-a^abs(2*n))+ cng*(b^abs(2*n)*a^abs(n)*r^(-abs(n))-a^abs(n)*r^abs(n))/(b^abs(2*n)-a^abs(2*n)), if emptyp(cpolesf) and emptyp(cpolesg) then( ans:cfouriercoeff_expand(ratsimp([[c0aux,cnaux],[]]),var,%pi,NN), return(ans) ) /*now let's deal with the poles*/ else( indxpf:map(first,cpolesf), indxpg:map(first,cpolesg), indxp:unique(append(indxpf,indxpg)), cpolesaux:map(lambda([e],[e,(b^e*r^e-a^(2*e)*b^e*r^(-e))*evalpole_scn(cnf,cpolesf,indxpf,e)/(b^(2*e)-a^(2*e))+ evalpole_scn(cng,cpolesg,indxpg,e)*(b^(2*e)*a^e*r^(-e)-a^e*r^e)/(b^(2*e)-a^(2*e))]),indxp), ans:cfouriercoeff_expand(ratsimp([[c0aux,cnaux],cpolesaux]),var,%pi,NN), return(ans) ) )$ dirichlet_laplace_wedge(R,a,fexpr,var,NN):=block( [bnf,polesf,bnaux,polesaux,n,p:a], [[bnf],polesf]:fouriersincoeff(fexpr,var,p), bnaux:bnf*(r/R)^(n*%pi/a), if emptyp(polesf) then return(fouriersincoeff_expand([[bnaux],[]],var,p,NN)) else( polesaux:map(lambda([e],[e[1],e[2]*(r/R)^(e[1]*%pi/a)]),polesf), return(fouriersincoeff_expand([[bnaux],polesaux],var,p,NN)) ) )$ neumann_laplace_wedge(R,a,fexpr,var,NN):=block( [bnf,polesf,bnaux,polesaux,n,p:a], [[bnf],polesf]:fouriersincoeff(fexpr,var,p), bnaux:bnf*(r/R)^(n*%pi/a)*(R*a/(n*%pi)), if emptyp(polesf) then return(fouriersincoeff_expand([[bnaux],[]],var,p,NN)) else( polesaux:map(lambda([e],[e[1],e[2]*(r/R)^(e[1]*%pi/a)*(R*a/(e[1]*%pi))]),polesf), return(fouriersincoeff_expand([[bnaux],polesaux],var,p,NN)) ) )$ dirichlet_laplace_rectangle(a,b,f1expr,f2expr,g1expr,g2expr,var1,var2,NN):=block( [n,f1n,polesf1,f2n,polesf2,g1n,polesg1,g2n,polesg2,indxpf1,indxpf2,indxpg1,indxpg2, indxp,Fnaux,Gnaux,polesFaux,polesGaux,ans], [[f1n],polesf1]:fouriersincoeff(f1expr,var1,a), [[f2n],polesf2]:fouriersincoeff(f2expr,var1,a), [[g1n],polesg1]:fouriersincoeff(g1expr,var2,b), [[g2n],polesg2]:fouriersincoeff(g2expr,var2,b), Fnaux:f1n*sinh(n*%pi*(b-var2)/a)/sinh(n*%pi*b/a)+f2n*sinh(n*%pi*var2/a)/sinh(n*%pi*b/a), Gnaux:g1n*sinh(n*%pi*(a-var1)/b)/sinh(n*%pi*a/b)+g2n*sinh(n*%pi*var1/b)/sinh(n*%pi*a/b), if emptyp(polesf1) and emptyp(polesf2) and emptyp(polesg1) and emptyp(polesg2) then( ans:fouriersincoeff_expand([[Fnaux],[]],var1,a,NN)+fouriersincoeff_expand([[Gnaux],[]],var2,b,NN), return(ans) ) else( indxpf1:map(first,polesf1), indxpf2:map(first,polesf2), indxpg1:map(first,polesg1), indxpg2:map(first,polesg2), indxp:unique(append(indxpf1,indxpf2,indxpg1,indxpg2)), polesFaux:map(lambda([e],[e,evalpole_scn(f1n,polesf1,indxpf1,e)*sinh(e*%pi*(b-var2)/a)/sinh(e*%pi*b/a) +evalpole_scn(f2n,polesf2,indxpf2,e)*sinh(e*%pi*var2/a)/sinh(e*%pi*b/a)]),indxp), polesGaux:map(lambda([e],[e,evalpole_scn(g1n,polesg1,indxpg1,e)*sinh(e*%pi*(a-var1)/b)/sinh(e*%pi*a/b) +evalpole_scn(g2n,polesg2,indxpg2,e)*sinh(e*%pi*var1/b)/sinh(e*%pi*a/b)]),indxp), ans:fouriersincoeff_expand([[Fnaux],polesFaux],var1,a,NN)+fouriersincoeff_expand([[Gnaux],polesGaux],var2,b,NN), return(ans) ) )$ neumann_laplace_rectangle(a,b,f1expr,f2expr,g1expr,g2expr,var1,var2,NN):=block( [n,f1n,polesf1,f2n,polesf2,g1n,polesg1,g2n,polesg2,indxpf1,indxpf2,indxpg1,indxpg2, indxp,Fnaux,Gnaux,polesFaux,polesGaux,C,ans,f10,g10,f20,g20], [[f10,f1n],polesf1]:fouriercoscoeff(f1expr,var1,a), [[f20,f2n],polesf2]:fouriercoscoeff(f2expr,var1,a), [[g10,g1n],polesg1]:fouriercoscoeff(g1expr,var2,b), [[g20,g2n],polesg2]:fouriercoscoeff(g2expr,var2,b), Fnaux:(-a/(n*%pi))*f1n*cosh(n*%pi*(b-var2)/a)/sinh(n*%pi*b/a)+(a/(n*%pi))*f2n*cosh(n*%pi*var2/a)/sinh(n*%pi*b/a), Gnaux:(-b/(n*%pi))*g1n*cosh(n*%pi*(a-var1)/b)/sinh(n*%pi*a/b)+(b/(n*%pi))*g2n*cosh(n*%pi*var1/b)/sinh(n*%pi*a/b), if emptyp(polesf1) and emptyp(polesf2) and emptyp(polesg1) and emptyp(polesg2) then( ans:fouriercoscoeff_expand([[0,Fnaux],[]],var1,a,NN)+fouriercoscoeff_expand([[0,Gnaux],[]],var2,b,NN), return(ans+C) ) else( indxpf1:map(first,polesf1), indxpf2:map(first,polesf2), indxpg1:map(first,polesg1), indxpg2:map(first,polesg2), indxp:unique(append(indxpf1,indxpf2,indxpg1,indxpg2)), polesFaux:map(lambda([e],[e,(-a/(n*%pi))*evalpole_scn(f1n,polesf1,indxpf1,e)*cosh(e*%pi*(b-var2)/a)/sinh(e*%pi*b/a) +(a/(n*%pi))*evalpole_scn(f2n,polesf2,indxpf2,e)*cosh(e*%pi*var2/a)/sinh(e*%pi*b/a)]),indxp), polesGaux:map(lambda([e],[e,(-b/(n*%pi))*evalpole_scn(g1n,polesg1,indxpg1,e)*cosh(e*%pi*(a-var1)/b)/sinh(e*%pi*a/b) +(b/(n*%pi))*evalpole_scn(g2n,polesg2,indxpg2,e)*cosh(e*%pi*var1/b)/sinh(e*%pi*a/b)]),indxp), ans:fouriercoscoeff_expand([[0,Fnaux],polesFaux],var1,a,NN)+fouriercoscoeff_expand([[0,Gnaux],polesGaux],var2,b,NN), return(ans+C) ) )$ /*Beggining of 2D Wave equation*/ BesselJZeros(m,k,[u]):=block([a,alpha,A,answ],local(a,alpha,A), if is(not(m>-1)) then return("The order m must be m>-1"), if (is(not(integerp(k))) or is(not(k > 0))) then return("k must be an integer k>0"), a[i]:=m+2*i, alpha[i,j]:=if is(equal(i,j)) then 2/((a[i]-1)*(a[i]+1)) elseif is(equal(j,i-1)) then 1/((a[i]-1)*sqrt((a[i]-2)*a[i])) elseif is(equal(j,i+1)) then 1/((a[i+1]-1)*sqrt((a[i+1]-2)*a[i+1])) else 0, A:genmatrix(alpha,ceiling(k/0.64)+10,ceiling(k/0.64)+10), answ:2/sqrt(first(eigens_by_jacobi(A,floatfield))), if (is(not(equal(length(u),0))) and is(equal(first(u),all))) then firstn(sort(answ),k) elseif (is(not(equal(length(u),0))) and is(not(equal(first(u),all)))) then return("Do you want the first k roots? In this case, the option is 'all'") else last(firstn(sort(answ),k)) )$ simpgen(l1,l2):=append(endcons(last(l1)+first(l2),rest(l1,-1)),rest(l2,1))$ pairing(L1,L2):=xreduce(append,map(lambda([x],map(lambda([y],[x,y]),L2)),L1))$ nintegrate(f,m,[I]):=block([n:length(args(lhs(apply(fundef,[f])))),tmp:xreduce(simpgen,makelist([1,4,1],j,1,m/2)),x,deltax,p,w,ww,xx],local(x,deltax,p,w), if is(not(evenp(m))) then return ("The number of subdivisions m must be even") elseif is(not(equal(length(I),n))) then return("The number of variables xi and intervals [ai,bi] must match") else ( for i:1 thru n do x[i]:args(lhs(apply(fundef,[f])))[i], for i:1 thru n do deltax[i]:(I[i][2]-I[i][1])/m, for j:1 thru n do p[j]:makelist(I[j][1]+k*deltax[j],k,0,m), for j:1 thru n do w[j]:(deltax[j]/3)*tmp, if is(equal(n,1)) then ww:w[1] else ww:map(lambda([x],apply("*",x)),map(flatten,xreduce(pairing,makelist(w[j],j,1,n)))), if is(equal(n,1)) then xx:p[1] else xx:map(flatten,xreduce(pairing,makelist(p[j],j,1,n))), if is(equal(n,1)) then bfloat(apply("+",ww*map(f,xx))) else bfloat(apply("+",ww*map(lambda([x],apply(f,x)),xx))) ) )$ a0(n,a,f):=block([r:args(lhs(apply(fundef,[f])))[1],theta:args(lhs(apply(fundef,[f])))[2],lam:BesselJZeros(0,n)/a,p:%pi,ff,k],local(ff), define(funmake(ff,[r,theta]),apply(f,[r,theta])*r*bessel_j(0,lam*r)), k:1/(p*a^2*(bessel_j(1,lam))^2), k*nintegrate(ff,14,[0,a],[0,2*p]) )$ a(m,n,a,f):=block([r:args(lhs(apply(fundef,[f])))[1],theta:args(lhs(apply(fundef,[f])))[2],lam:BesselJZeros(0,n)/a,p:%pi,ff,k],local(ff), define(funmake(ff,[r,theta]),apply(f,[r,theta])*r*bessel_j(0,lam*r)*cos(m*theta)), k:2/(p*a^2*(bessel_j(m+1,lam))^2), k*nintegrate(ff,14,[0,a],[0,2*p]) )$ b(m,n,a,f):=block( [r:args(lhs(apply(fundef,[f])))[1],theta:args(lhs(apply(fundef,[f])))[2],lam:BesselJZeros(0,n)/a,p:%pi,ff,k],local(ff), define(funmake(ff,[r,theta]),apply(f,[r,theta])*r*bessel_j(0,lam*r)*sin(m*theta)), k:2/(p*a^2*(bessel_j(m+1,lam))^2), k*nintegrate(ff,14,[0,a],[0,2*p]) )$ as0(n,a,c,f):=block([r:args(lhs(apply(fundef,[f])))[1],theta:args(lhs(apply(fundef,[f])))[2],lam:BesselJZeros(0,n)/a,p:%pi,ff,k],local(ff), define(funmake(ff,[r,theta]),apply(f,[r,theta])*r*bessel_j(0,lam*r)), k:1/(c*p*lam*a^2*(bessel_j(1,lam))^2), k*nintegrate(ff,14,[0,a],[0,2*p]) )$ as(m,n,a,c,f):=block([r:args(lhs(apply(fundef,[f])))[1],theta:args(lhs(apply(fundef,[f])))[2],lam:BesselJZeros(0,n)/a,p:%pi,ff,k],local(ff), define(funmake(ff,[r,theta]),apply(f,[r,theta])*r*bessel_j(0,lam*r)*cos(m*theta)), k:2/(c*p*lam*a^2*(bessel_j(m+1,lam))^2), k*nintegrate(ff,14,[0,a],[0,2*p]) )$ bs(m,n,a,c,f):=block([r:args(lhs(apply(fundef,[f])))[1],theta:args(lhs(apply(fundef,[f])))[2],lam:BesselJZeros(0,n)/a,p:%pi,ff,k],local(ff), define(funmake(ff,[r,theta]),apply(f,[r,theta])*r*bessel_j(0,lam*r)*sin(m*theta)), k:2/(c*p*lam*a^2*(bessel_j(m+1,lam))^2), k*nintegrate(ff,14,[0,a],[0,2*p]) )$ wave2d_disk(c,a,f,g,k,l):=block([A0,A,B,As0,As,Bs,lamb],local(A0,A,B,As0,As,Bs,lamb,t), A0[i]:=a0(i,a,f),A[i,j]:=a(i,j,a,f),B[i,j]:=b(i,j,a,f),As0[i]:=as0(i,a,c,g),As[i,j]:=as(i,j,a,c,g),Bs[i,j]:=bs(i,j,a,c,g), lamb[i,j]:=BesselJZeros(i,j)/a, sum(A0[n]*bessel_j(0,lamb[0,n]*r)*cos(c*lamb[0,n]*t),n,1,l) +sum(sum(bessel_j(m,lamb[m,n]*r)*(A[m,n]*cos(m*theta)+B[m,n]*sin(m*theta))*cos(c*lamb[m,n]*t),n,1,l),m,1,k) +sum(As0[n]*bessel_j(0,lamb[0,n]*r)*sin(c*lamb[0,n]*t),n,1,l) +sum(sum(bessel_j(m,lamb[m,n]*r)*(As[m,n]*cos(m*theta)+Bs[m,n]*sin(m*theta))*sin(c*lamb[m,n]*t),n,1,l),m,1,k) )$ Arect(f,a,b,m,n):=block([x:args(lhs(apply(fundef,[f])))[1],y:args(lhs(apply(fundef,[f])))[2]], define(funmake(ff,[x,y]),apply(f,[x,y])*sin(m*%pi*x/a)*sin(m*%pi*y/b)), 4*(nintegrate(ff,14,[0,a],[0,b]))/(a*b) )$ Brect(g,a,b,c,m,n):=block([x:args(lhs(apply(fundef,[f])))[1],y:args(lhs(apply(fundef,[f])))[2],lamb],local(lamb), lamb[m,n]:=((m/a)^2+(n/b)^2)/(%pi^2), define(funmake(gg,[x,y]),apply(g,[x,y])*sin(m*%pi*x/a)*sin(m*%pi*y/b)), 4*(nintegrate(gg,14,[0,a],[0,b]))/(a*b*c*sqrt(lam[m,n])) )$ wave2d_rectangle(c,a,b,f,g,k,l):=block([x:args(lhs(apply(fundef,[f])))[1],y:args(lhs(apply(fundef,[f])))[2],lamb,t],local(lamb,t), lamb[m,n]:=((m/a)^2+(n/b)^2)/(%pi^2), sum(sum((Arect(f,a,b,m,n)*cos(sqrt(lamb[0,n])*c*t)+Brect(g,a,b,c,m,n)*sin(sqrt(lamb[0,n])*c*t))*sin(m*%pi*x/a)*sin(n*%pi*y/b),m,1,k),n,1,l) )$