Application Center - Maplesoft
  Submit your work
Application Center Applications Fourier Method for Heat Equation Preview

App Preview:

Fourier Method for Heat Equation

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

Learn about Maple
Download Application

 

heat_eqf.mws

Fourier method for heat 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 heat equation using Fourier method

1 Example

Problem

>    restart;

>    eq:=diff(u(x,t),t)-a^2*diff(u(x,t),x,x)=0;#0<x,x<Pi,t>0;

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

>    init_c:=u(x,0)=phi(x);#0<x,x<Pi

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

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

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

>    phi:=x->Pi/2-abs(Pi/2-x);

phi := proc (x) options operator, arrow; 1/2*Pi-abs(1/2*Pi-x) end proc

>    'phi(x)'=convert(phi(x),piecewise);

phi(x) = PIECEWISE([x, x < 1/2*Pi],[Pi-x, 1/2*Pi <= x])

>    plot(phi,0..Pi,axes=boxed);

[Maple Plot]

>   

Eigenvalues and eigenfunctions

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

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

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

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

s1 := y(0) = 0

s2 := y(Pi) = 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)*Pi)+_C2*cos(sqrt(lambda)*Pi) = 0}

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

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

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

chl := -sin(sqrt(lambda)*Pi)

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

_Z1^2

>    indets(%) minus {l};

{_Z1}

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

>    ev:=unapply(%,k);

ev := proc (k) options operator, arrow; k^2 end proc

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

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

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

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

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

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

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

>    simplify(%,radical,symbolic):

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

Eigenvalues and eigenfunctions:

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

k^2

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

Solving problem

Solution we find in form:

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

spr := Sum(T[k](t)*sin(k*x)*2^(1/2)/Pi^(1/2),k = 1 .. infinity)

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

>    ode:={diff(u(t),t)+a^2*ev(k)*u(t)=0,

>    u(0)=int((phi(x))*ef(k,x),x=0..Pi)};

ode := {diff(u(t),t)+a^2*k^2*u(t) = 0, u(0) = 2*sin(1/2*Pi*k)*2^(1/2)/Pi^(1/2)/k^2}

>    dsolve(ode,u(t));

u(t) = 2*sin(1/2*Pi*k)*2^(1/2)/Pi^(1/2)/k^2*exp(-a^2*k^2*t)

Solution

>    sol:=subs(T[k](t)=rhs(%),spr);

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

>    sol_20:=value(subs(infinity=20,sol));

sol_20 := 4/Pi*exp(-a^2*t)*sin(x)-4/9*1/Pi*exp(-9*a^2*t)*sin(3*x)+4/25/Pi*exp(-25*a^2*t)*sin(5*x)-4/49*1/Pi*exp(-49*a^2*t)*sin(7*x)+4/81/Pi*exp(-81*a^2*t)*sin(9*x)-4/121*1/Pi*exp(-121*a^2*t)*sin(11*x)+...
sol_20 := 4/Pi*exp(-a^2*t)*sin(x)-4/9*1/Pi*exp(-9*a^2*t)*sin(3*x)+4/25/Pi*exp(-25*a^2*t)*sin(5*x)-4/49*1/Pi*exp(-49*a^2*t)*sin(7*x)+4/81/Pi*exp(-81*a^2*t)*sin(9*x)-4/121*1/Pi*exp(-121*a^2*t)*sin(11*x)+...
sol_20 := 4/Pi*exp(-a^2*t)*sin(x)-4/9*1/Pi*exp(-9*a^2*t)*sin(3*x)+4/25/Pi*exp(-25*a^2*t)*sin(5*x)-4/49*1/Pi*exp(-49*a^2*t)*sin(7*x)+4/81/Pi*exp(-81*a^2*t)*sin(9*x)-4/121*1/Pi*exp(-121*a^2*t)*sin(11*x)+...

>    plot3d(subs(a=1,sol_20),x=0..Pi,t=0..1);

[Maple Plot]

Test

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

>    simplify(%);

0 = 0

>    plot(subs(t=0,sol_20),x=0..Pi);

[Maple Plot]

>    plot(phi,0..Pi);

[Maple Plot]

>   

2 Example

Problem

>    restart;

>    eq:=diff(u(x,t),t)-a^2*diff(u(x,t),x,x)=0;#0<x,x<l,t>0;

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

>    init_c:=u(x,0)=phi(x);#0<x,x<l;

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

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

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

>    phi:=x->x*(l-x);

phi := proc (x) options operator, arrow; x*(l-x) end proc

>    #l>0;

Eigenvalues and eigenfunctions

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

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

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

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

s1 := y(0) = 0

s2 := 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 := {_C1*sin(sqrt(lambda)*l)+_C2*cos(sqrt(lambda)*l) = 0, _C2 = 0}

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

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

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

chl := sin(sqrt(lambda)*l)

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

Pi^2*_Z1^2/l^2

>    indets(%) minus {l};

{_Z1}

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

>    ev:=unapply(%,k);

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

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

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

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

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

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

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

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

>    simplify(%,radical,symbolic):

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

Eigenvalues and eigenfunctions:

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

Pi^2*k^2/l^2

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

Solving problem

Solution we find in form:

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

spr := Sum(T[k](t)*sin(Pi*k/l*x)*2^(1/2)/l^(1/2),k = 1 .. infinity)

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

>    ode:={diff(u(t),t)+a^2*ev(k)*u(t)=0,

>    u(0)=int((phi(x))*ef(k,x),x=0..l)};

ode := {u(0) = -2*l^(5/2)*2^(1/2)*((-1)^k-1)/Pi^3/k^3, diff(u(t),t)+a^2*Pi^2*k^2/l^2*u(t) = 0}

>    dsolve(ode,u(t));

u(t) = -2*l^(5/2)*2^(1/2)*((-1)^k-1)/Pi^3/k^3*exp(-a^2*Pi^2*k^2/l^2*t)

Solution

>    sol:=subs(T[k](t)=rhs(%),spr);

sol := Sum(-4*l^2*((-1)^k-1)/Pi^3/k^3*exp(-a^2*Pi^2*k^2/l^2*t)*sin(Pi*k/l*x),k = 1 .. infinity)

We convert this solution to other form. Index   k   change to   2*m+1  :

>    subs(k=2*m+1,op(1,%));assume(m,integer):

-4*l^2*((-1)^(2*m+1)-1)/Pi^3/(2*m+1)^3*exp(-a^2*Pi^2*(2*m+1)^2/l^2*t)*sin(Pi*(2*m+1)/l*x)

>    simplify(%):factor(%);

8*l^2/Pi^3/(2*m+1)^3*exp(-a^2*Pi^2*(2*m+1)^2/l^2*t)*sin(Pi*(2*m+1)/l*x)

>    Sum(%,m=0..infinity):

>    solution:=subs(a='a',k='k',l='l',m='m',%);

solution := Sum(8*l^2/Pi^3/(2*m+1)^3*exp(-a^2*Pi^2*(2*m+1)^2/l^2*t)*sin(Pi*(2*m+1)/l*x),m = 0 .. infinity)

>   

3 Example

Problem

[bic] 694

>    restart;

>    eq:=diff(u(x,t),t)+beta*u(x,t)-a^2*diff(u(x,t),x,x)=0;0<x,x<l,t>0;

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

0 < x, x < l, 0 < t

>    init_c:=u(x,0)=phi(x);

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

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

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

>    #l>0;

Eigenvalues and eigenfunctions

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

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

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

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

s1 := y(0) = 0

s2 := 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*sin(sqrt(lambda)*l)+_C2*cos(sqrt(lambda)*l) = 0}

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

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

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

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

>    #chl:=combine(%);   

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

Pi^2*_Z1^2/l^2

>    indets(%) minus {l};

{_Z1}

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

>    tr:=unapply(%,k);

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

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

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

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

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

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

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

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

>    simplify(%,radical,symbolic):

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

Eigenvalues and eigenfunctions:

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

Pi^2*k^2/l^2

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

Solving problem

Solution we find in form:

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

spr := Sum(T[k](t)*sin(Pi*k/l*x)*2^(1/2)/l^(1/2),k = 1 .. infinity)

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

>    uz:={diff(u(t),t)+beta*u(t)+a^2*tr(k)*u(t)=0,

>    u(0)=int((phi(xi))*tf(k,xi),xi=0..l)};

uz := {diff(u(t),t)+beta*u(t)+a^2*Pi^2*k^2/l^2*u(t) = 0, u(0) = int(phi(xi)*sin(Pi*k/l*xi)*2^(1/2)/l^(1/2),xi = 0 .. l)}

>    dsolve(uz,u(t));

u(t) = int(phi(xi)*sin(Pi*k/l*xi)*2^(1/2)/l^(1/2),xi = 0 .. l)*exp((-beta*l^2-k^2*Pi^2*a^2)*t/l^2)

Solution

>    sol:=simplify(subs(T[k](t)=rhs(%),spr));

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

or

>    subs(a='a',k='k',l='l',%);

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

>   

4 Example

Problem

[komec] 93 p.

>    restart;

>    eq:=diff(u(x,t),t)-a^2*diff(u(x,t),x,x)=0;0<x,x<l,t>0;

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

0 < x, x < l, 0 < t

>    init_c:=u(x,0)=phi(x);

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

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

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

>    a:=2;l:=3;phi:=x->x;

a := 2

l := 3

phi := proc (x) options operator, arrow; x end proc

Eigenvalues and eigenfunctions

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

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

>    assume(lambda>0);

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

s1 := y(0) = 0

s2 := D(y)(3) = 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 := {_C1*cos(3*sqrt(lambda))*sqrt(lambda)-_C2*sin(3*sqrt(lambda))*sqrt(lambda) = 0, _C2 = 0}

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

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

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

chl := cos(3*sqrt(lambda))*sqrt(lambda)

>    chl:=%/sqrt(lambda);   

chl := cos(3*sqrt(lambda))

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

1/36*Pi^2*(1+2*_Z1)^2

>    indets(%) minus {l};

{_Z1}

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

>    tr:=unapply(%,k);

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

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

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

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

>    dsolve(%,y(x));

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

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

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

>    solve(s1,_C2);

0

>    simplify(subs(_C2=%,y(x)));

_C1*sin(1/6*Pi*(1+2*k)*x)

>    select(has,%,x);

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

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

1/3*sin(1/6*Pi*(1+2*k)*x)*sqrt(6)

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

tf := proc (k, x) options operator, arrow; 1/3*sin(1/6*Pi*(1+2*k)*x)*sqrt(6) end proc

Eigenvalues and eigenfunctions:

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

1/36*Pi^2*(1+2*k)^2

1/3*sin(1/6*Pi*(1+2*k)*x)*sqrt(6)

Solving problem

Solution we find in form:

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

spr := Sum(1/3*T[k](t)*sin(1/6*Pi*(1+2*k)*x)*sqrt(6),k = 0 .. infinity)

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

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

>    u(0)=int((phi(xi))*tf(k,xi),xi=0..l)};

uz := {u(0) = 12*(-1)^k*6^(1/2)/Pi^2/(1+2*k)^2, diff(u(t),t)+1/9*Pi^2*(1+2*k)^2*u(t) = 0}

>    dsolve(uz,u(t));

u(t) = 12*(-1)^k*6^(1/2)/Pi^2/(1+4*k+4*k^2)*exp(-1/9*Pi^2*(1+2*k)^2*t)

Solution

>    simplify(subs(T[k](t)=rhs(%),spr));

Sum(24*(-1)^k/Pi^2/(1+4*k+4*k^2)*exp(-1/9*Pi^2*(1+2*k)^2*t)*sin(1/6*Pi*(1+2*k)*x),k = 0 .. infinity)

Solution

>    sol:=subs(k='k',factor(%));

sol := Sum(24*(-1)^k*exp(-1/9*Pi^2*(1+2*k)^2*t)*sin(1/6*Pi*(1+2*k)*x)/Pi^2/(1+2*k)^2,k = 0 .. infinity)

>   

5 Example

Problem

[komec] 98 p.

>    restart;

>    eq:=diff(u(x,t),t)-a^2*diff(u(x,t),x,x)=f(x,t);0<x,x<l,t>0;

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

0 < x, x < l, 0 < t

>    init_c:=u(x,0)=phi(x);

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

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

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

>    a:=4;l:=7;phi:=x->0;f:=(x,t)->2;

a := 4

l := 7

phi := 0

f := 2

Eigenvalues and eigenfunctions

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

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

>    assume(lambda>0);

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

s1 := D(y)(0) = 0

s2 := y(7) = 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 := {_C1*sqrt(lambda) = 0, _C1*sin(7*sqrt(lambda))+_C2*cos(7*sqrt(lambda)) = 0}

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

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

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

chl := sqrt(lambda)*cos(7*sqrt(lambda))

>    chl:=%/sqrt(lambda);   

chl := cos(7*sqrt(lambda))

>    _EnvAllSolutions := true:

>    solve(chl,lambda);

1/196*Pi^2*(1+2*_Z1)^2

>    indets(%) minus {l};

{_Z1}

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

>    tr:=unapply(%,k);

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

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

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

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

>    dsolve(%,y(x));

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

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

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

>    solve(s1,_C1);

0

>    simplify(subs(_C1=%,y(x)));

_C2*cos(1/14*Pi*(1+2*k)*x)

>    select(has,%,x);

cos(1/14*Pi*(1+2*k)*x)

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

1/7*cos(1/14*Pi*(1+2*k)*x)*sqrt(14)

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

tf := proc (k, x) options operator, arrow; 1/7*cos(1/14*Pi*(1+2*k)*x)*sqrt(14) end proc

Eigenvalues and eigenfunctions:

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

1/196*Pi^2*(1+2*k)^2

1/7*cos(1/14*Pi*(1+2*k)*x)*sqrt(14)

>    y:='y':

Solving problem

Solution we find in form:

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

spr := Sum(1/7*T[k](t)*cos(1/14*Pi*(1+2*k)*x)*sqrt(14),k = 0 .. infinity)

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

>    uz:={diff(u(t),t)+a^2*tr(k)*u(t)=int((f(xi))*tf(k,xi),xi=0..l),

>    u(0)=int((phi(xi))*tf(k,xi),xi=0..l)};

uz := {u(0) = 0, diff(u(t),t)+4/49*Pi^2*(1+2*k)^2*u(t) = 4*(-1)^k*14^(1/2)/Pi/(1+2*k)}

>    dsolve(uz,u(t));

u(t) = 49*(-1)^k*14^(1/2)/Pi^3/(1+2*k)^3-49*exp(-4/49*Pi^2*(1+2*k)^2*t)*(-1)^k*14^(1/2)/(Pi^3+6*Pi^3*k+12*Pi^3*k^2+8*Pi^3*k^3)

Solution:

>    simplify(subs(T[k](t)=rhs(%),spr));

Sum(98*(-1)^(1+k)*cos(1/14*Pi*(1+2*k)*x)*(-1+exp(-4/49*Pi^2*(1+2*k)^2*t))/Pi^3/(1+6*k+12*k^2+8*k^3),k = 0 .. infinity)

Solution

>    sol:=subs(k='k',factor(%));

sol := Sum(98*(-1)^(1+k)*cos(1/14*Pi*(1+2*k)*x)*(-1+exp(-4/49*Pi^2*(1+2*k)^2*t))/Pi^3/(1+2*k)^3,k = 0 .. infinity)

Limit of solution,  when   proc (t) options operator, arrow; infinity end proc  

>    eqn:=subs(u(x,t)=y(x),eq); s1,s2;

eqn := diff(y(x),t)-16*diff(y(x),`$`(x,2)) = 2

D(y)(0) = 0, y(7) = 0

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

y(x) = -1/16*x^2+49/16

>    limit(sol,t=infinity);

Sum(-98*(-1)^(1+k)*cos(1/14*Pi*(1+2*k)*x)/Pi^3/(1+2*k)^3,k = 0 .. infinity)

We test this at x=0:

>    subs(x=0,%);

Sum(-98*(-1)^(1+k)*cos(0)/Pi^3/(1+2*k)^3,k = 0 .. infinity)

>    value(%);

98/Pi^3*hypergeom([1, 1/2, 1/2, 1/2],[3/2, 3/2, 3/2],-1)

>    convert(%,StandardFunctions);

98/Pi^3*(-1/2*I*polylog(3,I)+1/2*I*polylog(3,-I))

>    simplify(%);

49/16

>   

6 Example

Problem

>    restart;with(linalg,laplacian):

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

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

0 < r, r < R, 0 < t

>    init_c:=u(r,0)=phi(r);

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

init_c := u(r,0) = phi(r)

bound_c := u(R,t) = 0

>    a:=2;R:=8;phi:=r->64-r^2;

a := 2

R := 8

phi := proc (r) options operator, arrow; 64-r^2 end proc

Solution representation formula

u(r,t) = Sum(A[k]*exp(-(a*mu[k])^2*t/(R^2))*BesselJ(0,mu[k]*r/R),k = 1 .. infinity)

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

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

Zeros of Bessel J(0,x) function

>    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

Note: see ?BesselJZeros

Solution

>    A:=k->2/(R*BesselJ(1,mu[k]))^2*int(r*phi(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*phi(r)*BesselJ(0,mu[k]/R*r),r = 0 .. R) end proc

>    k:='k':

>    sol:=sum(A(k)*exp(-(a*mu[k])^2*t/(R^2))*BesselJ(0,mu[k]*r/R),k=1..Order);

sol := 70.91342475*exp(-.3614491228*t)*BesselJ(0,.3006031948*r)-8.945760312*exp(-1.904453896*t)*BesselJ(0,.6900097638*r)+2.910494148*exp(-4.680437924*t)*BesselJ(0,1.081715989*r)-1.343417750*exp(-8.6900...
sol := 70.91342475*exp(-.3614491228*t)*BesselJ(0,.3006031948*r)-8.945760312*exp(-1.904453896*t)*BesselJ(0,.6900097638*r)+2.910494148*exp(-4.680437924*t)*BesselJ(0,1.081715989*r)-1.343417750*exp(-8.6900...
sol := 70.91342475*exp(-.3614491228*t)*BesselJ(0,.3006031948*r)-8.945760312*exp(-1.904453896*t)*BesselJ(0,.6900097638*r)+2.910494148*exp(-4.680437924*t)*BesselJ(0,1.081715989*r)-1.343417750*exp(-8.6900...
sol := 70.91342475*exp(-.3614491228*t)*BesselJ(0,.3006031948*r)-8.945760312*exp(-1.904453896*t)*BesselJ(0,.6900097638*r)+2.910494148*exp(-4.680437924*t)*BesselJ(0,1.081715989*r)-1.343417750*exp(-8.6900...
sol := 70.91342475*exp(-.3614491228*t)*BesselJ(0,.3006031948*r)-8.945760312*exp(-1.904453896*t)*BesselJ(0,.6900097638*r)+2.910494148*exp(-4.680437924*t)*BesselJ(0,1.081715989*r)-1.343417750*exp(-8.6900...
sol := 70.91342475*exp(-.3614491228*t)*BesselJ(0,.3006031948*r)-8.945760312*exp(-1.904453896*t)*BesselJ(0,.6900097638*r)+2.910494148*exp(-4.680437924*t)*BesselJ(0,1.081715989*r)-1.343417750*exp(-8.6900...

Plots

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

[Maple Plot]

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

[Maple Plot]

3d animation:

>    K := [seq(0.5*k,k=0..15)];

K := [0., .5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5]

>    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=.5`, `t=1.0`, `t=1.5`, `t=2.0`, `t=2.5`, `t=3.0`, `t=3.5`, `t=4.0`, `t=4.5`, `t=5.0`, `t=5.5`, `t=6.0`, `t=6.5`, `t=7.0`, `t=7.5`]
tit := [`t=0.`, `t=.5`, `t=1.0`, `t=1.5`, `t=2.0`, `t=2.5`, `t=3.0`, `t=3.5`, `t=4.0`, `t=4.5`, `t=5.0`, `t=5.5`, `t=6.0`, `t=6.5`, `t=7.0`, `t=7.5`]

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

>   

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