Application Center - Maplesoft

# Fourier Method for Heat Equation

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

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

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

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

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

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

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

 >

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

 > dsolve(L,y(x));

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

 > sist:={s1,s2};

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

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

 > #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..Pi));

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

Eigenvalues and eigenfunctions:

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

Solving problem

Solution we find in form:

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

 > dsolve(ode,u(t));

Solution

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

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

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

Test

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

 > simplify(%);

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

 > plot(phi,0..Pi);

 >

2 Example

Problem

 > restart;

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

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

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

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

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

 > dsolve(L,y(x));

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

 > sist:={s1,s2};

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

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

 > #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..l));

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

Eigenvalues and eigenfunctions:

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

Solving problem

Solution we find in form:

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

 > dsolve(ode,u(t));

Solution

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

We convert this solution to other form. Index     change to    :

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

 > simplify(%):factor(%);

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

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

 >

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

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

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

 > dsolve(L,y(x));

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

 > sist:={s1,s2};

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

 > chl:=linalg[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..l));

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

Eigenvalues and eigenfunctions:

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

Solving problem

Solution we find in form:

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

 > dsolve(uz,u(t));

Solution

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

or

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

 >

4 Example

Problem

[komec] 93 p.

 > restart;

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

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

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

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

Eigenvalues and eigenfunctions

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

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

 > assume(lambda>0);

 > dsolve(L,y(x));

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

 > sist:={s1,s2};

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

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

 > chl:=%/sqrt(lambda);

 > _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(%,y(x));

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

 > solve(s1,_C2);

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

 > select(has,%,x);

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

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

Eigenvalues and eigenfunctions:

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

Solving problem

Solution we find in form:

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

 > dsolve(uz,u(t));

Solution

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

Solution

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

 >

5 Example

Problem

[komec] 98 p.

 > restart;

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

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

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

 > a:=4;l:=7;phi:=x->0;f:=(x,t)->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);

 > dsolve(L,y(x));

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

 > sist:={s1,s2};

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

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

 > chl:=%/sqrt(lambda);

 > _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(%,y(x));

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

 > solve(s1,_C1);

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

 > select(has,%,x);

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

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

Eigenvalues and eigenfunctions:

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

 > y:='y':

Solving problem

Solution we find in form:

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

 > dsolve(uz,u(t));

Solution:

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

Solution

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

Limit of solution,  when

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

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

 > limit(sol,t=infinity);

We test this at x=0:

 > subs(x=0,%);

 > value(%);

 > convert(%,StandardFunctions);

 > simplify(%);

 >

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

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

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

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

Solution representation formula

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

Zeros of Bessel J(0,x) function

 > 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*phi(r)*BesselJ(0,mu[k]/R*r),r=0..R);

 > k:='k':

 > sol:=sum(A(k)*exp(-(a*mu[k])^2*t/(R^2))*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..10,frames=20);

3d animation:

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

 > m:=nops(%);

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

 > tit:=map(f,K);