6.83 avoid using brute force numerical methods (26.4.99)

6.83.1 Kevin Ulmes

I have a function that describes temperature throughout space, and a heat seeking particle that starts at a point and flies along the gradient to this function. I had to solve for the coordinates of the point when it was 5 units away from the origin. This problem was one that I had in my Math class, and I ended up solving it as shown below using brute force, which was what we were supposed to do. The function is T=1/2*x^2*y-x*z+x+y, and the initial starting point is (1,0,2). I was wondering if anyone could help me out, in solving this problem exactly, or if there is anyway that I could optimize the current while loop that I have so that it would run faster, and therefore calculate the resultant point more accurately. I am using Maple V r4 to do this:

> with(linalg): 
> with(plots): 
> mag:=proc(vect); RETURN(sqrt(dotprod(vect,vect))): end: 
> T:=1/2*x^2*y-x*z+x+y; 
 
                               2 
                     T := 1/2 x  y - x z + x + y 
 
> t_grad:=convert(grad(T,[x,y,z]),list): 
> initial_point:=[1,0,2]: 
> initial_direc:=subs({x=1,y=5,z=2},t_grad):
 

The following execuation loops determines the location of the particle when it is 5 units away from the origin. The accuracy of this result is determined by setting the err variable in the initial line. Digits are set to 20 to minimize roundoff error during the calculation.

Initilize variables so no errors returned

> Digits:=20: 
> point_current:=initial_point: 
> direc_current:=initial_direc: 
> dist:=mag(point_current): 
> er:=.001: 
> while dist < 5 do 
>   incre:=direc_current*fsolve(mag(scalarmul(convert(direc_current,vector),x)) 
    =er,x=0..infinity): 
>   point_current:=evalm(point_current+incre); 
>   direc_current:=subs({x=point_current[1],y=point_current[2], 
                         z=point_current[3]},t_grad): 
>   dist:=mag(point_current); 
> od: 
> current_point:=convert(point_current,list); 
> distance:=dist; 
> direction:=direc_current;
 

I am look at the current point, and taking the gradient vector, and taking a very small step in that direction. and evaluating the new point to check if it is more than 5 units away from the origin. I believe that I am doing this right, but it takes a very long time to solve, when err < .001 (like > 10 min). And as err > the time increases very rapidly. I’m running the Windows version of MapleVr4, on a PII 400, running Win98 (just in case that is important).

6.83.2 Robert Israel (28.4.99)

Consider the system of differential equations

 
dp/dt = grad(T(p)) 
where p = [x,y,z]. This describes a particle moving so that its velocity is 
equal to the gradient vector. It's more convenient for your purposes, however, 
to write p as a function of r (the distance from 0). By the Chain Rule, 
 
dp/dr = (dp/dt)/(dr/dt) = (grad(T(p)))/(p . grad(T(p))/r) 
      = r grad(T(p))/(p . grad(T(p)))
 

You have the initial condition p(sqrt(5)) = [1,0,2], and you want to find p(5). It’s unlikely that there would be an exact "closed-form" solution, but Maple can do it numerically.

> with(linalg): 
  T:=  1/2*x^2*y-x*z+x+y; 
  G:= grad(T,[x,y,z]); 
  Rs:=map(normal,evalm( sqrt(x^2+y^2+z^2) * G/dotprod([x,y,z],G))); 
  des:= {diff(x(r),r) = subs(x=x(r),y=y(r),z=z(r),Rs[1]), 
         diff(y(r),r) = subs(x=x(r),y=y(r),z=z(r),Rs[2]), 
         diff(z(r),r) = subs(x=x(r),y=y(r),z=z(r),Rs[3])}; 
  ics:= {x(sqrt(5))=1, y(sqrt(5))=0, z(sqrt(5))=2}; 
  F:=dsolve(des union ics, {x(r),y(r),z(r)},numeric); 
  F(5); 
 
[r = 5, y(r) = -1.952448370119240, x(r) = 3.252309163629305, 
 
        z(r) = 3.257365618028734]