• 1 Example
  • Problem
  • Solving eigenvalue problem
  • Solving problem
  • Solution
  • Checking the Solution
• 2 Example
  • Problem
  • Solving eigenvalue problem
  • Solving problem
  • Solution
• 3 Example
  • Problem
  • Solving eigenvalue problem
  • Particular solution
  • Solving problem
  • Solution
• 4 Example
  • Problem
  • Solving eigenvalue problem
  • Solving problem
  • Solution
• 5 Example
  • Problem
  • Solving
  • Solution
  • Checking the Solution
• 6 Example
  • Membrane Problem
  • Solution representation formula
  • Zeros of Bessel J(0,x) function
  • Solution
  • Plots
  • Other authors about Membrame Problem

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;#0<x,x<2,t>0;

>	ic1:=u(x,0)=phi(x);#0<x,x<2 ic2:=D[2](u)(x,0)=psi(x);#0<x,x<2

>	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));

>	simplify(%,radical,symbolic):

>	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

Vlad. 20.14.1

>	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;0<x,x<1,t>0;

>	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));

>	simplify(%,radical,symbolic):

>	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;0<x,x<l,t>0;

>	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));

>	#simplify(%,radical,symbolic):

>	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;0<x,x<l,t>0;

>	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

[vlad] 20.16.7.

>

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;0<r,r<R,t>0;

>	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:

>	addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]):

>	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);

>

Other authors about Membrame Problem

See:

http://www.mapleapps.com/categories/engineering/engmathematics/html/K11.8-membrane.html  by Alain Goriely:

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.

Back to contents