Application Center - Maplesoft

App Preview:

Fourier Method for Wave Equation

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

Learn about Maple
Download Application


 

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

eq := diff(u(x,t),`$`(t,2))-diff(u(x,t),`$`(x,2)) = 0

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

ic1 := u(x,0) = phi(x)

ic2 := D[2](u)(x,0) = psi(x)

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

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

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

phi := proc (x) options operator, arrow; piecewise(x < 1,x,x < 2,2-x) end proc

>    psi:=0;

psi := 0

>    plot(phi,0..2);

[Maple Plot]

>   

Solving eigenvalue problem

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

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

>    assume(lambda>0);

L := -diff(y(x),`$`(x,2)) = lambda*y(x)

s1 := y(0) = 0

s2 := y(2) = 0

>    dsolve(L,y(x));

y(x) = _C1*sin(sqrt(lambda)*x)+_C2*cos(sqrt(lambda)*x)

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

y := proc (x) options operator, arrow; _C1*sin(sqrt(lambda)*x)+_C2*cos(sqrt(lambda)*x) end proc

>    sist:={s1,s2};   

sist := {_C2 = 0, _C1*sin(2*sqrt(lambda))+_C2*cos(2*sqrt(lambda)) = 0}

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

matrix([[0, 1], [sin(2*sqrt(lambda)), cos(2*sqrt(lambda))]])

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

chl := -sin(2*sqrt(lambda)) = 0

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

1/4*Pi^2*_Z1^2

>    indets(%) minus {l};

{_Z1}

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

>    ev:=unapply(%,k);

ev := proc (k) options operator, arrow; 1/4*Pi^2*k^2 end proc

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

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

-diff(y(x),`$`(x,2)) = 1/4*Pi^2*k^2*y(x)

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

y(x) = _C1*sin(1/2*Pi*k*x)

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

1/2*_C1*sin(1/2*Pi*k*x)*2^(1/2)/(_C1^2)^(1/2)

>    simplify(%,radical,symbolic):

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

Eigenvalues and eigenfunctions:

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

1/4*Pi^2*k^2

1/2*sin(1/2*Pi*k*x)*sqrt(2)

>   

Solving problem

Solution of this problem is sought in in the form:

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

so := Sum(1/2*T[k](t)*sin(1/2*Pi*k*x)*sqrt(2),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)};

ode := {D(T[k])(0) = 0, T[k](0) = 4*2^(1/2)*sin(1/2*Pi*k)/Pi^2/k^2, diff(T[k](t),`$`(t,2))+1/4*Pi^2*k^2*T[k](t) = 0}

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

T[k](t) = 4*2^(1/2)*sin(1/2*Pi*k)/Pi^2/k^2*cos(1/2*Pi*k*t)

>   

Solution

>    sol:=subs(%,so);

sol := Sum(4*sin(1/2*Pi*k)/Pi^2/k^2*cos(1/2*Pi*k*t)*sin(1/2*Pi*k*x),k = 1 .. infinity)

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

sol_30 := 4/Pi^2*cos(1/2*Pi*t)*sin(1/2*Pi*x)-4/9*1/Pi^2*cos(3/2*Pi*t)*sin(3/2*Pi*x)+4/25/Pi^2*cos(5/2*Pi*t)*sin(5/2*Pi*x)-4/49*1/Pi^2*cos(7/2*Pi*t)*sin(7/2*Pi*x)+4/81/Pi^2*cos(9/2*Pi*t)*sin(9/2*Pi*x)-4...
sol_30 := 4/Pi^2*cos(1/2*Pi*t)*sin(1/2*Pi*x)-4/9*1/Pi^2*cos(3/2*Pi*t)*sin(3/2*Pi*x)+4/25/Pi^2*cos(5/2*Pi*t)*sin(5/2*Pi*x)-4/49*1/Pi^2*cos(7/2*Pi*t)*sin(7/2*Pi*x)+4/81/Pi^2*cos(9/2*Pi*t)*sin(9/2*Pi*x)-4...
sol_30 := 4/Pi^2*cos(1/2*Pi*t)*sin(1/2*Pi*x)-4/9*1/Pi^2*cos(3/2*Pi*t)*sin(3/2*Pi*x)+4/25/Pi^2*cos(5/2*Pi*t)*sin(5/2*Pi*x)-4/49*1/Pi^2*cos(7/2*Pi*t)*sin(7/2*Pi*x)+4/81/Pi^2*cos(9/2*Pi*t)*sin(9/2*Pi*x)-4...
sol_30 := 4/Pi^2*cos(1/2*Pi*t)*sin(1/2*Pi*x)-4/9*1/Pi^2*cos(3/2*Pi*t)*sin(3/2*Pi*x)+4/25/Pi^2*cos(5/2*Pi*t)*sin(5/2*Pi*x)-4/49*1/Pi^2*cos(7/2*Pi*t)*sin(7/2*Pi*x)+4/81/Pi^2*cos(9/2*Pi*t)*sin(9/2*Pi*x)-4...
sol_30 := 4/Pi^2*cos(1/2*Pi*t)*sin(1/2*Pi*x)-4/9*1/Pi^2*cos(3/2*Pi*t)*sin(3/2*Pi*x)+4/25/Pi^2*cos(5/2*Pi*t)*sin(5/2*Pi*x)-4/49*1/Pi^2*cos(7/2*Pi*t)*sin(7/2*Pi*x)+4/81/Pi^2*cos(9/2*Pi*t)*sin(9/2*Pi*x)-4...

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

[Maple Plot]

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

[Maple Plot]

Checking the Solution

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

>    simplify(%);

0 = 0

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

[Maple Plot]

>    plot(phi,0..2);

[Maple Plot]

>   

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;

eq := diff(u(x,t),`$`(t,2))-diff(u(x,t),`$`(x,2))+4*u(x,t) = 0

0 < x, x < 1, 0 < t

>    init_c:=u(x,0)=x^2-x,D[2](u)(x,0)=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;

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

L := -diff(y(x),`$`(x,2)) = lambda*y(x)

s1 := y(0) = 0

s2 := y(1) = 0

>    dsolve(L,y(x));

y(x) = _C1*sin(sqrt(lambda)*x)+_C2*cos(sqrt(lambda)*x)

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

y := proc (x) options operator, arrow; _C1*sin(sqrt(lambda)*x)+_C2*cos(sqrt(lambda)*x) end proc

>    sist:={s1,s2};   

sist := {_C2 = 0, _C1*sin(sqrt(lambda))+_C2*cos(sqrt(lambda)) = 0}

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

matrix([[0, 1], [sin(sqrt(lambda)), cos(sqrt(lambda))]])

>    chl:=det(%);      

chl := -sin(sqrt(lambda))

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

Pi^2*_Z1^2

>    indets(%) minus {l};

{_Z1}

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

>    tr:=unapply(%,k);

tr := proc (k) options operator, arrow; Pi^2*k^2 end proc

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

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

-diff(y(x),`$`(x,2)) = Pi^2*k^2*y(x)

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

y(x) = _C1*sin(Pi*k*x)

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

_C1*sin(Pi*k*x)*2^(1/2)/(_C1^2)^(1/2)

>    simplify(%,radical,symbolic):

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

Eigenvalues and eigenfunctions:

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

Pi^2*k^2

sin(Pi*k*x)*sqrt(2)

>   

Solving problem

Solution of this problem is sought in in the form:

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

spr := Sum(T[k](t)*sin(Pi*k*x)*sqrt(2),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)};

uz := {T[k](0) = 2*2^(1/2)*((-1)^k-1)/Pi^3/k^3, diff(T[k](t),`$`(t,2))+Pi^2*k^2*T[k](t)+4*T[k](t) = 0, D(T[k])(0) = 0}

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

T[k](t) = 2*2^(1/2)*((-1)^k-1)/Pi^3/k^3*cos(sqrt(Pi^2*k^2+4)*t)

>   

Solution

>    sol:=subs(%,spr);

sol := Sum(4*((-1)^k-1)/Pi^3/k^3*cos(sqrt(Pi^2*k^2+4)*t)*sin(Pi*k*x),k = 1 .. infinity)

>   

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;

eq := diff(u(x,t),`$`(t,2))-a^2*diff(u(x,t),`$`(x,2)) = A

0 < x, x < l, 0 < t

>    init_c:=u(x,0)=0,D[2](u)(x,0)=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;

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

L := -diff(y(x),`$`(x,2)) = lambda*y(x)

s1 := y(0) = 0

s2 := D(y)(l) = 0

>    dsolve(L,y(x));

y(x) = _C1*sin(sqrt(lambda)*x)+_C2*cos(sqrt(lambda)*x)

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

y := proc (x) options operator, arrow; _C1*sin(sqrt(lambda)*x)+_C2*cos(sqrt(lambda)*x) end proc

>    sist:={s1,s2};   

sist := {_C2 = 0, _C1*cos(sqrt(lambda)*l)*sqrt(lambda)-_C2*sin(sqrt(lambda)*l)*sqrt(lambda) = 0}

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

matrix([[0, 1], [cos(sqrt(lambda)*l)*sqrt(lambda), -sin(sqrt(lambda)*l)*sqrt(lambda)]])

>    chl:=det(%);      

chl := -cos(sqrt(lambda)*l)*sqrt(lambda)

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

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

1/4*Pi^2*(1+2*_Z1)^2/l^2

>    indets(%) minus {l};

{_Z1}

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

>    tr:=unapply(%,k);

tr := proc (k) options operator, arrow; 1/4*Pi^2*(1+2*k)^2/l^2 end proc

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

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

-diff(y(x),`$`(x,2)) = 1/4*Pi^2*(1+2*k)^2/l^2*y(x)

>    dsolve(%,y(x));

y(x) = _C1*sin(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)+_C2*cos(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)

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

y := proc (x) options operator, arrow; _C1*sin(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)+_C2*cos(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x) end proc

>    s1,s2;

_C2 = 0, 1/2*_C1*cos(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)-1/2*_C2*sin(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2) = 0
_C2 = 0, 1/2*_C1*cos(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)-1/2*_C2*sin(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2) = 0

>    subs(s1,y(x));

_C1*sin(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)

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

sin(1/2/l*Pi*(1+2*k)*x)

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

sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2)

>    #simplify(%,radical,symbolic):

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

Eigenvalues and eigenfunctions:

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

1/4/l^2*Pi^2*(1+2*k)^2

sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2)

>   

Particular solution

We find particular solution of problem:

>    eq; bound_c;

diff(u(x,t),`$`(t,2))-a^2*diff(u(x,t),`$`(x,2)) = A

u(0,t) = 0, D[1](u)(l,t) = 0

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

>   

diff(v(x),`$`(t,2))-a^2*diff(v(x),`$`(x,2)) = A

-a^2*diff(v(x),`$`(x,2)) = A

v(x) = -1/2*A/a^2*x^2+_C1*x+_C2

as := -1/2*A/a^2*x^2+_C1*x+_C2

l1 := _C2

l2 := -A/a^2*l+_C1

{_C2 = 0, _C1 = A/a^2*l}

psol := 1/2*x*A*(-x+2*l)/a^2

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;

eqn := diff(v(x,t),`$`(t,2))-a^2*diff(v(x,t),`$`(x,2)) = 0

0 < x, x < l, 0 < t

init_c := v(x,0) = -1/2*x*A*(-x+2*l)/a^2, 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);

spr := Sum(T[k](t)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),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};

uz := {D(T[k])(0) = 0, T[k](0) = -8*l^(5/2)/Pi^3/(1+2*k)^3*A/a^2*2^(1/2), diff(T[k](t),`$`(t,2))+1/4*a^2/l^2*Pi^2*(1+2*k)^2*T[k](t) = 0}

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

T[k](t) = -8*l^(5/2)*A*2^(1/2)/Pi^3/a^2/(1+6*k+12*k^2+8*k^3)*cos(1/2/l*a*Pi*(1+2*k)*t)

>    subs(%,spr);

Sum(-16*l^2*A/Pi^3/a^2/(1+6*k+12*k^2+8*k^3)*cos(1/2/l*a*Pi*(1+2*k)*t)*sin(1/2/l*Pi*(1+2*k)*x),k = 0 .. infinity)

>    factor(%);

Sum(-16*l^2*A*cos(1/2/l*a*Pi*(1+2*k)*t)*sin(1/2/l*Pi*(1+2*k)*x)/Pi^3/a^2/(1+2*k)^3,k = 0 .. infinity)

>   

Solution

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

sol := 1/2*x*A*(-x+2*l)/a^2+Sum(-16*l^2*A*cos(1/2/l*a*Pi*(1+2*k)*t)*sin(1/2/l*Pi*(1+2*k)*x)/Pi^3/a^2/(1+2*k)^3,k = 0 .. infinity)

>   

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;

eq := diff(u(x,t),`$`(t,2))+2*beta*diff(u(x,t),t)+beta^2*u(x,t)-a^2*diff(u(x,t),`$`(x,2))

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;

L := -a^2*diff(y(x),`$`(x,2)) = lambda*y(x)

s1 := y(0) = 0

s2 := D(y)(l) = 0

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

>    dsolve(L,y(x));

y(x) = _C1*sin(1/a*lambda^(1/2)*x)+_C2*cos(1/a*lambda^(1/2)*x)

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

y := proc (x) options operator, arrow; _C1*sin(1/a*lambda^(1/2)*x)+_C2*cos(1/a*lambda^(1/2)*x) end proc

>    sist:={s1,s2};   

sist := {_C1*cos(1/a*lambda^(1/2)*l)/a*lambda^(1/2)-_C2*sin(1/a*lambda^(1/2)*l)/a*lambda^(1/2) = 0, _C2 = 0}

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

matrix([[cos(1/a*lambda^(1/2)*l)/a*lambda^(1/2), -sin(1/a*lambda^(1/2)*l)/a*lambda^(1/2)], [0, 1]])

>    chl:=det(%);      

chl := cos(1/a*lambda^(1/2)*l)/a*lambda^(1/2)

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

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

1/4*Pi^2*(1+2*_Z1)^2*a^2/l^2

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

{_Z1}

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

>    tr:=unapply(%,k);

tr := proc (k) options operator, arrow; 1/4*Pi^2*(1+2*k)^2*a^2/l^2 end proc

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

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

-a^2*diff(y(x),`$`(x,2)) = 1/4*Pi^2*(1+2*k)^2*a^2/l^2*y(x)

>    dsolve(%,y(x));

y(x) = _C1*sin(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)+_C2*cos(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)

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

y := proc (x) options operator, arrow; _C1*sin(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)+_C2*cos(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x) end proc

>    s1,s2;

_C2 = 0, 1/2*_C1*cos(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)-1/2*_C2*sin(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2) = 0
_C2 = 0, 1/2*_C1*cos(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)-1/2*_C2*sin(1/2*sqrt(Pi^2+4*Pi^2*k+4*Pi^2*k^2))/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2) = 0

>    subs(s1,y(x));

_C1*sin(1/2/l*(Pi^2+4*Pi^2*k+4*Pi^2*k^2)^(1/2)*x)

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

sin(1/2/l*Pi*(1+2*k)*x)

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

sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2)

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

Eigenvalues and eigenfunctions:

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

1/4*Pi^2*(1+2*k)^2*a^2/l^2

sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2)

>   

Solving problem

Solution of this problem is sought in in the form:

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

spr := Sum(T[k](t)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),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)};

uz := {T[k](0) = Int(phi(x)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),x = 0 .. l), D(T[k])(0) = Int(-beta*phi(x)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),x = 0 .. l), diff(T[k](t),`$`(t,2))+2*beta*diff(T[...
uz := {T[k](0) = Int(phi(x)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),x = 0 .. l), D(T[k])(0) = Int(-beta*phi(x)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),x = 0 .. l), diff(T[k](t),`$`(t,2))+2*beta*diff(T[...
uz := {T[k](0) = Int(phi(x)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),x = 0 .. l), D(T[k])(0) = Int(-beta*phi(x)*sin(1/2/l*Pi*(1+2*k)*x)*2^(1/2)/l^(1/2),x = 0 .. l), diff(T[k](t),`$`(t,2))+2*beta*diff(T[...

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

>   

T[k](t) = 2^(1/2)/l^(1/2)*Int(phi(x)*sin(1/2/l*Pi*(1+2*k)*x),x = 0 .. l)*exp(-beta*t)*cos(1/2*Pi/l*a*(1+2*k)*t)

>    subs(%,spr);

Sum(2/l*Int(phi(x)*sin(1/2/l*Pi*(1+2*k)*x),x = 0 .. l)*exp(-beta*t)*cos(1/2*Pi/l*a*(1+2*k)*t)*sin(1/2/l*Pi*(1+2*k)*x),k = 0 .. infinity)

>    simplify(%);

Sum(2/l*Int(phi(x)*sin(1/2/l*Pi*(1+2*k)*x),x = 0 .. l)*exp(-beta*t)*cos(1/2*Pi/l*a*(1+2*k)*t)*sin(1/2/l*Pi*(1+2*k)*x),k = 0 .. infinity)

>    simplify(sort(%));

Sum(2/l*cos(1/2*Pi*(2*k+1)*a/l*t)*Int(phi(x)*sin(1/2*Pi*(2*k+1)/l*x),x = 0 .. l)*sin(1/2*Pi*(2*k+1)/l*x)*exp(-beta*t),k = 0 .. infinity)

>   

Solution

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

sol := Sum(2/l*cos(1/2*Pi*(2*k+1)*a/l*t)*Int(phi(x)*sin(1/2*Pi*(2*k+1)/l*x),x = 0 .. l)*sin(1/2*Pi*(2*k+1)/l*x)*exp(-beta*t),k = 0 .. infinity)

>   

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

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

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

lc := u(0,t) = 0

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

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

l := 1/2*Pi

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

ic := u(x,0) = 0

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

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

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

rs := sin(3*x)+sin(x)

Solving

With not Fourier method .   sol= sol1+sol2 .

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

sol1 := u(x,t) = sin(3*x)*w(t)

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

l1 := diff(sin(3*x)*w(t),`$`(t,2))-diff(sin(3*x)*w(t),`$`(x,2))-10*sin(3*x)*w(t) = sin(3*x)

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

w(t) = 1/2*exp(t)+1/2*exp(-t)-1

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

sol1 := sin(3*x)*(1/2*exp(t)+1/2*exp(-t)-1)

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

sol2 := u(x,t) = sin(x)*w(t)

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

l2 := diff(sin(x)*w(t),`$`(t,2))-diff(sin(x)*w(t),`$`(x,2))-10*sin(x)*w(t) = sin(x)

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

w(t) = 1/18*exp(3*t)+1/18*exp(-3*t)-1/9

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

sol2 := sin(x)*(1/18*exp(3*t)+1/18*exp(-3*t)-1/9)

Solution

>    sol:=sol1+sol2;

sol := sin(3*x)*(1/2*exp(t)+1/2*exp(-t)-1)+sin(x)*(1/18*exp(3*t)+1/18*exp(-3*t)-1/9)

>   

Checking the Solution

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

u := proc (x, t) options operator, arrow; sin(3*x)*(1/2*exp(t)+1/2*exp(-t)-1)+sin(x)*(1/18*exp(3*t)+1/18*exp(-3*t)-1/9) end proc

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

0

0

0

0

0

>   

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;

eq := diff(u(r,t),`$`(t,2))-a^2*(1/r*diff(u(r,t),r)+diff(u(r,t),`$`(r,2))) = 0

0 < r, r < R, 0 < t

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

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

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

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

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

bound_c := u(R,t) = 0

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

a := 2

R := 5

u0 := proc (r) options operator, arrow; 1/8-1/200*r^2 end proc

u1 := 0

Solution representation formula

u(r,t) = 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 .. infinity)

where   mu[k]  - Zeros of Bessel J(0,x) function   BesselJ(0,x)

A[k] = 2*int(r*u0(r)*BesselJ(0,mu[k]*r/R),r = 0 .. R)/(R^2*BesselJ(1,mu[k])^2), B[k] = 2*int(r*u1(r)*BesselJ(0,mu[k]*r/R),r = 0 .. R)/(R^2*BesselJ(1,mu[k])^2)

Zeros of Bessel J(0,x) function

>    Order:=10;

Order := 10

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

J := proc (x) options operator, arrow; BesselJ(0,x) end proc

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

[Maple Plot]

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

mu[1] := 2.404825558

mu[2] := 5.520078110

mu[3] := 8.653727913

mu[4] := 11.79153444

mu[5] := 14.93091771

mu[6] := 18.07106397

mu[7] := 21.21163663

mu[8] := 24.35247153

mu[9] := 27.49347913

mu[10] := 30.63460647

mu[11] := 33.77582021

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

A := proc (k) options operator, arrow; 2/R^2/BesselJ(1,mu[k])^2*int(r*u0(r)*BesselJ(0,mu[k]/R*r),r = 0 .. R) end proc

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

B := proc (k) options operator, arrow; 2/R^2/BesselJ(1,mu[k])^2*int(r*u1(r)*BesselJ(0,mu[k]/R*r),r = 0 .. R) end proc

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

sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...
sol := .1385027827*cos(.9619302232*t)*BesselJ(0,.4809651116*r)-.1747218817e-1*cos(2.208031244*t)*BesselJ(0,1.104015622*r)+.5684558837e-2*cos(3.461491165*t)*BesselJ(0,1.730745583*r)-.2623862777e-2*cos(4...

>   

Plots

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

[Maple Plot]

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

[Maple Plot]

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

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

m := 16

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

f := proc (x) options operator, arrow; cat(`t=`,convert(x,symbol)) end proc

>    tit:=map(f,K);

tit := [`t=0.`, `t=1.0`, `t=1.5`, `t=1.7`, `t=1.8`, `t=2.0`, `t=2.5`, `t=3.0`, `t=3.5`, `t=4.0`, `t=4.5`, `t=4.8`, `t=4.85`, `t=4.9`, `t=5.5`, `t=6.0`]
tit := [`t=0.`, `t=1.0`, `t=1.5`, `t=1.7`, `t=1.8`, `t=2.0`, `t=2.5`, `t=3.0`, `t=3.5`, `t=4.0`, `t=4.5`, `t=4.8`, `t=4.85`, `t=4.9`, `t=5.5`, `t=6.0`]

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

[Maple Plot]

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

[Maple Plot]

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

[Maple Plot]

>   

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