graph of function with metapost

Question

I have to create graph of function 3x^4 − 11x^3 + 12x^2 − 4x + 1 ; interval <0,1;1> .

I have just func.mp

 input gjktc.mp;
u:=1mm;

path p[];
beginfig(1);

draw (0,0)--(10u,10u);
draw (0,10u)--(10u,20u);

for i=1 upto 72:
draw (0u,0u+i*2u)--(10u,10u+i*2u);
endfor;

for i=1 upto 72:
draw (10u,10u+i*2u)--(20u,0u+i*2u);
endfor;

%penc 0.2pt;
%p2=(0,0) for i=1 upto 72:
%..(i*0.087*u,-sind(i*5)*u)
%endfor;

%draw p2 rotated 15 withcolor red;

endfig;

end;

I am sure, that is not well and not enought. Please help me. Graph has to look like this:

Answer 1

As I couldn't run you code, here are some general thoughts on drawing the graph of a function like 3x^4 − 11x^3 + 12x^2 − 4x + 1.

I would start by defining a macro that returns the value of this function:

vardef func (expr X) = 
  3X*X*X*X - 11X*X*X + 12X*X - 4X + 1
enddef;

Since MP can natively draw only cubics in one step, one needs to choose a subdivision of the interval into parts. Normally about 10 subintervals is plenty. Then, I prefer to compute the points on the graph and store them in an array. Taking the interval from 0 to 1 and 10 subintervals:

pair Z[]; numeric N;
N=10;
for j=0 upto N:
  Z[j] := ( j/N, func (j/N) );
endfor

Now create and draw the path:

path P;
P := Z0 for j=1 upto N: ..Z[j] endfor;
u := 10cm;
beginfig(1)
  draw P scaled u;
endfig

You will find that for some functions, the resulting path, though it passes through the required points, can get a little wobbly in between. That doesn't happen in this case, but when it does, it helps a lot if the derivative is used to supply direction information.

If you want shaded rectangles that match the graph as in your figure, the following will produce them. Create these before drawing the curve of the function, so it is visible on top of them. It is useful to first define a command that returns a rectangle path when the corners are given

vardef rect (expr LL, UR) = 
  LL -- (xpart UR, ypart LL) -- UR 
     -- (xpart LL, ypart UR) -- cycle
enddef;
beginfig(2)
  for j=1 upto N: 
    fill rect ((xpart Z[j-1],0), Z[j]) scaled u withcolor .8white;
    draw rect ((xpart Z[j-1],0), Z[j]) scaled u withcolor black;
  endfor
  draw P scaled u;
endfig
end

If you insist on shading with slanted lines, ask another question about that.

Answer 2

The graph of an arbitrary function cannot be rendered accurately using Bézier curves. For MetaPost, there is a bpolynomial package for accurately rendering polynomials of up to order 3. But since the polynomial considered here is of order 4, that package is not of much help. When drawing graphs of higher-order polynomials or non-polynomial functions, some kind of piecewise interpolation is required.

The naive approach is to construct a path by connecting some sample points using MetaPost's .. operator. This can lead to sufficient results. But even when interpolating well-behaving functions, the shape of such a path often looks strange due to arbitrary slopes and curvatures at the end-points.

Daniel Luecking's splines package provides a few macros for piecewise spline interpolation. When called with a false first parameter, the fcnspline macro computes a Bézier path such that there is no curvature at the end-points (second derivative equates to zero). This leads to a bit more ‘regular’ shapes. The package has some more features, e.g., sample points don't need to be distributed evenly or you can add your own equations to influence path shape. But the manual is really terse on all that ...

The following example shows two paths, calculated using the naive approach (violet) and calculated via fcnspline(false) (black). Note, how parentheses have been used in macro f, since * and ** operators have the same precedence in MetaPost. Drawing the hatched bars is left as an exercise. (There should be more than one hatching related MetaPost package.)

outputtemplate := "%j-%c.mps";
prologues := 2;

input splines

% Evaluate a function.
def f(expr x) =
  3(x**4) - 11(x**3) + 12(x**2) - 4x + 1
enddef;

% Construct naive Bézier path.
% Parameters are left and right interval boundaries and the number
% of sample points.
def naive_path(expr l,r,k)=
    (l, f(l))
    for i = 1 upto k-1:
      .. ((i/(k-1))[l,r], f((i/(k-1))[l,r]))
    endfor
enddef;

% Construct list of sample points for use with splines package.
% Parameters are left and right interval boundaries and the number
% of sample points.
def list(expr l,r,k)=
    (l, f(l))
    for i = 1 upto k-1:
      , ((i/(k-1))[l,r], f((i/(k-1))[l,r]))
    endfor
enddef;

% Show sample points.  Argument is a list of points.
def dodots (text t) =
  for _loc = t:
    drawdot (u*_loc) withcolor red;
  endfor
enddef;

% Output scale.
numeric u; u:=1cm;
% Adjust arrow length.
ahlength := 2bp;

path xaxis, yaxis;
xaxis := (-.1u,0)--(1.1u,0);
yaxis := (0,-.1u)--(0,1.1u);

% interval and sample points
numeric a,b,k;
a := 0.1;
b := 1;
k := 9;

beginfig(1);
  % Draw coordinate system.
  drawarrow xaxis;
  drawarrow yaxis;
  % Draw graph by naive spline construction.
  draw (naive_path(a,b,k)) scaled u withcolor .75(red+blue);
  % Draw graph using splines package.
  draw fcnspline(false)(list(a,b,k)) scaled u;
  % Draw sample points.
  dodots(list(a,b,k));
endfig;

end

Edit - To elaborate a bit on the example

For the function considered here, the angles of gradient at x = 0.1 and x = 1 are ca. -62 and -45 degrees, resp. With what I called the naive approach above we get a path with angles of -85 and -52 degrees, with macro fcnspline(false) we get angles of -53 and -39 degrees. While both solutions are quite poor at resembling the behaviour at the end points, the splines package's solution is clearly better at x = 0.1.

That doesn't mean the fcnspline macro always generates a path with superior shape. But the zero curvature condition tends to avoid surprising shapes. Other advantages of the splines package are, it provides ready-made solutions for spline interpolation, and with the fcnspline macro there's no need to refer to the derivative of a function. Only the actual function needs to be evaluated - with acceptable results.

One trick to improve path shape (there are many) - again without referring to the derivative - is to construct a path that covers a larger interval than needed, then cutting it to the considered interval via the subpath of operator. That doesn't work for all functions, e.g., near pole points, but you get the idea.

In the example below, I've increased the interval by Δx on both sides, adding two sample points, and finally cut off the first and last path segments via subpath (1,k) of <path>. For both approaches, paths have been drawn over those from the last example in a lighter shade. Angles of gradient for all four paths are written to the .log file.

outputtemplate := "%j-%c.mps";
prologues := 2;

input splines

% Evaluate a function.
def f(expr x) =
  3(x**4) - 11(x**3) + 12(x**2) - 4x + 1
enddef;

% Construct naive Bézier path.
% Parameters are left and right interval boundaries and the number
% of sample points.
def naive_path(expr l,r,k)=
    (l, f(l))
    for i = 1 upto k-1:
      .. ((i/(k-1))[l,r], f((i/(k-1))[l,r]))
    endfor
enddef;

% Construct list of sample points for use with splines package.
% Parameters are left and right interval boundaries and the number
% of sample points.
def list(expr l,r,k)=
    (l, f(l))
    for i = 1 upto k-1:
      , ((i/(k-1))[l,r], f((i/(k-1))[l,r]))
    endfor
enddef;

% Show sample points.  Argument is a list of points.
def dodots (text t) =
  for _loc = t:
    drawdot (u*_loc) withcolor .75(red+green);
  endfor
enddef;

% Return a string containing angles of gradient at end points
% for a given path.
def anglesofgradient(expr p) =
  decimal angle((postcontrol 0 of p) - (point 0 of p))
  & " and "
  & decimal angle((point length p of p) - (precontrol length p of p))
enddef;


% Output scale.
numeric u; u := 1cm;
% Adjust arrow length.
ahlength := 2bp;

% interval and sample points
numeric a,b,k,delta;
a := 0.1;
b := 1;
k := 9;
delta := (b-a)/(k-1);

beginfig(1);
  path pn, pntwo, ps, pstwo;
  % Graph using naive spline construction.
  pn := naive_path(a, b, k);
  pntwo := subpath (1,k) of (naive_path(a-delta, b+delta, k+2));
  % Graph using splines package.
  ps := fcnspline(false)(list(a, b, k));
  pstwo := subpath (1,k) of fcnspline(false)(list(a-delta, b+delta, k+2));
  % Draw graphs.
  draw pn scaled u withcolor .75(red+blue) withpen pencircle scaled 1bp;
  draw ps scaled u withcolor black withpen pencircle scaled .75bp;
  draw pntwo scaled u withcolor .5[.75(red+blue),white] withpen pencircle scaled .4bp;
  draw pstwo scaled u withcolor .5[black,white] withpen pencircle scaled .2bp;
  % Draw sample points.
  dodots(list(a,b,k));
  % Show angles.
  show "angles at end points";
  show "pn    (dark violet) : " & anglesofgradient(pn);
  show "ps    (black)       : " & anglesofgradient(ps);
  show "pntwo (light violet): " & anglesofgradient(pntwo);
  show "pstwo (gray)        : " & anglesofgradient(pstwo);
endfig;

end

These are the angles of gradient at the end points:

>> "pn    (dark violet) : -84.5836 and -52.30644"
>> "ps    (black)       : -53.27785 and -38.72513"
>> "pntwo (light violet): -64.68489 and -44.6875"
>> "pstwo (gray)        : -64.88264 and -46.19849"