Prof. Herod's Recipe for Cheese Cake

Here is a recipe for a cheese cake. If you have to worry about cholesterol, just stay away from this one. If you don't, this one is really great. Fix it a day before you want to impress guests and serve it with the best coffee you have. Wow! Guaranteed a success if you follow the directions.

I cook mine in a ten inch pan. It's about 1.5 inches thick. The color of the real cake has a range from brown, to orange, to yellow.

> with(plots):

> J:=plot3d([r,theta,1.5],r=0..5,theta=-Pi..Pi,coords=cylindrical,

> K:=plot3d([5,theta,z],theta=-Pi..Pi,z=0..1.5,coords=cylindrical,

> display3d({J,K});


Here is the receipe.

Ingredients: Herod's Cheese Cake

2 1/12 lb. cream cheese ( five 8 oz. pkg. )

1 3/4 cup sugar

3 t. flour

1 1/2 t. grated orange rind

1 1/2 t. grated lemon rind

1/4 t. vanilla

5 eggs

2 egg yolks

1/4 c. heavy cream.

Put cheese in mixer bowl. Beat at low speed. Add sugar gradually, then the remainder of ingredients in order . (Eggs should be added one at a time.) When blended and smooth, pour in a lined pan and place in a pre-heated 550 degree oven. Bake 12 - 15 minutes. Reduce heat to 200 degrees. ( Cool oven quickly by leaving the door open and fanning as necessary.) Continue baking for 1 hour. Cool before cutting. The cake is better after refrigeration. Lasts indefinitely ... if you can resist it ... in the refrigerator.

The Mathematical Question

Now, as an applied mathematician, here is your question: Eggs congeal at about 140 degrees. All the ingredients of the cheese cake before cooking are about 46 degrees. When the ingredients are mixed and placed in the hot oven, how long does it take to get the center of the cooking cake to 140 degrees?

The Mathematical Model.

This is a heat diffusion problem. We model it as

[Maple Math] = c [Maple Math] u

side boundary: u(t, 5, [Maple Math] , z) = 550, t > 0, [Maple Math] [Maple Math] [Maple Math] [Maple Math] [Maple Math] , 0 < z < 1.5,

top and bottom boundary: u(t, r, [Maple Math] , 0) = 550 = u(t, r, [Maple Math] , 1.5),

initial condition: u(0, r, [Maple Math] , z) = 46.

That is, we assume the heat diffuses into the cheese cake as modeled by the diffusions equation for a cylinder, that the cylinder is initially at 46 degrees, and the surrounding temperature is 550 degrees.

We ask: at what time is the temperature in the middle 140 degrees? That is, compute t so that

u(t, 0, 0, 3/4 ) = 140.


Recall that when the Laplacian operator is written out for a cylinder the equation is

[Maple Math] + [Maple Math] = [Maple Math] [Maple Math] .

The number c above is the rate of diffusion for heat through the cream cheese/egg mixture. For this problem, since the initial condition and boundary condition are independent of theta, the solution will be independent of theta. Thus, we can rewrite the partial differential equation as

[Maple Math] [Maple Math] = [Maple Math] + [Maple Math] .

Agree that the steady state solution for the equation will be 550 degrees. (As a cook, you must not let the cake achieve this state!) Thus, we solve the problem with homogeneous boundary conditions and add 550. For ease of typing, take a = 5, the radius of the cake, and b = 1.5, the height of the cake.

Perform a separation of variables on u. That is, assume that u(t, r, z) = T(t) R(r) Z(z). We get that

[Maple Math] (r R ') ' Z T + Z '' R T = [Maple Math] T ' R Z.

Dividing by R Z T gives

[Maple Math] (r R ') ' / R + Z '' / Z = [Maple Math] T ' / T.

As usual, we make the argument that each of the terms of the sum on the left side of the above equation must be constant. We have that

(1) [Maple Math] (r R ') ' = [Maple Math] R , R(a) = 0.

(2) Z '' = [Maple Math] Z, Z(0) = Z(b) = 0.


(3) T ' = - c ( [Maple Math] ) T.

We handle these one at a time.

Analysis of equation (1):

We rewrite equation (1) as

(r R ') ' = [Maple Math] r R , R(a) = 0.


[Maple Math] R '' + r R ' + [Maple Math] [Maple Math] R = 0, with R(a) = 0.

We verify that the solution for this last equation is among the Bessel functions.

> R:=r->BesselJ(0,mu*r);


Analysis of equation (2):

Equation (2) is recognized as an equation, with boundary conditions, which defines the sine functions.

> Z:=z->sin(lambda*z);


Analysis of equation (3):

Equation (3) is recognized as an equation for the exponential function.

> T:=t->exp(-c*(mu^2+lambda^2)*t);


Products of solutions of (1), (2), and (3).

Products of solutions for equations (1), (2), and (3) should form a solution for the original equation.

> u:=(t,r,z)->BesselJ(0,mu*r)*sin(lambda*z)*exp(-c*(mu^2+lambda^2)*t);

> 1/r*diff(r*diff(u(t,r,z),r),r)+diff(u(t,r,z),z,z)-1/c*diff(u(t,r,z),t);

> simplify(%);


General solution

We add constants times products of solutions for (1), (2), and (3) to make the general solution of the original equation:

u(t, r, z) = [Maple Math]

We verify that this sum is a solution. For simplicity, take only 9 terms.

> u:=(t,r,z)->sum(sum(A[n,m]*BesselJ(0,mu[m]*r)*sin(lambda[n]*z)*
exp(-c*(mu[m]^2+lambda[n]^2)*t),m = 1 .. 3),n = 1 .. 3);

> 1/r*diff(r*diff(u(t,r,z),r),r)+diff(u(t,r,z),z,z)-1/c*diff(u(t,r,z),t):

> simplify(%);


Computing [Maple Math] 's and [Maple Math] 's.

We identify the eigenvalues [Maple Math] and [Maple Math] ,

For [Maple Math] , we have that sin( [Maple Math] b) = 0, so that [Maple Math] b = n [Maple Math] .

> for n from 1 to 40 do


For [Maple Math] , we have that BesselJ(0, [Maple Math] a) = 0, so that [Maple Math] is the [Maple Math] zero of BesselJ(0,x) divided by a.

> for m from 1 to 40 do


Computing coefficients

We choose the coefficients [Maple Math] so that when t = 0,

46 - 550 = [Maple Math] .

These coefficients will come from the Fourier coefficients:

[Maple Math] = [Maple Math] .

To speed the calculations, we do all the computations for n and then all the computations for m and multiply these. Here are the computations for the sine terms, for the terms involving n.

> a:=5; b:=1.5;

> for n from 1 to 40 do
T[n]:=int(sin(lambda[n]*z),z = 0 .. 3/2)/
int(sin(lambda[n]*z)^2,z = 0 .. 3/2):


Here are the computation for the Bessel terms, for the terms involving m.

> for m from 1 to 30 do
B[m]:=int(BesselJ(0,mu[m]*r)*r,r = 0 .. 5)/
int(BesselJ(0,mu[m]*r)^2*r,r = 0 .. 5):


Here is multiplying these together to get [Maple Math] .

> for n from 1 to 40 do
for m from 1 to 30 do
od: od:


From here on, the program takes a lot of space. Some how, this is not a surprise: if you eat much of this cheese cake, you will take a lot of space! I will indicate what the computations do. You might choose to skip some.

A graph for an approximation of the initial value

Here is the definition for the solution. This is needed in what follows.

> u:=(t,r,z)->sum(sum(A[n,m]*BesselJ(0,mu[m]*r)*sin(lambda[n]*z)* exp(-c*(mu[m]^2+lambda[n]^2)*t),n=1..40),m=1..30);


We are ready to find when the middle of the cake achieves 140. To do this, you need to know c for this cheese cake. It is

> c:=0.0077;

A graphical method for finding when the center is 140 degrees.

A numerical procedure for finding t.

Assignment: Suppose you worked slower and all your ingredients had warmed to 60 degrees when you put the cake in the 550 degree oven. How much sooner would your cheese cake get to 140 in the middle?