# [Maxima] Find the angle between two tangents by iterating en2

William Porter wmporter.omegapar at comcast.net
Fri Sep 3 11:22:38 CDT 2010

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20100903/b9f8d1bf/attachment.html>
```