Problem 8:

A square plate [-1, 1]  x [-1, 1]  is at temperature u = 0 .  At time t = 0  the temperature is increased to u = 5  along one of the four sides while being held at u = 0  along the other three sides, and heat then flows into the plate according to u[t] = Delta*u .  When does the temperature reach u = 1  at the center of the plate?

Solution:

This is a standard heat-equation problem with homogeneous boundary conditions.  It will be slightly more convenient to take the plate as [0,2] x [0,2], and the side with u = 5  will be taken as x = 2 .  The solution is obtained as the sum of a steady-state solution v (x, y) = sum(a[j]*sinh(j*Pi*x/2)*sin(j*Pi*y/2),j = 1 .. infinity)  with a[j] = int(5/sinh(j*Pi)*sin(j*Pi*y/2),y = 0 .. 2)  so that v(2,y) = 5 , and a solution of the equation with homogeneous boundary conditions, w(x,y,t) = sum(sum(b[i,j]*sin(i*Pi*x/2)*sin(j*Pi*y/2)*exp(-Pi^2*(i^2+j^2)/4*t),i = 1 .. infinity),j = 1 .. infinity)  with b[i,j] = -a[j]*int(sin(i*Pi*x/2)*sinh(j*Pi*x/2),x = 0 .. 2)  so that w(x,y,0) = -v(x,y) .

>    aj:= int(5/sinh(j*Pi)*sin(j*Pi*y/2),y = 0 .. 2) assuming j::posint;

aj := -10*((-1)^j-1)/sinh(j*Pi)/j/Pi

>    bij:= -aj*int(sin(i*Pi*x/2)*sinh(j*Pi*x/2),x = 0 .. 2) assuming i::posint,j::posint;

bij := 5*((-1)^j-1)/sinh(j*Pi)/j/Pi^2*(j*((-1)^i)^2*I+exp(2*j*PPi)*((-1)^i)^2*I+i*exp(2*j*Pi)*((-1)^i)^2-I*j-i-I*exp(2*j*Pi)*j)/(j*I+i)/(j*I-i)*exp(-j*Pi)/((-1)^i)
bij := 5*((-1)^j-1)/sinh(j*Pi)/j/Pi^2*(j*((-1)^i)^2*I+exp(2*j*Pi)*i-i*((-1)^i)^2+j*exp(2*j*Pi)*((-1)^i)^2*I+i*exp(2*j*Pi)*((-1)^i)^2-I*j-i-I*exp(2*j*Pi)*j)/(j*I+i)/(j*I-i)*exp(-j*Pi)/((-1)^i)
bij := 5*((-1)^j-1)/sinh(j*Pi)/j/Pi^2*(j*((-1)^i)^2*I+exp(2*j*Pi)*i-i*((-1)^i)^2+j*exp(2*j*Pi)*((-1)^i)^2*I+i*exp(2*j*Pi)*((-1)^i)^2-I*j-i-I*exp(2*j*Pi)*j)/(j*I+i)/(j*I-i)*exp(-j*Pi)/((-1)^i)

>    bij:= simplify(evalc(convert(%,exp))) assuming i::posint,j::posint;

bij := -20*((-1)^(i+j)-(-1)^i)*i/(i^2+j^2)/Pi^2/j

We want to evaluate this at x=1, y=1.  By symmetry it is obvious that v(1,1) = 5/4.

As for w(x,y,t) , the factor exp(-Pi^2*(i^2+j^2)/4*t) < 10^(-31)  for 1/4 < t  (some preliminary investigation assured us that the answer we want should be somewhere between about 1/4 and 1) if 116 < i^2+j^2 .  So this approximation to u(1,1,t)  should be accurate to approximately 30 decimals when 1/4 < t  :

>    U:= 5/4 + add(add(bij*sin(i*Pi*1/2)*sin(j*Pi*1/2)*exp(-Pi^2*(i^2+j^2)/4*t),j=1..floor(sqrt(116-i^2))),i=1..10);

U := 5/4-20/Pi^2*exp(-1/2*Pi^2*t)+40/3/Pi^2*exp(-5/2*Pi^2*t)-8/Pi^2*exp(-13/2*Pi^2*t)+172/35/Pi^2*exp(-25/2*Pi^2*t)-40/9*1/Pi^2*exp(-41/2*Pi^2*t)-20/9*1/Pi^2*exp(-9/2*Pi^2*t)+8/3/Pi^2*exp(-17/2*Pi^2*t)...
U := 5/4-20/Pi^2*exp(-1/2*Pi^2*t)+40/3/Pi^2*exp(-5/2*Pi^2*t)-8/Pi^2*exp(-13/2*Pi^2*t)+172/35/Pi^2*exp(-25/2*Pi^2*t)-40/9*1/Pi^2*exp(-41/2*Pi^2*t)-20/9*1/Pi^2*exp(-9/2*Pi^2*t)+8/3/Pi^2*exp(-17/2*Pi^2*t)...
U := 5/4-20/Pi^2*exp(-1/2*Pi^2*t)+40/3/Pi^2*exp(-5/2*Pi^2*t)-8/Pi^2*exp(-13/2*Pi^2*t)+172/35/Pi^2*exp(-25/2*Pi^2*t)-40/9*1/Pi^2*exp(-41/2*Pi^2*t)-20/9*1/Pi^2*exp(-9/2*Pi^2*t)+8/3/Pi^2*exp(-17/2*Pi^2*t)...
U := 5/4-20/Pi^2*exp(-1/2*Pi^2*t)+40/3/Pi^2*exp(-5/2*Pi^2*t)-8/Pi^2*exp(-13/2*Pi^2*t)+172/35/Pi^2*exp(-25/2*Pi^2*t)-40/9*1/Pi^2*exp(-41/2*Pi^2*t)-20/9*1/Pi^2*exp(-9/2*Pi^2*t)+8/3/Pi^2*exp(-17/2*Pi^2*t)...

Now to solve u(1,1,t) = 1 :

>    Digits:= 40: fsolve(U=1,t=1/4 .. 1);

.4240113870336883637974336685932564512478

This actually has 39 correct digits.