4.6.3 2)    heat3-1.mws     Heat Equation     non-homogeneous equation and boundary conditions (Dirichlet-Dirichlet)

>

restart;

>	with(plots):

Warning, the name changecoords has been redefined

>	L:=2;g1:=1;g2:=3;b:=-2;a:=0.5;

Initial Condition:

>	u[0](x):=x*(2*L-x)^2;

Steady State solution:

>	us(x):=b*x^2/2+((g2-g1)/L-b*L/2)*x+g1;

>	plot({us(x),u[0](x)},x=0..L,axes=boxed);

Egenfunctions:

>	X[n]:=sin(n*Pi/L*x);

>	T[n]:=exp(-n^2*Pi^2/a^2/L^2*t);

Fourier Coefficients:

>	d[n]:=2/L*int((u[0](x)-us(x))*X[n],x=0..L);

>	d[n]:=simplify(factor(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},d[n])));

Solution:

>	u(x,t):=us(x)+sum(d[n]*X[n]*T[n],n=1..20):

>	plot3d(u(x,t),x=0..L,t=0..0.3,axes=boxed,projection=0.9,style=wireframe,color=black);

>	animate({u[0](x),u(x,t),us(x)},x=0..L,t=0..0.5,frames=500,axes=boxed);

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

>	u1:=subs(t=0.01,u(x,t)):

>	u2:=subs(t=0.05,u(x,t)):

>	u3:=subs(t=0.2,u(x,t)):

>	plot({us(x),u[0](x),u0,u1,u2,u3},x=0..L,axes=boxed,color=black);

Accuracy of truncated solution

Fo=alpha*t/L^2    Fourier Number (non-dimensional time), alpha = thermodiffusivity

>	T[n]:=exp(-n^2*Pi^2*Fo);

>	UT1(x,t):=us(x)+sum(d[n]*X[n]*T[n],n=1..1):

>	UT4(x,t):=us(x)+sum(d[n]*X[n]*T[n],n=1..4):

>	UT(x,t):=us(x)+sum(d[n]*X[n]*T[n],n=1..100):

>	u0:=subs(Fo=0,UT(x,t)):u10:=subs(Fo=0,UT1(x,t)):u40:=subs(Fo=0,UT4(x,t)):t0:=(0.0)*a^2*L^2;

>	u1:=subs(Fo=0.05,UT(x,t)):u11:=subs(Fo=0.05,UT1(x,t)):u41:=subs(Fo=0.05,UT4(x,t)):t1:=(0.05)*a^2*L^2;

>	u2:=subs(Fo=0.2,UT(x,t)):u12:=subs(Fo=0.2,UT1(x,t)):u42:=subs(Fo=0.2,UT4(x,t)):t2:=(0.2)*a^2*L^2;

>	u3:=subs(Fo=0.30,UT(x,t)):u13:=subs(Fo=0.30,UT1(x,t)):u43:=subs(Fo=0.30,UT4(x,t)):t3:=(0.30)*a^2*L^2;

>	plot({us(x),u[0](x),u0,u1,u2,u3},x=0..L,axes=boxed,color=black,numpoints=200);

>	plot({us(x),u[0](x),u10,u11,u12,u13},x=0..L,axes=boxed,color=black);

>	plot({us(x),u[0](x),u40,u41,u42,u43},x=0..L,axes=boxed,color=black);

>	plot({us(x),u[0](x),u0,u10,u40,u1,u11,u41,u2,u12,u42,u3,u13,u43},x=0..L,axes=boxed,color=black,numpoints=200);

>	animate({us(x),u[0](x),u12,UT(x,t),UT1(x,t),UT4(x,t)},x=0..L,Fo=0..0.5,frames=300,axes=boxed,numpoints=200,color=black);

>

>