MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

RE: Sovling non-linear eq-sys (beginners question)

  • To: mathgroup at smc.vnet.net
  • Subject: [mg71233] RE: [mg71207] Sovling non-linear eq-sys (beginners question)
  • From: "David Park" <djmp at earthlink.net>
  • Date: Sat, 11 Nov 2006 03:38:19 -0500 (EST)

You can solve many conic section problems with the free package from my web
site below. Here is an example of finding a general conic section from 5
points on the curve.

Needs["ConicSections`ConicSections`"]

The general equation for a conic section is

generaleqn = a*x^2 + b*x*y + c*y^2 + d*x + e*y == 1;

We need 5 points for a conic in general position. Here is a test set.

pt1 = {x -> 2, y -> 5};
pt2 = {x -> -1, y -> 4};
pt3 = {x -> -2, y -> -2};
pt4 = {x -> 4, y -> 5};
pt5 = {x -> 6, y -> -7};

We substitute into the general equation and solve the resulting equations.

{generaleqn /. pt1, generaleqn /. pt2, generaleqn /. pt3, generaleqn /. pt4,
    generaleqn /. pt5};
sol = First@Solve[%]
{a -> 1113/38866, b -> 867/38866, c -> 1001/19433,
  d -> -(11013/38866), e -> -(228/19433)}

We then substitute the solution into the general equation to obtain the
specific conic equation. I switched to approximate numbers to keep the
following expressions short.

coniceqn = generaleqn /. sol // N
-0.28335820511501053*x + 0.02863685483455977*x^2 -
   0.01173261977049349*y + 0.022307415221530385*x*y +
   0.051510317501157823*y^2 == 1.

Now, I finally use some routines from the ConicSections package. The first
routine, ParseConic, will parse the equation and return the a semi-axis
size, the eccentricity, and a parametrized curve for the conic. P is an
orthogonal matrix; T is a vector translation; R is either an IdentityMatrix
or a matrix that reflects in the line x == y. They tell how to go from the
conic in standard position to a conic in general position. I don't show all
the output.

{{avalue, eccentricity}, curve[t_], {P, T, R}} = ParseConic[coniceqn];

The curve parametrization that was calculated is

curve[t] =
{5.35468 + 7.89892 Cos[t] + 2.10747 Sin[t], -1.04558 -
    3.21402 Cos[t] + 5.17941 Sin[t]}

ParametricPlot[Evaluate[curve[t]], {t, 0, 2Pi},
    AspectRatio -> Automatic,
    Frame -> True];

We can obtain information about the conic in standard position with the
command...

information = StandardConic[{avalue, eccentricity}]

{conictype -> "Ellipse",
 conicequation ->
   0.013750816629065068*x^2 + 0.031981904691988035*y^2 == 1,
 coniccurve -> {8.527775420556369*Cos[t],
    5.591751174503861*Sin[t]}, coniccurvedomain -> {-Pi, Pi},
 coniccenter -> {0, 0},
 conicfocus -> {{6.438576894460455, 0},
    {-6.438576894460455, 0}},
 conicdirectrix ->
   {x == -11.294880035682084, x == 11.294880035682084},
 conicvertex -> {{8.527775420556369, 0}, {-8.527775420556369, 0}}}

For the corresponding information about the conic in its actual position we
can use...

TransformEllipseRules[P, T, R][information] // Simplify

{conictype -> "Ellipse",
 conicequation ->
   1.*x^2 + x*(-9.89487870619946 + 0.7789757412398916*y) +
     1.7987421383647793*(-4.5214428589005795 + y)*
      (4.293670631128352 + y) == 0,
 coniccurve -> {5.354680644977494 + 7.898924351991414*Cos[t] +
     2.1074687543949118*Sin[t], -1.045581448350521 -
     3.214023600560889*Cos[t] + 5.179406978295339*Sin[t]},
 coniccurvedomain -> {-Pi, Pi},
 coniccenter ->
   {5.354680644977494, -1.045581448350521},
 conicfocus -> {{11.31846713274418, -3.472209387325002},
    {-0.6091058427891927, 1.3810464906239601}},
 conicdirectrix -> {6.413969458840636 + 1.*x ==
     0.4068938322913037*y, 1.*x == 17.97421203377969 +
      0.4068938322913037*y}, conicvertex ->
   {{13.253604996968907, -4.259605048911411},
    {-2.54424370701392, 2.168442152210368}}}

This should work for any type of conic section in any position.


David Park
djmp at earthlink.net
http://home.earthlink.net/~djmp/



From: 17538 at student.hhs.se [mailto:17538 at student.hhs.se]
To: mathgroup at smc.vnet.net

I wonder how I can use Mathematica to solve a system of non-linear
equations.

If I want to find the parameters of an ellipse, given coordinates of
some points on it, then I have a system of four non-linear equations.
Where should I go from here in order to make Mathematica help me with
the substitutions necessary to express each parameter as a function of
the coordinates of the points?

Here follows an example:

The formula of an ellipse (with an axis parallell to the x-axis of the
coordinate system) is:

(x-ox)^2 / s^2 + (y-oy)^2 / b^2 = 1

where
s: half length of small axis.
b: half length of big axis.
ox: x-coordinate of the center of the ellipse.
oy: y-coordinate of the center of the ellipse.
x,y: coordinates of a point on perimeter of the ellipse.

I solve for each of the 4 parameters s, b, ox, oy and insert unique
{x,y}-values in each equation:
s = f(x1,y1,b,ox,oy)
b = f(x2,y2,s,ox,oy)
ox = f(x3,y3,b,s,oy)
oy = f(x4,y4,b,s,ox)

[Of course, there are two functions for each parameter, one describing
the upper half and the other describing the lower half, so to speak.]

Now I should be able to substitute my way to a solution like this:

Inserting b = f(x2,y2,s,ox,oy) in s = f(x1,y1,b,ox,oy) gives:
s = f(x1,y1,f(x2,y2,s,ox,oy),ox,oy)

Solving for s gives:
s = f(x1,y1,x2,y2,ox,oy)

And then continue with inserting oy = f(...) in s' = f(...) above and
again solving for s, then the same with ox so that I arrive at s =
f(x1,y1,x2,y2,x3,y3,x4,y4). Then the procedure is repeated with each of
the other parameters b, ox, oy.

The thing is that it gets quite tedious even for a pretty simple
problem like this. And what about doing it for an ellipse with a tilted
axis (not parallell with the x-axis of the coordinate system)? That
would introduce a new parameter, a new equation in the system and
significantly longer equations too (maybe it is pointless to work with
algebraic solutions in that case anyway).

So, what tricks can Mathematica perform in order to do the work for me?

I'm grateful for any comments or reference that could help me!



  • Prev by Date: RE: 2 dimension Newton Raphson
  • Next by Date: Re: I could do with a list of geometric shape names can anyone help????
  • Previous by thread: Re: Sovling non-linear eq-sys (beginners question)
  • Next by thread: Merge of Matrices 2