Find the angle between two tangents by iterating en2

  William:

I didn't spend much time trying to understand your math, but I did 
notice a few problems.
After I stared at the first 9 input lines, I noticed that I could get 
the same output you got by using:

x^2+d*y*4+d^2*(-4);
e1:rhs(first(solve(x^2+d*y*4+d^2*(-4),y)));
/*Added this line from 'f1()' below*/
e2:subst([x=en1, d=en2],[first(diff(e1,x)), diff(e1,x)]);
en1: 10;
en2: 15;
/*Using %phi conflicts with a constant name maybe?? but didn't change it*/
%phi[1] : first(e2);
/*%phi[1] : first(e2); Removed this line*/
%phi[2] : second(e2);
[(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))];
[(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))];
((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*(-1)));
  %,numer;

I then stared at the function 'f1()' and changed it to use a different 
input parameter name and placed the math statements within the body of 
the for loop:

x^2+d*y*4+d^2*(-4)$
e1:  rhs(first(solve(x^2+d*y*4+d^2*(-4),y)));
en1: 10$   /* set x = 10   */

/*f1(ang) := block([n, en2:5],  Use a different parameter name */
f1(nIn) := block([ang, n:nIn-1, en2:1],    /*add initialization for n, 
and modify en2 value*/
    /*for i:1 thru n while ang # 0.178092938 do ( */
    for i:0 thru n do (
       e2:subst([x=en1, d=en2],[first(diff(e1,x)), diff(e1,x)]),
       %phi[1] : first(e2),
       %phi[2] : second(e2),
       [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))],
       [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))],
       ang : 
((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*(-1))),
       en2 : en2 + 1 /*Need variable increment here that can be less 
than 1.0*/
    ),
    return( float(ang) )
);

The for loop increment is not fine enough to achieve the value you 
sought without more modifications (like a variable increment value that 
gets smaller as you get closer to your target with compensation for 
overrun), but maybe this might get you going in the right direction.

f1(56);
    returns a value of:
0.178099165269

I didn't do any experimentation with the 'newton()' function because I 
thought you might be able to make it work based on what I did with 'f1()'.

Maybe someone knows a nice direct calculation that avoids the iterative 
approach.

I hope this helps,

Paul Bowyer


On 09/03/2010 09:22 AM, William Porter wrote:
>
> Find the angle between two tangents by iterating the distant between 
> the focal point and the vertex of the parabola, d=en2. I have the 
> following parabolic equation x^2+d*y*4+d^2*(-4) and want to find d=en2 
> when the angle between two tangents is 0.1674317785. I know manually 
> that...
>
> (%i1) x^2+d*y*4+d^2*(-4)$
>
> (%i2) e1:  rhs(first(solve(x^2+d*y*4+d^2*(-4),y)));
>
> (%i3) en1: 10$
>
> (%i4) en2: 5$
>
> (%i5) %phi[1] : first(e2)$
>
> (%i6) %phi[1] : first(e2)$
>
> (%i7) %phi[2] : second(e2)$
>
> (%i8) [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))]$
>
>           [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))]$
>
>           
> ((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*(-1)));
>
> (%i9) %,numer;
>
> This gives me the angle (%i8) between the two tangents as 
> 1.5707963268. Now if I keep incrementing d=en2 by 10 it will get me 
> very close, i.e.,
>
> en2: 15 gives me ang: 0.6435011088
>
> en2: 25 gives me ang: 0.3947911197
>
> etc. en2...
>
>  en2: 55 gives me ang: 0.1813197744
>
> Now if I switch from a for loop method to say the Newton method then 
> this should get closer to a value with say 10 decimal places. Problem 
> is I do not know how to set this up. I think it is something like the 
> following, but I am not sure.
>
> x^2+d*y*4+d^2*(-4)$
>
> e1:  rhs(first(solve(x^2+d*y*4+d^2*(-4),y)));
>
> en1: 10$                        /* set x = 10   */
>
> f1(ang) := block([n, en2:5],
>
>               e2: subst([x=en1, d=en2],[first(diff(e1,x)), diff(e1,x)]),
>
>               %phi[1] : first(e2),
>
>               %phi[2] : second(e2),
>
>               ang : [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))],
>
>                        [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))],
>
>                        
> ((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*(-1))),
>
>               for i:1 thru n while ang # 0.178092938 do en2 : en2+i,
>
>               return (en2));         /*Should return en2:55. I do not 
> understand what's going on here. */
>
> newton(exp,var,x0,eps):=
>
>    block([en2:5,e2,xn,s,numer:true],
>
>      s:diff(exp,var),
>
>      xn:x0,
>
>    loop,
>
>      if abs(subst(xn,var,exp))<eps then return(xn),
>
>      xn:xn-subst(xn,var,exp)/subst(xn,var,s),
>
>      e2: subst([x=en1, d=en2],[first(diff(e1,x)),diff(e1,x)]),
>
>      %phi[1] : first(e2),
>
>      %phi[2] : second(e2),
>
>      ang : [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))],
>
>            [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))],
>
>            
> ((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*(-1))),
>
>      go(loop))$
>
> h:expand(e1);   /*I am not sure what goes in here? In either case it 
> does not work.  */
>
> newton(h,x,3.5,1.0e-10);         /*Last value should be the angle 
> 0.1674317785 */
>
> Can someone help me with the algorithms?
>
> William Porter
>
> (425) 343-7711
>
> wporter at omegapar.com
>
>
>
> __________ Information from ESET NOD32 Antivirus, version of virus 
> signature database 5420 (20100903) __________
>
> The message was checked by ESET NOD32 Antivirus.
>
> http://www.eset.com
>
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima