MapleHeatSeek6.mws

Heat Seeking Particle

The Billy Johnson Contour Map

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.  The worksheet investigates finding this path.  The temperature on the plate is given by the function

T(x,y) = (-1/18)*(15*x^6+15*y^6-80*x^4-80*y^4+121*x^2+121*y^2-72)

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

>    with(plots): 

Warning, the name changecoords has been redefined

>    x:='x';

>    y:='y';

x := 'x'

y := 'y'

>    f:=(-1/18)*(15*x^6+15*y^6-80*x^4-80*y^4+121*x^2+121*y^2-72);

f := -5/6*x^6-5/6*y^6+40/9*x^4+40/9*y^4-121/18*x^2-121/18*y^2+4

Determining the Direction of Motion (Gradient)

>    fx:=diff(f,x);

fx := -5*x^5+160/9*x^3-121/9*x

>    fy:=diff(f,y);

fy := -5*y^5+160/9*y^3-121/9*y

Approximating the Hottest Point on the Plate

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

HotPt := {x = 0., y = 0.}

>    Tstart:=eval(f,{x=0,y=0});

Tstart := 4

Approximating Some of the Relative Hot and Cold Points on the Plate

>    RelColdPt1:=fsolve({fx=0,fy=0},{x,y},{x=-1.5..-0.5,y=0.5..1.5});

RelColdPt1 := {x = -1.044529939, y = 1.044529939}

>    TC1:=eval(f,RelColdPt1);

TC1 := -2.251940952

>    RelColdPt2:=fsolve({fx=0,fy=0},{x,y},{x=0.5..1.5,y=0.5..1.5});

RelColdPt2 := {y = 1.044529939, x = 1.044529939}

>    TC2:=eval(f,RelColdPt2);

TC2 := -2.251940952

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

RelHotPt1 := {x = 0., y = 1.044529939}

>    TH1:=eval(f,RelHotPt1);

TH1 := .874029524

>    RelHotPt2:=fsolve({fx=0,fy=0},{x,y},{x=-1.7..-1.4,y=1.4..1.7});

RelHotPt2 := {x = -1.569876671, y = 1.569876671}

>    TH2:=eval(f,RelHotPt2);

TH2 := -.9282356e-1

Picking a Starting Point

>    x1:=1.0;

x1 := 1.0

>    y1:=1.95;

y1 := 1.95

Temperature at the Starting Point

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

T1 := -6.22708717

Contour Plot of the Temperature of the Plate

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

[Maple Plot]

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

In this plot the path is rather "jagged" because the step size (timestep = 0.15) is too large.  Also, the object reaches its relatively hot destination before the last iteration of the "for" loop so it winds up "wobbling about" in the relatively hot region.

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

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

gx := -5*P[1]^5+160/9*P[1]^3-121/9*P[1]

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

gy := -5*P[2]^5+160/9*P[2]^3-121/9*P[2]

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

point2d1 := Array(%id = 2579416)

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

route2d1 := Array(%id = 2659984)

>    timestep:=0.15;

timestep := .15

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

point2d1[1] := Vector(%id = 2637408)

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

[Maple Plot]

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

>    x1:=1.0;

x1 := 1.0

>    y1:=-1.0;

y1 := -1.0

Temperature at the Starting Point

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

T1 := -2.222222223

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

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

gx := -5*P[1]^5+160/9*P[1]^3-121/9*P[1]

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

gy := -5*P[2]^5+160/9*P[2]^3-121/9*P[2]

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

point2d2 := Array(%id = 13722704)

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

route2d2 := Array(%id = 14044028)

>    timestep:=0.15;

timestep := .15

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

point2d2[1] := Vector(%id = 13573408)

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

[Maple Plot]

Both Paths

>    display(LevelContour,path2d1,path2d2);

[Maple Plot]

Below we see the temperature function (T) graphed as a surface and also graphed as a surface along with the plane z = 0.  The intersection of the temperature function surface and the plane would give the level curve corresponding to a T value of 0 (which would be an "equal temperature" curve where the temperature was 0 everywhere on the curve.

>    surface:=plot3d(f,x=-2..2,y=-2..2):

>    plane:=plot3d(0,x=-2..2,y=-2..2):

>    display(surface);

[Maple Plot]

>    display(surface,plane);

[Maple Plot]

Below is the graph of the intersection of the surface and plane shown above drawn on the contour map.  This would be the equal temperature curve corresponding to T = 0 described above.

>    EqualTemp0:=implicitplot(15*x^6+15*y^6-80*x^4-80*y^4+121*x^2+121*y^2-72=0,x=-2..2,y=-2..2,color=black,thickness=3,scaling=constrained):

>    display(LevelContour,EqualTemp0);

[Maple Plot]

Below we see more paths constructed using a smaller step size (timestep) and zooming in on the upper right quarter of the contour (heat) map.  Note that the first path is now not so jagged with the smaller step size.

>    x1:=1;

>    y1:=1.95;

x1 := 1

y1 := 1.95

>    LevelContour2:=contourplot(f,x=-0.1..2,y=-0.1..2,contours=32,filled=true,scaling=constrained):

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

gx := -5*P[1]^5+160/9*P[1]^3-121/9*P[1]

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

gy := -5*P[2]^5+160/9*P[2]^3-121/9*P[2]

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

point2d1 := Array(%id = 19529880)

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

route2d1 := Array(%id = 15395672)

>    timestep:=0.03;

timestep := .3e-1

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

point2d1[1] := Vector(%id = 14622112)

>    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(LevelContour2,path2d1);

[Maple Plot]

>    x1:=0.5;

>    y1:=1;

x1 := .5

y1 := 1

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

gx := -5*P[1]^5+160/9*P[1]^3-121/9*P[1]

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

gy := -5*P[2]^5+160/9*P[2]^3-121/9*P[2]

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

point2d1 := Array(%id = 24858400)

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

route2d1 := Array(%id = 13727004)

>    timestep:=0.028;

timestep := .28e-1

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

point2d1[1] := Vector(%id = 13724716)

>    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)]:

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

>    display(LevelContour2,path2d2);

[Maple Plot]

>    x1:=1.2;

>    y1:=1;

x1 := 1.2

y1 := 1

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

gx := -5*P[1]^5+160/9*P[1]^3-121/9*P[1]

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

gy := -5*P[2]^5+160/9*P[2]^3-121/9*P[2]

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

point2d1 := Array(%id = 13947744)

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

route2d1 := Array(%id = 13947384)

>    timestep:=0.028;

timestep := .28e-1

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

point2d1[1] := Vector(%id = 13946704)

>    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)]:

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

>    display(LevelContour2,path2d3);

[Maple Plot]

>    display(LevelContour2,path2d1,path2d2,path2d3);

[Maple Plot]

>    display(LevelContour,path2d1,path2d2,path2d3);

[Maple Plot]

>   

>   

>