Application Center - Maplesoft

# Fourier Method for Wave Equation

You can switch back to the summary page by clicking here.

wave_eqf.mws

Fourier method for wave equation

by Aleksas Domarkas

Vilnius University, Faculty of Mathematics and Informatics,

Naugarduko 24, Vilnius, Lithuania

aleksas@ieva.mif.vu.lt

NOTE: In this session we find  solutions of boundary value problems for wave equation using Fourier method.

1 Example

Problem

 > restart;

 > eq:=diff(u(x,t),t,t)-diff(u(x,t),x,x)=0;#00;

 > ic1:=u(x,0)=phi(x);#0

 > bound_c:=u(0,t) = 0,u(2,t)=0;

 > phi:=x->piecewise(x<1,x,x<2,2-x);

 > psi:=0;

 > plot(phi,0..2);

 >

Solving eigenvalue problem

 > L:=-diff(y(x),x,x)=lambda*y(x);

 > s1:=y(0)=0; s2:=y(2)=0;

 > assume(lambda>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > linalg[genmatrix](sist,{_C1,_C2});

 > chl:=linalg[det](%)=0;

 > #chl:=combine(%);

 > _EnvAllSolutions := true:

 > solve(chl,lambda);

 > indets(%) minus {l};

 > subs(%[1]='k',%%):

 > ev:=unapply(%,k);

 > y:='y':assume(k,posint);

 > subs(lambda=ev(k),L);

 > dsolve({%,s1,s2},y(x));

 > rhs(%)/sqrt(int(rhs(%)^2,x=0..4));

 > ef:=unapply(%,(k,x)):

Eigenvalues and eigenfunctions:

 > ev(k);ef(k,x);

 >

Solving problem

Solution of this problem is sought in in the form:

 > so:=Sum(T[k](t)*ef(k,x),k=1..infinity);

For coefficients  T[k](t) we have problem:

 > ode:={diff(T[k](t),t,t)+ev(k)*T[k](t)=0,

 > T[k](0)=int(phi(x)*ef(k,x),x=0..4),

 > D(T[k])(0)=int(psi(x)*ef(k,x),x=0..4)};

 > simplify(dsolve(ode,T[k](t)));

 >

Solution

 > sol:=subs(%,so);

 > sol_30:=value(subs(infinity=30,sol));

 > plot3d(sol_30,x=0..2,t=0..10,axes=boxed,orientation=[-30,75]);

 > plots[animate](sol_30,x=0..2,t=0..4,frames=80);

Checking the Solution

 > subs(u(x,t)=sol_30,eq):

 > simplify(%);

 > plot(subs(t=0,sol_30),x=0..2);

 > plot(phi,0..2);

 >

2 Example

Problem

 > restart;with(linalg):

```Warning, the protected names norm and trace have been redefined and unprotected
```

 > eq:=diff(u(x,t),t,t)-diff(u(x,t),x,x)+4*u(x,t)=0;00;

 > init_c:=u(x,0)=x^2-x,D[2](u)(x,0)=0;

 > bound_c:=u(0,t) = 0,u(1,t)=0;

 > ic1:=rhs(init_c[1]):ic2:=rhs(init_c[2]):dp:=rhs(eq):

Solving eigenvalue problem

 > L:=-diff(y(x),x,x)=lambda*y(x);

 > s1:=y(0)=0; s2:=y(1)=0;

 > assume(lambda>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > genmatrix(sist,{_C1,_C2});

 > chl:=det(%);

 > #chl:=combine(%);

 > _EnvAllSolutions := true:

 > solve(chl,lambda);

 > indets(%) minus {l};

 > subs(%[1]='k',%%):

 > tr:=unapply(%,k);

 > y:='y':assume(k,posint);

 > subs(lambda=tr(k),L);

 > dsolve({%,s1,s2},y(x));

 > rhs(%)/sqrt(int(rhs(%)^2,x=0..1));

 > tf:=unapply(%,(k,x)):

Eigenvalues and eigenfunctions:

 > tr(k);tf(k,x);

 >

Solving problem

Solution of this problem is sought in in the form:

 > spr:=Sum(T[k](t)*tf(k,x),k=1..infinity);

For coefficients  T[k](t) we have problem:

 > uz:={diff(T[k](t),t,t)+tr(k)*T[k](t)+4*T[k](t)=int(dp*tf(k,x),x=0..1),

 > T[k](0)=int(ic1*tf(k,x),x=0..1),

 > D(T[k])(0)=int(ic2*tf(k,x),x=0..1)};

 > simplify(dsolve(%,T[k](t)));

 >

Solution

 > sol:=subs(%,spr);

 >

3 Example

Problem

 > restart;with(linalg):

```Warning, the protected names norm and trace have been redefined and unprotected
```

 > eq:=diff(u(x,t),t,t)-a^2*diff(u(x,t),x,x)=A;00;

 > init_c:=u(x,0)=0,D[2](u)(x,0)=0;

 > bound_c:=u(0,t)=0,D[1](u)(l,t)=0;

 > #l>0, A -- constant;

Solving eigenvalue problem

 > L:=-diff(y(x),x,x)=lambda*y(x);

 > s1:=y(0)=0; s2:=D(y)(l)=0;

 > assume(lambda>0,l>0,a>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > genmatrix(sist,{_C1,_C2});

 > chl:=det(%);

 > #chl:=combine(%);

 > _EnvAllSolutions := true:

 > solve(chl/sqrt(lambda),lambda);

 > indets(%) minus {l};

 > subs(%[1]='k',%%):

 > tr:=unapply(%,k);

 > y:='y':assume(k,posint);

 > subs(lambda=tr(k),L);

 > dsolve(%,y(x));

 > y:=unapply(rhs(%),x);

 > s1,s2;

 > subs(s1,y(x));

 > simplify(select(has,%,x));

 > %/sqrt(int(%^2,x=0..l));

 > tf:=unapply(%,(k,x)):

Eigenvalues and eigenfunctions:

 > tr(k);tf(k,x);

 >

Particular solution

We find particular solution of problem:

 > eq; bound_c;

 > subs(u(x,t)=v(x),eq);

 > simplify(%);

 > dsolve(%,v(x));

 > as:=rhs(%);

 > l1:=subs(x=0,%);

 > l2:=subs(x=l,diff(%%,x));

 > solve({l1,l2},{_C1,_C2});

 > psol:=factor(subs(%,as));

 >

Solving problem

We make change: u(x,t)=v(x,t)+psol. Then for v(x,t) we have problem:

 > eqn:=diff(v(x,t),t,t)-a^2*diff(v(x,t),x,x)=0;00;

 > init_c:=v(x,0)=-psol,D[2](v)(x,0)=0;

 > bound_c:=v(0,t)=0,D[1](v)(l,t)=0;

 >

Solution of this problem is sought in in the form:

 > spr:=Sum(T[k](t)*tf(k,x),k=0..infinity);

For coefficients  T[k](t) we have problem:

 > uz:={diff(T[k](t),t,t)+a^2*tr(k)*T[k](t)=0,

 > T[k](0)=int(-psol*tf(k,x),x=0..l),

 > D(T[k])(0)=0};

 > simplify(dsolve(%,T[k](t)));

 > subs(%,spr);

 > factor(%);

 >

Solution

 > sol:=psol+subs(a='a',k='k',l='l',h='h',c='c',%);

 >

4 Example

Problem

 > restart;with(linalg):

```Warning, the protected names norm and trace have been redefined and unprotected
```

 > eq:=diff(u(x,t),t,t)+2*beta*diff(u(x,t),t)

 > +beta^2*u(x,t)-a^2*diff(u(x,t),x,x);

 > init_c:=u(x,0)=phi(x),D[2](u)(x,0)=-beta*phi(x);

 > bound_c:=u(0,t)=0,D[1](u)(l,t)=0;

Solving eigenvalue problem

 > L:=-a^2*diff(y(x),x,x)=lambda*y(x);

 > s1:=y(0)=0; s2:=D(y)(l)=0;

 > assume(lambda>0,l>0,a>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > genmatrix(sist,{_C1,_C2});

 > chl:=det(%);

 > #chl:=combine(%);

 > _EnvAllSolutions := true:

 > solve(chl/sqrt(lambda),lambda);

 > indets(%) minus {l,a};

 > subs(%[1]='k',%%):

 > tr:=unapply(%,k);

 > y:='y':assume(k,posint);

 > subs(lambda=tr(k),L);

 > dsolve(%,y(x));

 > y:=unapply(rhs(%),x);

 > s1,s2;

 > subs(s1,y(x));

 > simplify(select(has,%,x));

 > %/sqrt(int(%^2,x=0..l));

 > tf:=unapply(%,(k,x)):

Eigenvalues and eigenfunctions:

 > tr(k);tf(k,x);

 >

Solving problem

Solution of this problem is sought in in the form:

 > spr:=Sum(T[k](t)*tf(k,x),k=0..infinity);

For coefficients  T[k](t) we have problem:

 > uz:={diff(T[k](t),t,t)+2*beta*diff(T[k](t),t)+beta^2*T[k](t)+tr(k)*T[k](t)=0,

 > T[k](0)=Int(phi(x)*tf(k,x),x=0..l),D(T[k])(0)=Int(-beta*phi(x)*tf(k,x),x=0..l)};

 > simplify(dsolve(%,T[k](t)));

 >

 > subs(%,spr);

 > simplify(%);

 > simplify(sort(%));

 >

Solution

 > sol:=subs(a='a',k='k',l='l',h='h',c='c',%);

 >

5 Example

Problem

 > restart;

 > eq:=diff(u(x,t),t,t)-diff(u(x,t),x,x)-10*u(x,t)=2*sin(2*x)*cos(x);

 > lc:=u(0,t)=0;

 > rc:=D[1](u)(Pi/2,t)=0;l:=Pi/2;

 > ic:=u(x,0)=0;

 > ict:=D[2](u)(x,0)=0;

 > rs:=combine(rhs(eq));

Solving

With not Fourier method .   sol= sol1+sol2 .

 > sol1:=u(x,t)=op(1,rs)*w(t);

 > l1:=subs(sol1,lhs(eq)=op(1,rs));

 > dsolve({l1,w(0)=0,D(w)(0)=0},w(t));

 > sol1:=subs(%,rhs(sol1));

 > sol2:=u(x,t)=op(2,rs)*w(t);

 > l2:=subs(sol2,lhs(eq)=op(2,rs));

 > dsolve({l2,w(0)=0,D(w)(0)=0},w(t));

 > sol2:=subs(%,rhs(sol2));

Solution

 > sol:=sol1+sol2;

 >

Checking the Solution

 > u:=unapply(sol,(x,t));

 > simplify(lhs(eq)-rhs(eq));

 > simplify(lhs(lc)-rhs(lc));

 > simplify(lhs(rc)-rhs(rc));

 > simplify(lhs(ic)-rhs(ic));

 > simplify(lhs(ict)-rhs(ict));

 >

6   Example

Membrane   Problem

 > restart;with(linalg,laplacian):

 > eq:=diff(u(r,t),t,t)-a^2*expand(laplacian(u(r,t),[r,phi],coords=polar))=0;00;

 > ic:=u(r,0)=u0(r);

 > ict:=D[2](u)(r,0)=u1(r);

 > bound_c:=u(R,t)=0;

 > a:=2;R:=5;u0:=r->(1-(r/5)^2)/8;u1:=0;

Solution representation formula

where    - Zeros of Bessel J(0,x) function

Zeros of Bessel J(0,x) function

 > Order:=10;

 > J:=x->BesselJ(0,x);

 > plot(J(x),x=0..30); #plot(J(x),x=10..30);

 > mu[1]:=fsolve(J(x),x,2..3);

 > for k from 1 to Order do

 > mu[k+1]:=fsolve(J(x),x,%..%+2*Pi);od;

Note: see ?BesselJZeros

Solution

 > A:=k->2/(R*BesselJ(1,mu[k]))^2*int(r*u0(r)*BesselJ(0,mu[k]/R*r),r=0..R);

 > B:=k->2/(R*BesselJ(1,mu[k]))^2*int(r*u1(r)*BesselJ(0,mu[k]/R*r),r=0..R);

 > k:='k':

 > sol:=sum((A(k)*cos(a*mu[k]*t/R)+B(k)*sin(a*mu[k]*t/R)) *BesselJ(0,mu[k]*r/R),k=1..Order);

 >

Plots

 > plot(simplify(subs(r=0,sol)),t=0..10);

 > plots[animate](sol,r=-R..R,t=0..6,frames=20);

 > K := [0.,1.0,1.5,1.7,1.8,2.0,2.5,3.0,3.5,4.0,4.5,4.8,4.85,4.9,5.5,6.0];

 > m:=nops(%);

 > f:=x->cat(`t=`,convert(x,symbol));

 > tit:=map(f,K);

 > for k to m do P||k:=plot(subs(t=K[k],sol),r=-R..R,title=tit[k]):end do:

 > plots[display]([seq(P||k,k=1..m)],insequence=true);

 > plots[display](%,axes=none,labels=[` `,` `]);

3d animation:

 > for k to m do Pc||k:=plot3d(subs(t=K[k],sol),r=0..R,theta=0..2*Pi, coords=z_cylindrical, title=tit[k]):end do:

 > plots[display3d]([seq(Pc||k,k=1..m)],insequence=true);

 >

See:

or

Study of Membrane Oscillation Using Maple

by D. Frenkel , L. Golebiowski and R. Portugal

 >

 >

While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the conevibutors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.