Axiheat1.html

Worksheet: Axiheat.mw

Solution to the Axisymmetric Heat Equation

In this worksheet we solve the radial heat equation for u(r,t):

DE: = + r<1

BC: u(1,t)=0

IC: u(r,0)=f(r)

> restart:with(plots):

Cylindrical coordinates

> addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);

I don't like typing BesselJZeros! Maple has a nice alias feature:

> alias(Z=BesselJZeros):

Radial Eigenfunctions

> M:=20; M is the number of terms in the expansion

> R:=seq(BesselJ(0,Z(0,m)*r),m=1..M): R[m] is the m-th eigenfunction

> T:=seq(exp(-(Z(0,m))^2*t),m=1..M): T[m] is the m-th time function

> j:=3:Tmax:=evalf(4/(Z(0,j)^2));Let's graph a solution R[j]*T[j]

> animate3d(R[j]*T[j],r=0..1,theta=0..2*Pi,t=0..Tmax,shading=zhue,style=patchnogrid,coords=z_cylindrical);

Plot_2d

Orthogonal Expansion of the Initial Condition

Now, let's expand the solution for a given initial condition f(r)

> f:=(r)->5*r^3*(1-r);

> plot(f(r),r=0..1,thickness=2);

Plot_2d

Compute the coefficients:

> A:=n->int(r*f(r)*R[n],r=0..1)/int(r*R[n]^2,r=0..1);

proc (n) options operator, arrow; `/`(`*`(int(`*`(r, `*`(f(r), `*`(R[n]))), r = 0 .. 1)), `*`(int(`*`(r, `*`(`^`(R[n], 2))), r = 0 .. 1))) end proc

> M1:=7; M is the number of terms in the expansion

> a:=seq(evalf(A(n)),n=1..M1);

Define the solution

> F:=sum(a[n]*R[n],n=1..M1);

To see if we're on track lets plot the real initial data with its solution approximation (at t=0)

> plot({f(r),F},r=0..1,y=0..(.6),thickness=2);

Plot_2d

Solution to the Heat Equation

> u:=sum(a[n]*T[n]*R[n],n=1..M1);

> animate3d(u,r=0..1,theta=0..2*Pi,t=0..1/2,shading=zhue,style=patchnogrid,coords=z_cylindrical,frames=40);

Plot_2d

>