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 temmperature.  The task is to determine the path this "heat seeking" particle will take.  This worksheet investigates finding this path.

The temperature at a point on the plate is given by f.  The particle starts at the point (x1,y1).

>	with(plots):

>

x:='x';

>

y:='y';

>	f:=100-x^3-2*y^2;

Determining the Direction of Motion (Gradient)

>	fx:=diff(f,x);

>	fy:=diff(f,y);

Approximating the Hottest Point on the Plate

>	HotPt:=fsolve({fx=0,fy=0},{x,y},{x=-1..1,y=-1..1});

>	assign(%);

Picking a Starting Point

>

x1:=x+4;

>

y1:=y+6;

>

x:='x';

>

y:='y';

Temperature at the Starting Point

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

Contour Plot of the Temperature of the Plate

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

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

>	LevelContour:=contourplot(f,x=-5..5,y=-2..8,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.15;

>	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;

>

y1:=4;

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

>	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(LevelContour,path2d2);

Both Paths

>	display(LevelContour,path2d1,path2d2);

Plotting the Analytical Solutions in Green

>	AnalPath2:=plot([1/(3*t-1),4*exp(-4*t),t=0..0.3],x=-5..5,y=-2..8,color=green,thickness=3):

>	AnalPath1:=plot([1/(3*t+1/4),6*exp(-4*t),t=0..4],x=-5..5,y=-2..8,color=green,thickness=3):

>	display(LevelContour,path2d1,path2d2,AnalPath1,AnalPath2);

>