Application Center - Maplesoft

# Cauchy problems for wave equations (updated to Maple 7)

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

wave_eqc.mws

Solving Cauchy problems for  wave 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  wave equations.

Cauchy problem   for 1d wave equation (1)

Problem

from M.Kawski

ftp://math.la.asu.edu/pub/kawski/MAPLE/362/dAlembert2hats.mws

 > restart;

 > eq:=diff(u(x,t),t,t)=4*diff(u(x,t),x,x);#x>-infinity,x0;

 > ic:=u(x,0)=u0(x);ict:=D[2](u)(x,0)=u1(x);#x>-infinity,x

 > hat2:=x->1-abs(abs(x)-2);

 > ramp:=x->(x+abs(x))/2;

 > u0:=ramp@hat2;

 > u1:=0;

 > plot(u0,-4..4,-1..2,scaling=constrained,axes=boxed);

 > convert(u0(x),piecewise);

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

 >

Solving(1 method)

We use  d'Alembert's  formula for solution:

 > a:=2:

 > sol1:=(u0(x-a*t)+u0(x+a*t))/2;

 >

Solving(2 method)

 > assume(t>0);

 > with(inttrans):

 > fourier(eq,x,w);

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

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

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

 >

Checking the Solution

sol1=sol2 ???

 > simplify(sol1-sol2);

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

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

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

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

Plots

 > plot3d(sol1,x=-10..10,t=0..5,orientation=[-60,20],axes=framed);

 > with(plots):with(plottools):

```Warning, the name changecoords has been redefined
```

```Warning, the name arrow has been redefined
```

Animation, created by M.Kawski:

 > display([    seq(       display([          textplot([1.4,1.4,cat(`t = `,convert(tt,string))],              'align=RIGHT',font=[TIMES,ROMAN,18]),          plot(subs(t=tt,sol1),x=-6..6,-1..2,              thickness=3)          ]),       tt=[seq(evalf(k/20,3),k=0..100)])    ],insequence=true,    titlefont=[TIMES,BOLD,16],    title=`Animation: d'Alembert's solution`);

 >

Cauchy problem   for 1d wave equation (2)

Problem

 > restart;

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

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

 > u0:=unapply(randpoly([x]),x);

 > u1:=unapply(randpoly([x]),x);

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

Solving(1 method)

We use  d'Alembert's  formula for solution:

 > (u0(x-a*t)+u0(x+a*t))/2+int(u1(xi),xi = x-a*t .. x+a*t)/(2*a);

 > sol1:=simplify(%);

 >

Solving(2 method)

 > L:=f->diff(f,x,x);

 > d:=max(degree(u0(x)),degree(u1(x)));

 > Sum(a^(2*k)*t^(2*k)/(2*k)!*(L@@k)(u0(x))+ a^(2*k)*t^(2*k+1)/(2*k+1)!*(L@@k)(u1(x)),k=0..d);

 > sol2:=simplify(value(%));

 > is(sol1=sol2);

 >

Solving(3 method)

 > with(inttrans):

 > fourier(eq,x,w);

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

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

 > sol3:=invfourier(rhs(%),w,x);

 > is(sol1=sol3);

Checking the Solution

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

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

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

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

 >

Cauchy problem  for 2d wave equation

Problem( 2d wave equation)

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

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

 > ic:=u(x,y,0)=u0(x,y);ict:=D[3](u)(x,y,0)=u1(x,y);

 > u0:=unapply(randpoly([x,y]),x,y);

 > u1:=unapply(randpoly([x,y]),x,y);

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

Solving(1 method)

We use  Poisson formula for solution:

where B={( ):  }.

For computing integrals we use polar coordinates:

 > T:={xi=x+r*cos(phi),eta=y+r*sin(phi)};

 > subs(T,1/(2*Pi*a)* Doubleint(u1(xi,eta)/sqrt((a*t)^2-r^2)*r, phi=0..2*Pi,r=0..a*t)+1/(2*Pi*a)*Diff( Doubleint(u0(xi,eta)/sqrt((a*t)^2-r^2)*r, phi=0..2*Pi,r=0..a*t),t));

 > sol1:=simplify(value(%));

Solving(2 method)

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

 > d:=max(degree(u0(x,y)),degree(u1(x,y)));

 > Sum(a^(2*k)*t^(2*k)/(2*k)!*(L@@k)(u0(x,y))+ a^(2*k)*t^(2*k+1)/(2*k+1)!*(L@@k)(u1(x,y)),k=0..d);

 > sol2:=simplify(value(%));

 > is(sol1=sol2);

 >

Solving(3 method)

 > with(inttrans);

 > fourier(eq,x,w1);

 > fourier(%,y,w2);

 > subs(fourier(fourier(u(x,y,t),x,w1),y,w2)=s(t),%);

 > dsolve({%,s(0)=fourier(fourier(u0(x,y),x,w1),y,w2), D(s)(0)=fourier(fourier(u1(x,y),x,w1),y,w2)},s(t));

 > invfourier(rhs(%),w2,y);

 > sol3:=invfourier(%,w1,x);

 > is(sol1=sol3);

Checking the Solution

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

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

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

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

 >

Cauchy problem  for 3d wave equation

Problem

 > restart;with(linalg):with(student):

```Warning, the protected names norm and trace have been redefined and unprotected
```

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

 > ic:=u(x,y,z,0)=u0(x,y,z);ict:=D[4](u)(x,y,z,0)=u1(x,y,z);

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

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

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

Solving(1 method)

We use  Kirchhoff  formula for solution:

where  S={( ):  }.

 > sol:=Int(u1,S)/(4*Pi*a^2*t)+Diff(Int(u0,S)/t,t)/(4*Pi*a^2);

For computing integrals we use spherical coordinates:

 > T:={xi[1]=x+r*cos(phi)*sin(theta),xi[2]=y+r*sin(phi)*sin(theta), xi[3]=z+r*cos(theta)};

 > linalg[jacobian]([x+r*cos(phi)*sin(theta),y+r*sin(phi)*sin(theta), z+r*cos(theta)],[phi,r,theta]);

 > J:=simplify(det(%));

 > subs(Int(u1,S)=Doubleint(subs(T,u1(xi[1],xi[2],xi[3]))*J, theta=0..Pi,phi=0..2*Pi), Int(u0,S)=Doubleint(subs(T,u0(xi[1],xi[2],xi[3]))*J, theta=0..Pi,phi=0..2*Pi), r=a*t,sol);

 > sol1:=simplify(value(%));

Solving(2 method)

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

 > d:=max(degree(u0(x,y,z)),degree(u1(x,y,z)));

 > Sum(a^(2*k)*t^(2*k)/(2*k)!*(L@@k)(u0(x,y,z))+ a^(2*k)*t^(2*k+1)/(2*k+1)!*(L@@k)(u1(x,y,z)),k=0..d);

 > sol2:=simplify(value(%));

 > is(sol1=sol2);

 >

Solving(3 method)

 > with(inttrans):

```Warning, the name hilbert has been redefined
```

 > fourier(eq,x,w1);

 > fourier(%,y,w2);

 > fourier(%,z,w3);

 > ode:=subs(fourier(fourier(fourier(u(x,y,z,t),x,w1),y,w2),z,w3)=s(t),%);

 > dsolve({ode,s(0)=fourier(fourier(fourier(u0(x,y,z),x,w1),y,w2),z,w3), D(s)(0)=fourier(fourier(fourier(u1(x,y,z),x,w1),y,w2),z,w3)},s(t));

 > invfourier(rhs(%),w3,z):

 > invfourier(%,w2,y):

 > sol3:=invfourier(%,w1,x);

 > is(sol1=sol3);

 >

Checking the Solution

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

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

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

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

 >

Non-homogenous wave equation

Problem

 > restart;with(linalg,laplacian):

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

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

 > ict:=D[4](u)(x,y,z,0)=0;#-infty

The function F(x,y,z,t)  may be polynomial  respect to x,y,z . 

 > randpoly([x,y,z,t,t*sin(t),cos(t),exp(t)]):

 > F:=unapply(%,[x,y,z,t]);

procedure

 > wsol:=proc(u) option `Copyright Aleksas Domarkas, 2002`; local dsol; dsol:=proc(F) local  g,f,d,L; if has(F, t) and not has(F,[x,y,z]) then g:=unapply(F,t); f:=1; elif not has(F,t)  then g:=1;f:=F; else  f:=select(has,F,[x,y,z]); g:=unapply(F/f,t); end if; d:=max(degree(f),1); L:=f->linalg[laplacian](f,[x,y,z]); sum(a^(2*k)/(2*k+1)!*(L@@k)(f)*Int((t-tau)^(2*k+1)*g(tau),tau=0..t),k=0..d);value(%); end proc: if type(u,`+`) then map(dsol,u) else dsol(u) end if: end proc:

solution

 > sol:=wsol(F(x,y,z,t));

 >

test

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

 > simplify(rhs(eq)-lhs(eq)):if %<>0 then combine(%,trig) else % fi;

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

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

 >

Decomposition method

 > restart;with(linalg,laplacian):with(inttrans,fourier,invfourier):

Proble m

From V.S.Vladimirov(ed.), Exercises book on Equations of Mathematical  Physics, Nauka, Moscow, 1982(in Russian) 12.38.7.

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

 > ic:=u(x,y,z,0)=x^2*exp(y+z);#-infty

 > ict:=D[4](u)(x,y,z,0)=sin(x)*exp(y+z);#-infty

I

 > eq1:=eq;

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

 > ict1:=D[4](u)(x,y,z,0)=0;

Solving of I

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

 > simplify(%/rhs(eq1));

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

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

 >

II

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

 > ic2:=u(x,y,z,0)=0;

 > ict2:=D[4](u)(x,y,z,0)=sin(x)*exp(y+z);

Solving of II

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

 > simplify(%/rhs(ict2));

 > dsolve({%,v(0)=0,D(v)(0)=1},v(t));convert(%,trig): combine(%,trig);

 > sol2:=rhs(%)*rhs(ict2);

 >

III

 > eq3:=eq2;

 > ic3:=u(x,y,z,0)=x^2*exp(y+z);

 > ict3:=D[4](u)(x,y,z,0)=0;

Solving of III

 > subs(u(x,y,z,t)=v(x,t)*exp(y+z),eq3);

 > eqn:=simplify(%/exp(y+z));

 > fourier(eqn,x,w);

 > subs(fourier(v(x,t),x,w)=s(t),%);

 > dsolve({%,s(0)=fourier(x^2,x,w),D(s)(0)=0},s(t));

 > subs(s(t)=fourier(v(x,t),x,w),%);

 > invfourier(%,w,x);

 > convert(%,trig): combine(%,trig);

 > sol3:=rhs(%)*exp(y+z);

 >

Solution

 > sol:=sol1+sol2+sol3;

Checking the Solution

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

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

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

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

 >

Radial-symmetric solution

For radial-symmetric solutions of linear elliptic equations  see   radsymell.mws

Problem

 > restart;with(linalg,laplacian):

[vlad] 12.38.9.

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

 > ic:=u(x,y,z,0)=cos((x^2+y^2+z^2)^(1/2));

 > ict:=D[4](u)(x,y,z,0)=cos((x^2+y^2+z^2)^(1/2));

Solving problem

 > eqn:=diff(u(r,t),t,t)-a^2*expand(laplacian(u(r,t),[r,theta,phi],coords=spherical));

Solution we search in form    :

 > subs(u(r,t)=v(r,t)/r,eqn);

 > eq1:=simplify(%*r);

 > ic1:=v(r,0)=r*cos(r);

 > ict1:=D[2](v)(r,0)=r*cos(r);

 > phi:=unapply(rhs(ic1),r);

 > psi:=unapply(rhs(ict1),r);

By  d'Alembert's  formula :

 > v(r,t)=(phi(r+a*t)+phi(r-a*t))/2+Int(psi(xi),xi=r-a*t .. r+a*t)/2/a;

 > value(rhs(%))/r;

 > expand(%);

 > collect(%,{sin(r),cos(r)});

 > applyop(collect,[1,1],%,r):applyop(collect,[2,1],%,r);

Solution

 > sol:=subs(r=sqrt(x^2+y^2+z^2),%);

 >

Checking the Solution

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

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

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

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

 >

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.