Plotting relations

[BigMoosie]

I am trying to write a piece of software that will plot a relation (or implicit function). Is there a know easy solution to doing this or will I have to use approximation methods? At the moment I have it trying every number from -100 to 100 for both x and for y and shading the pixel according to the absolute error, it is very slow and as you can see from these casini curves the line width varies.

 

Please let me know if you have any suggestions or know of any tutorial that might help.

[woelen]

I am trying to write a piece of software that will plot a relation (or implicit function). Is there a know easy solution to doing this or will I have to use approximation methods? At the moment I have it trying every number from -100 to 100 for both x and for y and shading the pixel according to the absolute error' date=' it is very slow and as you can see from these casini curves the line width varies.

 

Please let me know if you have any suggestions or know of any tutorial that might help.

I've once written such a program too. It required one point to be entered on the curve by hand and subsequent points are computed numerically.

 

A curve is described by an implicit relation

 

f(x,y) = 0

 

So, you have to compute one pair (x0, y0), which satisfies this equation. That is your starting point.

 

Next, taking differentials of this equation results in the following:

 

∂f/∂x * dx + ∂f/∂y * dy = 0.

 

Computation of ∂f/∂x and ∂y/∂x in the point (x0, y0) can be done numerically for a generic function f, but if f has a special known form, it is best to compute the partial derivatives by hand and use that in your program.

 

Now, given your point (x0, y0) and the value of the partial derivatives, you can find a vector (dx dy), which is along the curve: any vector

a*(-∂f/∂y ∂f/∂x) satisfies the differential equation, given above.

 

Now you compute x1 = x0 + dx and y1 = y0 + dy and you plot a line segment from (x0 y0) to (x1 y1).

 

Suppose, you want to step along the curve with steps of size h, then you can compute your factor a accordingly. A faster algorithm is to use a 1-norm instead of a 2-norm for determining a. It is not important whether you follow the curve with exactly h or with steps which at most deviate from h by a factor of sqrt(2). So, use the 1-norm |∂f/∂y| + |∂f/∂x| as the length and compute h from that.

 

An even better approach is to make the size of step h dependent on the radius of curvature. A strongly bent curve in point (x0, y0) requires h to be small, a curve, which almost is like a straight line in (x0, y0) allows a large h to be chosen.

 

A formula for the radius of curvature can be found here for a curve, described as x=x(t), y=y(t). This is easy to transform to an implicit form in terms of f(x,y), its partial derivatives ∂f/∂x and ∂f/∂y and its second partial derivatives ∂2f/∂x2, ∂2f/∂y2 and ∂2f/∂x∂y. That is left as an exercise for you.

 

http://mathworld.wolfram.com/RadiusofCurvature.html

 

You see, it is not that difficult. Most of the effort will be to transform the math to a good program, especially if you program in a language like C or C++. Computing first and second order derivatives is cumbersome numerically, due to round-off noise. Another problem may be if you want to be able to define your functions in a flexible way, without the need to recompile your program all the time. For the latter thing I've made a solution. I've developed an easy to use string parsing routine with general variables and all kinds of mathematical functions, which allows you to type expressions like "sin(x)+sin(y)+x*y" on stdin and evaluate this as a function with variables x and y, which you set. Have a look on my web site:

 

http://woelen.scheikunde.net/science/software/misc.html

 

You need to take the module parser.tgz. It is written for UNIX/Linux, but I expect it to compile on Windows without any problems with Visual Studio, Pelles C and DevC++.

 

A way of drifting control is to add -k*h*grad(f*f) to the vector (dx dy). This vector steers your next point (x1, y1) towards the curve if it drifts away. You have this vector almost for nothing, you need to compute f(x,y) for each step and you already have the two partial derivatives ∂f/∂x and ∂f/∂y. This drifting control really is necessary, otherwise your drawn curve slowly moves away from the real curve. This correctionvector -k*h*grad(f*f) is zero, when your point (x,y) is on the curve, but it steers towards the curve if the point is somewhat besides the curve. If you assure that your step h is never larger than a pixel size, then this drift control assures that you will remain very close to the curve and the drift is limited such that you will not see it in your graphs.

 

I implemented this algorithm and it is VERY robust, especially when combined with step-size control. If your curves are reasonably smooth without kinks and sharp bends then it also is very robust even without step-size control. This method is so robust, that even when your initial point is not on the real curve, but a little besides it, then you see that the drawn curve nicely starts in your point (x0 y0) and moves towards the real curve while moving forwards with small steps h.

[BigMoosie]

Woelen, thankyou for your very detailed explaination, it is quite novel and I appreciate it very much. I will do my best to follow your guide.

 

Unfortunately I cannot use your parser as I am not writing this in C++, I am using JavaScript in conjunction with the canvas support found in Firefox 1.5. I am lucky though because JavaScript has a built in function eval() that will execute strings as code, my only hurdle is to stop errors occuring and supporting ^ as exponent.

 

Once again thankyou, I will let you know when I finish.