Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
215 lines (167 sloc) 7.81 KB
/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 14.12.0 ] */
/* [wxMaxima: input start ] */
kill(all)$
load("lm_curves.mac");
arrayinfo(group2);
mycurve: GERONO $
/* [wxMaxima: input end ] */
/* [wxMaxima: title start ]
The Rational Parametrization of Plane Curves
[wxMaxima: title end ] */
/* [wxMaxima: comment start ]
This Maxima file parametrizes plane curves from their implicit equations using the method
for those of relative degree two as described in Dave Boyles' article:
"Finding Rational Parametric Curves Of Relative Degree One or Two", Dave Boyles.
MAA College Mathematics Journal November 2010
I have kept Boyles' naming convention for D, P, Q, R, tau and tau0.
The lemniscates of Bernoulli and Gerono have sub-optimal parametrizations. I have derived
cleaner ones using (respectively) a pencil of circles and a pencil of parabolas.
Similarly, the Double Folium and Trifolium parametrizations could be better as these are
implicit equations of relative degree one and need not have this algorithm applied. The
method does still work though using a Q polynomial equal to 1.
Some curves are defined elsewhere with different multiplying constants. These are scaling
factors and don't affect the shape of the curve. Occasionally, I have made constants equal
to one, to clean up the implicit equations. In the case of the Besace, I have left the 'a'
constant in as it is key to the Y scaling and therefore the shape of the curve. If you
plug this in to GeoGebra I suggest using a slider for 'a'.
I am not a Maxima expert. Less still any kind of mathematician! I'm sure there are better
ways to implement the algorithm (especially the P and Q selection step), but I've checked
all ten curves and the parametrizations do all match their associated implicit couterparts.
Usage:
Set one (only) of the curve names to true.
Ctrl-R to run.
Expressions for x(t) and y(t) are the last of the outputs at the end of the file.
-=:logicmonkey:=-
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
eq: group2[mycurve]@eqn $
ex: lhs( eq ) - rhs( eq )$
ex: ratexpand( ex );
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
ext: ev(ex, x=x*t, y=y*t) $
hi_order: hipow( ext, t ) $
rel_order: hi_order - lopow( ext, t )$
if ((rel_order<1) or (rel_order>2)) then (error("Relative order",rel_order,"!={1,2}"))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
extau: ex, y=x*tau;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
extau: ratsimp( extau / x^(hi_order-2) );
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
A: ratcoef( extau, x, 2 );
B: ratcoef( extau, x, 1 );
C: ratcoef( extau, x, 0 );
D: factor( B^2-4*A*C );
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
We want discriminant D = B^2-4AC = P^2*Q (a square * a polynomial Q)
Where Q = alpha + beta*tau + gamma*tau^2
Case by case manual rearrangement of D to determine P and Q:
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
if( mycurve = BERNOULLI )then (P: 2*(1+tau^2), Q: 2*(1 - tau^2) )$
if( mycurve = GERONO )then (P: 2, Q: (1 - tau^2) )$
if( mycurve = BESACE )then (P: 2, Q: (a^2+1 - tau^2) )$
if( mycurve = TRICEVA )then (P: 2*(tau^2-3)*(1+tau^2), Q: 1+tau^2 )$
if( mycurve = WALKER )then (P: tau, Q: (12*tau^2 + 1) )$
if( mycurve = RAMPHOID )then (P: tau, Q: tau*(4- 3*tau) )$
if( mycurve = QUADRIFOLIUM )then (P: 4*tau*(1+tau^2), Q: 1+tau^2 )$
if( mycurve = CARDIOID )then (P: 4*(1+tau^2), Q: 1+tau^2 )$
if( mycurve = CIRCLE )then (P: 2, Q: 1+tau^2 )$
if( mycurve = TACNODE )then (P: 4*tau, Q: 2+tau^2 )$
if( mycurve = VANHOEIJ )then (P: tau-1, Q: (tau+9)*(9*tau+1))$
if( mycurve = PIRIFORM )then (P: 1, Q: 1-4*tau^2 )$
if( mycurve = SCARABAEUS )then (P: 2*b*(tau+1)*(tau-1)*(tau^2+1), Q:(tau^2+1) )$
if( mycurve = CAPRICORNOID )then (P: 2*a*tau, Q: a*b*(tau^2+1) )$
if( mycurve = LINKS )then (P: 4*tau, Q: (2*tau^2+1) )$
P: ratsimp( P );
Q: ratsimp( Q );
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
Examine the coefficients of Q to get tau0:
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
alpha: ratcoef( Q, tau, 0 ) $
beta: ratcoef( Q, tau, 1 ) $
gamma: ratcoef( Q, tau, 2 ) $
if (alpha=0) then
tau0: beta / (t^2 - gamma)
else
tau0: 4*alpha*t/(t^2 - 2*beta*t + beta^2 -4*alpha*gamma)
;
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
R^2 = Q(tau0)
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
(if (alpha = 0) then (
/* Qtau0: (t * beta / (t^2 - gamma))^2 */
R : t * beta / (t^2 - gamma)
) else (
/* Qtau0: alpha * (t^2 - t*beta^2 + 4*alpha*gamma)^2 / (t^2 - 2*beta*t + beta^2 -4*alpha*gamma)^2 */
R : sqrt(alpha) * (t^2 - t*beta^2 + 4*alpha*gamma) / (t^2 - 2*beta*t + beta^2 -4*alpha*gamma)
));
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
P: P, tau=tau0;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
display2d: true$
xt: factor(-B + P*R )/factor( 2*A ), tau=tau0 $
yt: factor( xt*tau0 )$
group2[mycurve]@text;
rat_plane_curve: [xt, yt];
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
Surfaces of Revolution
Rotation of the parametrization (0 < v < 2pi):
[f1, f2] around the y axis: Sy : [f1(u).cos v, f2(u), f1(u).sin v]
[f1, f2] around the x axis: Sx : [f1(u), f2(u).cos v, f2(u).sin v]
Having gone to the effort of generating a rational plane curve parametrization, I prefer the rational surface of rotation.
Cos and sin are replaced with the rational parametrization for a circle [(1-v^2)/(1+v^2), 2v/(1+v^2)]
Rational rotation of the parametrization (-inf < v < +inf):
[f1, f2] around the y axis: Sy : [f1(u).(1-v^2)/(1+v^2), f2(u), f1(u).2v/(1+v^2)]
[f1, f2] around the x axis: Sx : [f1(u), f2(u).(1-v^2)/(1+v^2), f2(u).2v/(1+v^2)]
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
rat_surface_of_rotation_y: [xt * (1-v^2)/(1+v^2), yt, xt * 2*v/(1+v^2)], t=u;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
rat_surface_of_rotation_x: [xt, yt* (1-v^2)/(1+v^2), yt * 2*v/(1+v^2)], t=u;
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
And the transcendental rotations for completeness...
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
tran_surface_of_rotation_y:[xt*cos(v), yt, xt*sin(v)], t=u;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
tran_surface_of_rotation_x:[xt, yt*cos(v), yt*sin(v)], t=u;
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
/*
* Fxyz output in MathMod format
* (not all formulas use a and b)
*/
display2d: false$
my_string: sconcat( "\"Fx\":\[\"", tran_surface_of_rotation_x[1], "\"\],"), a=0.8, b=1.6$
print(my_string)$
my_string: sconcat( "\"Fy\":\[\"", tran_surface_of_rotation_x[2], "\"\],"), a=0.8, b=1.6$
print(my_string)$
my_string: sconcat( "\"Fz\":\[\"", tran_surface_of_rotation_x[3], "\"\],"), a=0.8, b=1.6$
print(my_string)$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
my_string: sconcat( "\"Fx\":\[\"", tran_surface_of_rotation_y[1], "\"\],"), a=0.8, b=1.6$
print(my_string)$
my_string: sconcat( "\"Fy\":\[\"", tran_surface_of_rotation_y[2], "\"\],"), a=0.8, b=1.6$
print(my_string)$
my_string: sconcat( "\"Fz\":\[\"", tran_surface_of_rotation_y[3], "\"\],"), a=0.8, b=1.6$
print(my_string)$
/* [wxMaxima: input end ] */
/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$
You can’t perform that action at this time.