yacas 1.3.3-2 / usr / share / yacas / plots.rep / plot2d.ys

//////////////////////////////////////////////////
/// Plot2D --- adaptive two-dimensional plotting
//////////////////////////////////////////////////

/// definitions of backends
Use("plots.rep/backends-2d.ys");

/*
	Plot2D is an interface for various backends (Plot2D'...). It calls
Plot2D'get'data to obtain the list of points and values, and then it calls
Plot2D'<backend> on that data.

	Algorithm for Plot2D'get'data: 
	1) Split the given interval into Div(points+3, 4) subintervals, and split each subinterval into 4 parts.
	2) For each of the parts: evaluate function values and call Plot2D'adaptive
	3) concatenate resulting lists and return
*/

	LocalSymbols(var, func, range, option, options'list, delta, options'hash, c, fc, all'values, dummy)
[

// declaration of Plot2D with variable number of arguments
Function() Plot2D(func);
Function() Plot2D(func, range);
Function() Plot2D(func, range, options, ...);

/// interface routines
1 # Plot2D(_func) <-- ("Plot2D" @ {func, -5:5});
2 # Plot2D(_func, _range) <-- ("Plot2D" @ {func, range, {}});
3 # Plot2D(_func, _range, option_IsFunction) _ (Type(option) = "=" Or Type(option) = "==") <-- ("Plot2D" @ {func, range, {option}});

/// Plot a single function
5 # Plot2D(_func, _range, options'list_IsList)_(Not IsList(func)) <-- ("Plot2D" @ {{func}, range, options'list});

/// Top-level 2D plotting routine:
/// plot several functions sharing the same xrange and other options
4 # Plot2D(func'list_IsList, _range, options'list_IsList) <--
[
	Local(var, func, delta, options'hash, c, fc, all'values, dummy);
	all'values := {};
	options'hash := "ListToHash" @ {options'list};
	// this will be a string - name of independent variable
	options'hash["xname"] := "";
	// this will be a list of strings - printed forms of functions being plotted
	options'hash["yname"] := {};
	// parse range
	If (
		Type(range) = "=" Or Type(range) = "==",	// variable also specified -- ignore for now, store in options
		[
			// store alternative variable name
			options'hash["xname"] := String(range[1]);
			range := range[2];
		]
	);
	If(
		Type(range) = ":",	// simple range
		range := N(Eval({range[1], range[2]}))
	);
	// set default option values
	If(
		options'hash["points"] = Empty,
		options'hash["points"] := 23
	);
	If(
		options'hash["depth"] = Empty,
		options'hash["depth"] := 5
	);
	If(
		options'hash["precision"] = Empty,
		options'hash["precision"] := 0.0001
	);
	If(
		options'hash["output"] = Empty Or IsString(options'hash["output"]) And Plot2D'outputs()[options'hash["output"]] = Empty,
		options'hash["output"] := Plot2D'outputs()["default"]
	);
	// a "filename" parameter is required when using data file
	If(
		options'hash["output"] = "datafile" And options'hash["filename"] = Empty,
		options'hash["filename"] := "output.data"
	);
	// we will divide each subinterval in 4 parts, so divide number of points by 4 now
	options'hash["points"] := N(Eval(Div(options'hash["points"]+3, 4)));
	// in case it is not a simple number but an unevaluated expression
	options'hash["precision"] := N(Eval(options'hash["precision"]));
	// store range in options
	options'hash["xrange"] := {range[1], range[2]};
	// compute the separation between grid points
	delta := N(Eval( (range[2] - range[1]) / (options'hash["points"]) ));
	// check that the input parameters are valid (all numbers)
	Check(IsNumber(range[1]) And IsNumber(range[2]) And IsNumber(options'hash["points"]) And IsNumber(options'hash["precision"]),
		"Plot2D: Error: plotting range '"
		:(ToString()Write(range))
		:"' and/or the number of points '"
		:(ToString()Write(options'hash["points"]))
		:"' and/or precision '"
		:(ToString()Write(options'hash["precision"]))
		:"' is not numeric"
	);
	// loop over functions in the list
	ForEach(func, func'list)
	[
		// obtain name of variable
		var := VarList(func);        // variable name in a one-element list
    	Check(Length(var)<=1,
			"Plot2D: Error: expression is not a function of one variable: "
			:(ToString()Write(func))
		);
		// Allow plotting of constant functions
		If(Length(var)=0, var:={dummy});
		// store variable name if not already done so
		If(
			options'hash["xname"] = "",
			options'hash["xname"] := String(VarList(var)[1])
		);
		// store function name in options
		DestructiveAppend(options'hash["yname"], ToString()Write(func));
		// compute the first point to see if it's okay
		c := range[1];
		fc := N(Eval(Apply({var, func}, {c})));
		Check(IsNumber(fc) Or fc=Infinity Or fc= -Infinity Or fc=Undefined,
			"Plot2D: Error: cannot evaluate function '"
			:(ToString()Write(func))
			:"' at point '"
			:(ToString()Write(c))
			:"' to a number, instead got '"
			:(ToString()Write(fc))
			:"'"
		);
		// compute all other data points
		DestructiveAppend(all'values, Plot2D'get'data(func, var, c, fc, delta, options'hash));
		If(InVerboseMode(), Echo({"Plot2D: using ", Length(all'values[Length(all'values)]), " points for function ", func}), True);
	];
	// call the specified output backend
	Plot2D'outputs()[options'hash["output"]] @ {all'values, options'hash};
];

//HoldArg("Plot2D", range);
//HoldArg("Plot2D", options);
HoldArgNr("Plot2D", 2, 2);
HoldArgNr("Plot2D", 3, 2);
HoldArgNr("Plot2D", 3, 3);

/// this is the middle-level plotting routine; it generates the initial
/// grid, calls the adaptive routine, and gathers data points.
/// func must be just one function (not a list)
Plot2D'get'data(_func, _var, _x'init, _y'init, _delta'x, _options'hash) <--
[
	Local(i, a, fa, b, fb, c, fc, result);
	// initialize list by first points (later will always use Tail() to exclude first points of subintervals)
	result := { {c,fc} := {x'init, y'init} };
	For(i:=0, i<options'hash["points"], i++)
	[
		{a,fa} := {c, fc};	// this is to save time but here a = x'init + i*delta'x
		// build subintervals
		{b, c} := N(Eval({x'init + (i+1/2)*delta'x, x'init + (i+1)*delta'x}));	// this is not computed using "a" to reduce roundoff error
		{fb, fc} := N(Eval(MapSingle({var, func}, {b, c})));
		result := Concat(result,
			Tail(Plot2D'adaptive(func, var, {a,b,c}, {fa, fb, fc}, options'hash["depth"],
				// since we are dividing into "points" subintervals, we need to relax precision
				options'hash["precision"]*options'hash["points"] )));
	];
	result;
];

//////////////////////////////////////////////////
/// Plot2D'adaptive --- core routine to collect data
//////////////////////////////////////////////////
/*
	Plot2D'adaptive returns a list of pairs of coordinates { {x1,y1}, {x2,y2},...}
	All arguments except f() and var must be numbers. var is a one-element list containing the independent variable. The "a,b,c" and "fa, fb, fc" arguments are values of the function that are already computed -- we don't want to recompute them once more.
	See documentation (Algorithms.chapt.txt) for the description of the algorithm.
*/

Plot2D'adaptive(_func, _var, {_a,_b,_c}, {_fa, _fb, _fc}, _depth, _epsilon) <--
[
	Local(a1, b1, fa1, fb1);

	a1 := N(Eval((a+b)/2));
	b1 := N(Eval((b+c)/2));
	{fa1, fb1} := N(Eval(MapSingle({var, func}, {a1, b1})));
	If(
		depth<=0 Or
		(
		  // condition for the values not to oscillate too rapidly
		  sign'change(fa, fa1, fb) + sign'change(fa1, fb, fb1) + sign'change(fb, fb1, fc) <= 2
		  And
		  // condition for the values not to change too rapidly
		  N(Eval(Abs( (fa-5*fa1+9*fb-7*fb1+2*fc)/24 ) ))	// this is the Simpson quadrature for the (fb,fb1) subinterval (using points b,b1,c), subtracted from the 4-point Newton-Cotes quadrature for the (fb,fb1) subinterval (using points a, a1, b, b1)
		  	<= N(Eval( epsilon*(	// the expression here will be nonnegative because we subtract the minimum value
		    //(fa+fc+2*fb+4*(fa1+fb1))/12	// this is 1/4 of the Simpson quadrature on the whole interval
			(5*fb+8*fb1-fc)/12	// this is the Simpson quadrature for the (fb,f1) subinterval
			- Min({fa,fa1,fb,fb1,fc}) ) ) )
		),
		// okay, do not refine any more
		{{a,fa}, {a1,fa1}, {b,fb}, {b1,fb1}, {c,fc}},
		// not okay, need to refine more
		Concat(
			// recursive call on two halves of the interval; relax precision by factor of 2
			Plot2D'adaptive(func, var, {a, a1, b}, {fa, fa1, fb}, depth-1, epsilon*2),	// Tail() omits the pair {b, fb}
			Tail(Plot2D'adaptive(func, var, {b, b1, c}, {fb, fb1, fc}, depth-1, epsilon*2))
		)
	);
];

];	// LocalSymbols()