Heat136No76.html

Heat Seeking Particle

Adapted from the Maple 8 Getting Started Guide

A particle is placed on a heated plate and begins moving. The motion will always be in the direction of the greatest increase in temperature. The task is to determine the path this "heat seeking" particle will take. This worksheet investigates approximating this path.

The temperature at a point on the plate is given by f. The particle starts at the point (x1,y1). In this example the temperature at any point on the plate is given by the function in Section 13.6#76 in the 8th edition of Calculus by Larson, Hostetler, Edwards. Initially both x and y are allowed to be negative. At the end of the worksheet the questions asked in Section 13.6#76 are answered.

> with(plots):

Warning, the name changecoords has been redefined

> x:='x';

> y:='y';

> f:=400*exp(-(x^2+y)/2);

Determining the Direction of Motion (Gradient)

> fx:=diff(f,x);

> fy:=diff(f,y);

Here is a Contour Plot of the Temperature of the Plate. The "Contour Curves" are curves along which the temperature remains the same. Lighter color indicates warmer temperature so yellow is hot and red is cold.

> contourplot(f,x=-2..2,y=-2..2,contours=8,filled=true);

Picking a Starting Point

> x1:=1;

> y1:=1.5;

Temperature at the Starting Point

> T1:=eval(f,{x=x1,y=y1});

Determining the Path of the Heat Seeking Particle on the Contour Plot

> LevelContour:=contourplot(f,x=-2..2,y=-2..2,contours=8,filled=true):

> gx:=eval(fx,{x=P[1],y=P[2]});

> gy:=eval(fy,{x=P[1],y=P[2]});

> point2d1:=Array(1..45);

> route2d1:=Array(1..45);

> timestep:=0.08;

> point2d1[1]:=<x1,y1>;

> for i from 1 to 44 do
route2d1[i]:=LinearAlgebra[Normalize](eval(<gx,gy>,P=point2d1[i]));
point2d1[i+1]:=eval(<P[1],P[2]>,P=point2d1[i]+timestep*route2d1[i]);
end do:

> listpoints2d1:=[seq(convert(point2d1[i],list),i=1..45)]:

> path2d1:=pointplot(listpoints2d1,style=line,color=blue,thickness=3):

> display(LevelContour,path2d1);

The Path of the Heat Seeking Particle on the Contour Plot with a Different Starting Point

> x1:=-1.6;

> y1:=1;

Temperature at the Starting Point

> T1:=eval(f,{x=x1,y=y1});

Determining the New Path of the Heat Seeking Particle on the Contour Plot

> LevelContour:=contourplot(f,x=-2..2,y=-2..2,contours=8,filled=true):

> gx:=eval(fx,{x=P[1],y=P[2]});

> gy:=eval(fy,{x=P[1],y=P[2]});

> point2d2:=Array(1..45);

> route2d2:=Array(1..45);

> timestep:=0.08;

> point2d2[1]:=<x1,y1>;

> for i from 1 to 44 do
route2d2[i]:=LinearAlgebra[Normalize](eval(<gx,gy>,P=point2d2[i]));
point2d2[i+1]:=eval(<P[1],P[2]>,P=point2d2[i]+timestep*route2d2[i]);
end do:

> listpoints2d2:=[seq(convert(point2d2[i],list),i=1..45)]:

> path2d2:=pointplot(listpoints2d2,style=line,color=blue,thickness=3):

> display(LevelContour,path2d2);

Both Paths

> display(LevelContour,path2d1,path2d2);

Looking at the point (3,5)

> x1:=3;

> y1:=5;

Temperature at the Starting Point

> T1:=eval(f,{x=x1,y=y1});

Determining the New Path of the Heat Seeking Particle on the Contour Plot

> LevelContour2:=contourplot(f,x=0..6,y=0..6,contours=13,filled=true):

> gx:=eval(fx,{x=P[1],y=P[2]});

> gy:=eval(fy,{x=P[1],y=P[2]});

> point2d2:=Array(1..45);

> route2d2:=Array(1..45);

> timestep:=0.15;

> point2d2[1]:=<x1,y1>;

> for i from 1 to 44 do
route2d2[i]:=LinearAlgebra[Normalize](eval(<gx,gy>,P=point2d2[i]));
point2d2[i+1]:=eval(<P[1],P[2]>,P=point2d2[i]+timestep*route2d2[i]);
end do:

> listpoints2d2:=[seq(convert(point2d2[i],list),i=1..45)]:

> path2d2:=pointplot(listpoints2d2,style=line,color=blue,thickness=3):

> display(LevelContour2,path2d2);

> XdirectionMAX:=eval(fx,{x=x1,y=y1});

> YdirectionMAX:=eval(fy,{x=x1,y=y1});

The direction of greatest increase in heat from the point (3,5) would be the direction of the gradient at (3,5) which would be in the direction

of the vector < -6 , -1 >.

The directions of no change in heat on the plate at (3,5) would be directions orthogonal to the gradient at (3,5), i.e., directions tangent to the

level curve at (3,5) which would be in the direction of the vector < 1 , -6 > or the vector < -1 , 6 >.

The pictures below show the plate blown up around the point (3,5) along with the heat-seeking path, direction of greatest increase in heat (in red) and the direction of no change in heat (in green).

> LevelContour3:=contourplot(f,x=2.8..3.2,y=4.8..5.2,contours=13,filled=true,scaling=constrained):

> display(LevelContour3);

> gx:=eval(fx,{x=P[1],y=P[2]});

> gy:=eval(fy,{x=P[1],y=P[2]});

> point2d2:=Array(1..45);

> route2d2:=Array(1..45);

> timestep:=0.0044;

> point2d2[1]:=<x1,y1>;

> for i from 1 to 44 do
route2d2[i]:=LinearAlgebra[Normalize](eval(<gx,gy>,P=point2d2[i]));
point2d2[i+1]:=eval(<P[1],P[2]>,P=point2d2[i]+timestep*route2d2[i]);
end do:

> listpoints2d2:=[seq(convert(point2d2[i],list),i=1..45)]:

> path2d2:=pointplot(listpoints2d2,style=line,color=blue,thickness=3):

> display(LevelContour3,path2d2);

> vector1:=arrow(<3,5>,<-0.2,-0.2/6>,width=[0.045,relative],color=red,thickness=3,scaling=constrained):

> vector2:=arrow(<3,5>,<-0.2/6,0.2>,width=[0.045,relative],color=green,thickness=3,scaling=constrained):

> display(vector1,vector2);

> display(LevelContour3,vector1,vector2);

> display(LevelContour3,path2d2,vector1,vector2);

The surface graphed below is a graph of the temperature distribution function. In the second picture we also see the plane corresponding to the level curve at the (x,y) point (3,5) graphed in green.

> Surface:=plot3d(f,x=2.8..3.2,y=4.8..5.2,axes=boxed,labels=[x,y,T],scaling=constrained):

> display(Surface);

> PlaneLevelCurve:=plot3d(400*exp(-7),x=2.8..3.2,y=4.8..5.2,axes=boxed,labels=[x,y,T],scaling=constrained,color=green):

> display(Surface,PlaneLevelCurve);

Below we see the surface graphed over larger regions of the xy-coordinate plane.

> Surface2:=plot3d(f,x=0..4,y=0..6,axes=boxed,labels=[x,y,T]):

> display(Surface2);

> Surface3:=plot3d(f,x=-5..5,y=-5..5,axes=boxed,labels=[x,y,T]):

> display(Surface3);