m22            heat5dn-2.mws      2-D Heat Equation       DD-NN

>

restart;

>	with(plots):

Warning, the name changecoords has been redefined

>	L:=2;M:=4;alpha:=0.5;

>

f(y):=1;

>	plot(f(y),y=0..M,axes=boxed);

>	g(x,y):=x*(x-L)*y*(y-M);

>	plot3d(g(x,y),x=0..L,y=0..M,axes=boxed);

Steady State Solution:

>	a[0]:=int(f(y),y=0..M)/L/M;

>	a[m]:=2*int(f(y)*cos(m*Pi*y/M),y=0..M)/sinh(m*Pi*L/M)/M;

>	us[m](x,y):=a[m]*sinh(m*Pi*x/M)*cos(m*Pi*y/M):

>	us(x,y):=a[0]*x+sum(us[m](x,y),m=1..2):

>	plot3d(us(x,y),x=0..L,y=0..M,axes=boxed,projection=0.92);

>	plot3d({us(x,y),g(x,y)},x=0..L,y=0..M,axes=boxed,style=wireframe,projection=0.92);

Supplemental Problem:

>	w(x,y):=g(x,y)-us(x,y);

>	A[n,0]:=2*int(int(w(x,y),y=0..M)*sin(n*Pi*x/L),x=0..L)/L/M:

>	A[n,m]:=4*int(int(w(x,y)*cos(m*Pi*y/M),y=0..M)*sin(n*Pi*x/L),x=0..L)/L/M:

>	U[n,0](x,y,t):=A[n,0]*sin(n*Pi*x/L)*exp(-n^2/L^2*Pi^2*t/alpha^2):

>	U[n,m](x,y,t):=A[n,m]*sin(n*Pi*x/L)*cos(m*Pi*y/M)*exp(-(m^2/M^2+n^2/L^2)*Pi^2*t/alpha^2):

>	U(x,y,t):=sum(U[n,0](x,y,t),n=1..20)+sum(sum(U[n,m](x,y,t),m=1..25),n=1..20):

>	U(x,y,0):=subs(t=0,U(x,y,t)):

Solution of IBVP:

>	u(x,y,t):=us(x,y)+U(x,y,t):

>	u(x,y,0):=subs(t=0,u(x,y,t)):

>	animate3d(u(x,y,t),x=0..L,y=0..M,t=0..0.3,frames=100,axes=boxed,projection=0.92);

>

>	animate3d({us(x,y),g(x,y),u(x,y,0),u(x,y,t)},x=0..L,y=0..M,t=0..0.3,frames=100,axes=boxed,style=wireframe,projection=0.92);

>	plot3d({us(x,y),g(x,y),u(x,y,0)},x=0..L,y=0..M,axes=boxed,style=wireframe,projection=0.92);

>	plot3d({g(x,y),u(x,y,0)},x=0..L,y=0..M,axes=boxed,style=wireframe,projection=0.92);