Application Center - Maplesoft
  Submit your work
Application Center Applications Cauchy problems for heat equations (updated to Maple 7) Preview

App Preview:

Cauchy problems for heat equations (updated to Maple 7)

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

Learn about Maple
Download Application

 

heat_eqc.mws

Solving Cauchy problems for  heat equations

by Aleksas Domarkas

Vilnius University, Faculty of Mathematics and Informatics,

Naugarduko 24, Vilnius, Lithuania

aleksas@ieva.mif.vu.lt

NOTE: In this session we solve Cauchy problems for  heat equations

1 Example

Problem

>    restart;

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

>    ic:=u(x,0)=phi(x);#x>-infinity,x<infinity;

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

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

>    a:=1:
phi:=x->piecewise(x<0,0,x<1,x,x<2,2-x,x<3,x-2,x<4,4-x,0);

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

>    'phi(x)'=phi(x);

phi(x) = PIECEWISE([0, x < 0],[x, x < 1],[2-x, x < 2],[x-2, x < 3],[4-x, x < 4],[0, otherwise])

>    plot(phi,-2..6,-1..2,axes=boxed);

[Maple Plot]

Solving problem

>    with(inttrans):

>    phi:=unapply(convert(phi(x),Heaviside),x);

phi := proc (x) options operator, arrow; x*Heaviside(x)-2*x*Heaviside(x-1)+2*Heaviside(x-1)-4*Heaviside(x-2)+2*x*Heaviside(x-2)-2*x*Heaviside(x-3)+6*Heaviside(x-3)-4*Heaviside(x-4)+x*Heaviside(x-4) end...
phi := proc (x) options operator, arrow; x*Heaviside(x)-2*x*Heaviside(x-1)+2*Heaviside(x-1)-4*Heaviside(x-2)+2*x*Heaviside(x-2)-2*x*Heaviside(x-3)+6*Heaviside(x-3)-4*Heaviside(x-4)+x*Heaviside(x-4) end...

>    assume(t>0);

>    fourier(eq,x,w);

diff(fourier(u(x,t),x,w),t)+w^2*fourier(u(x,t),x,w) = 0

>    ode:=subs(fourier(u(x,t),x,w)=s(t),%);

ode := diff(s(t),t)+w^2*s(t) = 0

>    dsolve({ode,s(0)=fourier(phi(x),x,w)},s(t));

s(t) = -(1-2*exp(-I*w)-2*exp(-I*w)*w^2*Pi*Dirac(w)+4*exp(-2*I*w)*w^2*Pi*Dirac(w)+2*exp(-2*I*w)-2*exp(-3*I*w)-6*exp(-3*I*w)*w^2*Pi*Dirac(w)+4*exp(-4*I*w)*w^2*Pi*Dirac(w)+exp(-4*I*w))/w^2*exp(-w^2*t)
s(t) = -(1-2*exp(-I*w)-2*exp(-I*w)*w^2*Pi*Dirac(w)+4*exp(-2*I*w)*w^2*Pi*Dirac(w)+2*exp(-2*I*w)-2*exp(-3*I*w)-6*exp(-3*I*w)*w^2*Pi*Dirac(w)+4*exp(-4*I*w)*w^2*Pi*Dirac(w)+exp(-4*I*w))/w^2*exp(-w^2*t)

>    invfourier(rhs(%),w,x):
sol:=collect(map(simplify,expand(%)),erf);

sol := 1/2*x*erf(1/2*x/t^(1/2))+(3-x)*erf(1/2*(x-3)/t^(1/2))+(x-2)*erf(1/2*(x-2)/t^(1/2))+(1-x)*erf(1/2*(x-1)/t^(1/2))+(1/2*x-2)*erf(1/2*(x-4)/t^(1/2))+1/Pi^(1/2)*t^(1/2)*exp(-1/4*(x-4)^2/t)+1/Pi^(1/2)...
sol := 1/2*x*erf(1/2*x/t^(1/2))+(3-x)*erf(1/2*(x-3)/t^(1/2))+(x-2)*erf(1/2*(x-2)/t^(1/2))+(1-x)*erf(1/2*(x-1)/t^(1/2))+(1/2*x-2)*erf(1/2*(x-4)/t^(1/2))+1/Pi^(1/2)*t^(1/2)*exp(-1/4*(x-4)^2/t)+1/Pi^(1/2)...
sol := 1/2*x*erf(1/2*x/t^(1/2))+(3-x)*erf(1/2*(x-3)/t^(1/2))+(x-2)*erf(1/2*(x-2)/t^(1/2))+(1-x)*erf(1/2*(x-1)/t^(1/2))+(1/2*x-2)*erf(1/2*(x-4)/t^(1/2))+1/Pi^(1/2)*t^(1/2)*exp(-1/4*(x-4)^2/t)+1/Pi^(1/2)...

>    plot3d(sol,x=0..4,t=0..1,orientation=[50,40],numpoints=2000,axes=framed);

[Maple Plot]

Checking the Solution

>    subs(u(x,t)=sol, eq):
simplify(lhs(%)-rhs(%));

0

>    plot(limit(sol,t=0,right), x=0..4);

[Maple Plot]

>   

2 Example

Problem

[vlad] 13.5.8

>    restart;

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

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

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

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

>    a:=1/2;phi:=x->sin(x)*exp(-x^2);

a := 1/2

phi := proc (x) options operator, arrow; sin(x)*exp(-x^2) end proc

>   

Solving by Poisson formula

By Poisson formula:

u(x,t) = Int(phi(xi)*exp(-(x-xi)^2/(4*a^2*t)),xi = -infinity .. infinity)/2/a/sqrt(Pi*t)

>    sol := Int(phi(xi)*exp(-(x-xi)^2/(4*a^2*t)),
xi = -infinity .. infinity)/2/a/sqrt(Pi*t);

sol := Int(sin(xi)*exp(-xi^2)*exp(-(x-xi)^2/t),xi = -infinity .. infinity)/(Pi*t)^(1/2)

>    assume(t>0);

>    convert(%,exp);

Int(-1/2*I*(exp(xi*I)-1/exp(xi*I))*exp(-xi^2)*exp(-(x-xi)^2/t),xi = -infinity .. infinity)/(Pi*t)^(1/2)

>    expand(%);

-1/2*I*Int(1/exp(xi^2)/exp(1/t*x^2)*exp(1/t*x*xi)^2/exp(1/t*xi^2)*exp(xi*I)-1/exp(xi^2)/exp(1/t*x^2)*exp(1/t*x*xi)^2/exp(1/t*xi^2)/exp(xi*I),xi = -infinity .. infinity)/(Pi*t)^(1/2)

>    value(%);

-1/2*I*(exp(-1/4*(4*x^2-4*I*x+t)/(t+1))*t/((t+1)*t)^(1/2)*Pi^(1/2)-exp(-1/4*(4*x^2+4*I*x+t)/(t+1))*t/((t+1)*t)^(1/2)*Pi^(1/2))/(Pi*t)^(1/2)

>    evalc(%);

exp((-x^2-1/4*t)/(t+1))*sin(x/(t+1))*t/(t^2+t)^(1/2)*Pi^(1/2)/(Pi*t)^(1/2)

>    sol1:=simplify(%);

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

Solving using Fourier transform

>    with(inttrans):

>    fourier(eq,x,w);

diff(fourier(u(x,t),x,w),t)+1/4*w^2*fourier(u(x,t),x,w) = 0

>    ode:=subs(fourier(u(x,t),x,w)=s(t),%);

ode := diff(s(t),t)+1/4*w^2*s(t) = 0

>    dsolve({ode,s(0)=fourier(phi(x),x,w)},s(t));

s(t) = (-1/2*I*sqrt(Pi)*exp(-1/4*(w-1)^2)+1/2*I*sqrt(Pi)*exp(-1/4*(w+1)^2))*exp(-1/4*w^2*t)

>    sol2:=simplify(evalc(invfourier(rhs(%),w,x)));

sol2 := exp(-1/4*(4*x^2+t)/(t+1))*sin(x/(t+1))/(t+1)^(1/2)

>    is(sol1=sol2);

true

>   

Checking the Solution

>    u:=unapply(sol1,(x,t)):

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

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

>    u:='u':

0

0

>   

3 Example

Problem

[vlad] 13.6.2

>    restart;with(linalg,laplacian);

[laplacian]

>    eq:=diff(u(x,y,t),t)-a^2*laplacian(u(x,y,t),[x,y])=sin(t)*sin(x)*sin(y);

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

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

ic := u(x,y,0) = 1

>    a:=1;

a := 1

>    f:=(x,y,t)->sin(t)*sin(x)*sin(y);

f := proc (x, y, t) options operator, arrow; sin(t)*sin(x)*sin(y) end proc

>    phi:=(x,y)->1;

phi := 1

Solving by decomposition method

I

>    eq1:=eq;

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

eq1 := diff(u(x,y,t),t)-diff(u(x,y,t),`$`(x,2))-diff(u(x,y,t),`$`(y,2)) = sin(t)*sin(x)*sin(y)

ic1 := u(x,0) = 0

II

>    eq2:=lhs(eq)=0;

>    ic2:=ic;

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

ic2 := u(x,y,0) = 1

solving I

>    subs(u(x,y,t)=v(t)*sin(x)*sin(y),eq1);

diff(v(t)*sin(x)*sin(y),t)-diff(v(t)*sin(x)*sin(y),`$`(x,2))-diff(v(t)*sin(x)*sin(y),`$`(y,2)) = sin(t)*sin(x)*sin(y)

>    dsolve({%,v(0)=0},v(t));

>    s1:=rhs(%)*sin(x)*sin(y);

v(t) = -1/5*cos(t)+2/5*sin(t)+1/5*exp(-2*t)

s1 := (-1/5*cos(t)+2/5*sin(t)+1/5*exp(-2*t))*sin(x)*sin(y)

>   

solving II

>    subs(u(x,y,t)=v(t)*rhs(ic2),eq2);

>    simplify(%/rhs(ic2));

>    dsolve({%,v(0)=1},v(t));

>    s2:=rhs(%)*rhs(ic2);

diff(v(t),t)-diff(v(t),`$`(x,2))-diff(v(t),`$`(y,2)) = 0

diff(v(t),t) = 0

v(t) = 1

s2 := 1

>   

Solution

>    sol:=s1+s2;

sol := (-1/5*cos(t)+2/5*sin(t)+1/5*exp(-2*t))*sin(x)*sin(y)+1

>   

Checking the Solution

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

>    simplify(convert(lhs(eq)-rhs(eq),trig));

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

>    u:='u':

0

0

>   

4 Example

problem

[vlad] 13.6.5.

>    restart;with(linalg,laplacian):with(student,Doubleint):

>    eq:=diff(u(x,y,t),t)-laplacian(u(x,y,t),[x,y])/8=1/8;

>    ic:=u(x,y,0)=exp(-(x-y)^2);a:=sqrt(1/8):

eq := diff(u(x,y,t),t)-1/8*diff(u(x,y,t),`$`(x,2))-1/8*diff(u(x,y,t),`$`(y,2)) = 1/8

ic := u(x,y,0) = exp(-(x-y)^2)

>   

I

>    eq1:=eq;

>    ic1:=u(x,y,0)=0;

eq1 := diff(u(x,y,t),t)-1/8*diff(u(x,y,t),`$`(x,2))-1/8*diff(u(x,y,t),`$`(y,2)) = 1/8

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

II

>    eq2:=lhs(eq)=0;

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

eq2 := diff(u(x,y,t),t)-1/8*diff(u(x,y,t),`$`(x,2))-1/8*diff(u(x,y,t),`$`(y,2)) = 0

ic2 := u(x,0) = exp(-(x-y)^2)

solving I

>    subs(u(x,y,t)=v(t),eq1);

diff(v(t),t)-1/8*diff(v(t),`$`(x,2))-1/8*diff(v(t),`$`(y,2)) = 1/8

>    dsolve({%,v(0)=0},v(t));

>    s1:=rhs(%);

v(t) = 1/8*t

s1 := 1/8*t

>   

solving II

>    phi:=unapply(rhs(ic2),(x,y));

phi := proc (x, y) options operator, arrow; exp(-(x-y)^2) end proc

By Poisson formula

>    1/(sqrt(4*a^2*Pi*t))^2*Doubleint(phi(xi,eta)*exp((-(x-xi)^2-(y-eta)^2)/(4*a^2*t)),xi=-infinity..infinity,eta=-infinity..infinity);

2/Pi/t*Int(Int(exp(-(xi-eta)^2)*exp(2*(-(x-xi)^2-(y-eta)^2)/t),xi = -infinity .. infinity),eta = -infinity .. infinity)

>    assume(t>0);

>    value(%);

exp((-x^2-y^2+2*y*x)/(t+1))/(t+1)^(1/2)

>    s2:=simplify(%);

s2 := exp(-(x-y)^2/(t+1))/(t+1)^(1/2)

>   

Solution

>    sol:=s1+s2;

sol := 1/8*t+exp(-(x-y)^2/(t+1))/(t+1)^(1/2)

>   

Checking the Solution

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

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

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

0

0

>   

5 Example

problem

[vlad] 13.7.4

>    restart;with(linalg,laplacian): with(student,Tripleint):

>    eq:=diff(u(x,y,z,t),t)-laplacian(u(x,y,z,t),[x,y,z])=cos(x-y+z);

>    ic:=u(x,y,z,0)=exp(-(x+y-z)^2);a:=1:

eq := diff(u(x,y,z,t),t)-diff(u(x,y,z,t),`$`(x,2))-diff(u(x,y,z,t),`$`(y,2))-diff(u(x,y,z,t),`$`(z,2)) = cos(x-y+z)

ic := u(x,y,z,0) = exp(-(x+y-z)^2)

>    assume(t>0);

I

>    eq1:=eq;

>    ic1:=u(x,y,0)=0;

eq1 := diff(u(x,y,z,t),t)-diff(u(x,y,z,t),`$`(x,2))-diff(u(x,y,z,t),`$`(y,2))-diff(u(x,y,z,t),`$`(z,2)) = cos(x-y+z)

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

II

>    eq2:=lhs(eq)=0;

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

eq2 := diff(u(x,y,z,t),t)-diff(u(x,y,z,t),`$`(x,2))-diff(u(x,y,z,t),`$`(y,2))-diff(u(x,y,z,t),`$`(z,2)) = 0

ic2 := u(x,0) = exp(-(x+y-z)^2)

solving I

>    subs(u(x,y,z,t)=v(t)*rhs(eq1),eq1);

diff(v(t)*cos(x-y+z),t)-diff(v(t)*cos(x-y+z),`$`(x,2))-diff(v(t)*cos(x-y+z),`$`(y,2))-diff(v(t)*cos(x-y+z),`$`(z,2)) = cos(x-y+z)
diff(v(t)*cos(x-y+z),t)-diff(v(t)*cos(x-y+z),`$`(x,2))-diff(v(t)*cos(x-y+z),`$`(y,2))-diff(v(t)*cos(x-y+z),`$`(z,2)) = cos(x-y+z)

>    dsolve({%,v(0)=0},v(t));

>    s1:=rhs(%)*rhs(eq1);

v(t) = 1/3-1/3*exp(-3*t)

s1 := (1/3-1/3*exp(-3*t))*cos(x-y+z)

>   

solving II

>    phi:=unapply(rhs(ic2),(x,y,z));

phi := proc (x, y, z) options operator, arrow; exp(-(x+y-z)^2) end proc

>    R:=-infinity..infinity;

R := -infinity .. infinity

By Poisson formula

>    1/(sqrt(4*a^2*Pi*t))^3*Tripleint(phi(xi,eta,tau)*exp((-(x-xi)^2-(y-eta)^2-(z-tau)^2)/(4*a^2*t)),xi=R,eta=R,tau=R);

1/8/(Pi*t)^(3/2)*Int(Int(Int(exp(-(xi+eta-tau)^2)*exp(1/4*(-(x-xi)^2-(y-eta)^2-(z-tau)^2)/t),xi = -infinity .. infinity),eta = -infinity .. infinity),tau = -infinity .. infinity)

>    value(%);

1/(Pi*t)^(3/2)*exp((-x^2-y^2-2*y*x-z^2+2*y*z+2*x*z)/(12*t+1))*Pi^(3/2)*t^2*(1+8*t)^(1/2)/((12*t+1)*t*(1+8*t))^(1/2)

>    s2:=simplify(%);

s2 := exp(-(x+y-z)^2/(12*t+1))/(12*t+1)^(1/2)

solution

>    sol:=s1+s2;

sol := (1/3-1/3*exp(-3*t))*cos(x-y+z)+exp(-(x+y-z)^2/(12*t+1))/(12*t+1)^(1/2)

>   

Checking the Solution

>    u:=unapply(sol,(x,y,z,t)):

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

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

>    u:='u':

0

0

>   

6 Example

Problem

[vlad] 13.7.5.

>    restart;with(linalg,laplacian):with(student,Tripleint):

>    eq:=diff(u(x,y,z,t),t)-laplacian(u(x,y,z,t),[x,y,z])=0;#t>0,-infinity<x,y,z<infinity

>    ic:=u(x,y,z,0)=cos(x*y)*sin(z);#-infinity<x,y,z<infinity

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

ic := u(x,y,z,0) = cos(x*y)*sin(z)

>    assume(t>0);

Solving problem

>    phi:=unapply(rhs(ic),(x,y,z));

phi := proc (x, y, z) options operator, arrow; cos(x*y)*sin(z) end proc

>    R:=-infinity..infinity;a:=1:

R := -infinity .. infinity

By Poisson formula

>    1/(sqrt(4*a^2*Pi*t))^3*Tripleint(phi(xi,eta,tau)*exp((-(x-xi)^2-(y-eta)^2-(z-tau)^2)/(4*a^2*t)),xi=R,eta=R,tau=R);

1/8/(Pi*t)^(3/2)*Int(Int(Int(cos(xi*eta)*sin(tau)*exp(1/4*(-(x-xi)^2-(y-eta)^2-(z-tau)^2)/t),xi = -infinity .. infinity),eta = -infinity .. infinity),tau = -infinity .. infinity)

>    convert(%,exp):

>    expand(%):

>    value(%):

Solution

>    u(x,y,z,t)=simplify(evalc(%));

u(x,y,z,t) = exp(-t*(x^2+y^2+1+4*t^2)/(1+4*t^2))*cos(x*y/(1+4*t^2))*sin(z)/(1+4*t^2)^(1/2)

Checking the Solution

>    u:=unapply(rhs(%),(x,y,z,t)):

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

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

>    u:='u':

0

0

>   

7 Example

Problem

>    restart;

>    with(linalg,laplacian):

>    eq:=diff(u(x,y,z,t),t)-a^2*laplacian(u(x,y,z,t),[x,y,z])=f;#-infty<x,y,z<infty,t>0

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

>    ic:=u(x,y,z,0)=phi;#-infty<x,y,z<infty

ic := u(x,y,z,0) = phi

>    phi:=randpoly([x, y, z]);

phi := 79*y+56*x^2*z^2+49*x*y^2*z+63*y^2*z^2+57*x^3*y^2-59*x^2*y^3

>    f:=randpoly([x, y, z]);

f := 77*y*z^2+66*x^3*y+54*x*y*z^2-5*x*z^3+99*y^2*z^2-61*x*z^4

>   

Solving problem

>    L:=f->laplacian(f,[x,y,z]);

L := proc (f) options operator, arrow; laplacian(f,[x, y, z]) end proc

>    d:=max(degree(phi),degree(f),1);

d := 5

>    sol:=int(sum(a^(2*k)*(t-tau)^k/k!*(L@@k)(f),k=0..d),tau=0..t)+
sum(a^(2*k)*t^k/k!*(L@@k)(phi),k=0..d):

Solution:

>    sol:=simplify(%);

sol := 252*t^2*a^2*x*y-15*t^2*a^2*x*z-366*t^2*a^2*x*z^2+342*a^2*t*x*y^2+98*a^2*t*x*z-354*a^2*t*x^2*y+54*x*y*z^2*t+476*a^4*t^2+77*y*z^2*t+66*x^3*y*t-5*x*z^3*t-61*x*z^4*t+132*a^4*t^3-244*a^4*t^3*x+77*t^2...
sol := 252*t^2*a^2*x*y-15*t^2*a^2*x*z-366*t^2*a^2*x*z^2+342*a^2*t*x*y^2+98*a^2*t*x*z-354*a^2*t*x^2*y+54*x*y*z^2*t+476*a^4*t^2+77*y*z^2*t+66*x^3*y*t-5*x*z^3*t-61*x*z^4*t+132*a^4*t^3-244*a^4*t^3*x+77*t^2...
sol := 252*t^2*a^2*x*y-15*t^2*a^2*x*z-366*t^2*a^2*x*z^2+342*a^2*t*x*y^2+98*a^2*t*x*z-354*a^2*t*x^2*y+54*x*y*z^2*t+476*a^4*t^2+77*y*z^2*t+66*x^3*y*t-5*x*z^3*t-61*x*z^4*t+132*a^4*t^3-244*a^4*t^3*x+77*t^2...
sol := 252*t^2*a^2*x*y-15*t^2*a^2*x*z-366*t^2*a^2*x*z^2+342*a^2*t*x*y^2+98*a^2*t*x*z-354*a^2*t*x^2*y+54*x*y*z^2*t+476*a^4*t^2+77*y*z^2*t+66*x^3*y*t-5*x*z^3*t-61*x*z^4*t+132*a^4*t^3-244*a^4*t^3*x+77*t^2...

>   

Checking the Solution

>    u:=unapply(sol,(x,y,z,t)):

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

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

0

0

>   

Other examples

Other examples see in inttrans.mws.

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