Code Example: Piecewise Cubic Approximator

#include "pw_cubic_approximator.h"  	// for the PWCubicApproximator class 

	typedef PWCubicApproximator<double, double>		PWCA_dd; 

// prototypes 

	template <class IT, class OT> 
	void	approx_test_midpoints(
		const char	   *title, 
		const char	   *f_name, 
		OT		  (* f)(IT), 
		const char	   *a_name, 
		/* const */ PWCA_dd  & approx ); 

int 
main()
{
	// Now fit a piecewise cubic to exp(). 
	// Find out how close the fit is, and how close together the knots are. 

	PWCA_dd		ax_pwc_dd5(exp, 1.0, 0.01); 

	ax_pwc_dd5.extrapol_param = 0.5; 
	ax_pwc_dd5.map_range(-1, 4.5); 

	approx_test_midpoints("Test of deviation between the knots:", 
	                      "exp", exp, "ax_pwc_dd5", ax_pwc_dd5); 

	printf("Number of function evaluations: %d\n\n", ax_pwc_dd5.nofun); 
}

template <class IT, class OT> 
void 
approx_test_midpoints(
	const char	   *title, 
	const char	   *f_name, 
	OT		  (* f)(IT), 
	const char	   *a_name, 
	/* const */ PWCA_dd  & approx ) 
{
	cout << title << endl 
	     << endl ; 

	PWCA_dd::member_interpolator::Dataset::const_iterator  i, 
	                                                       inext; 

	if (approx.interp.dataset.size())
		for (i = approx.interp.dataset.begin(), inext = i, ++ inext; 
		     inext != approx.interp.dataset.end(); 
		     ++ i, ++ inext) 
			{
			double		x = ((*inext).first + (*i).first) / 2; 

			if (   approx.get_mapped_range().find(x) 
			    != approx.get_mapped_range().end()) 
				{
				double  y = f(x), 
				        y_star = approx(x); 
				double  dy = y_star - y; 
							
				printf("%s(%6g) = ", a_name, x); 
				cout << y_star << "; " << f_name << "(x) = " << y ; 
				cout << "; dy = " << dy << endl ; 
				}
			}
		
	putchar('\n'); 
}
	

Output:


Test of deviation between the knots:

ax_pwc_dd5(  -2.5) = 0.089973; exp(x) = 0.082085; dy = 0.00788796
ax_pwc_dd5( -1.75) = 0.174423; exp(x) = 0.173774; dy = 0.000648699
ax_pwc_dd5( -1.25) = 0.286076; exp(x) = 0.286505; dy = -0.000428518
ax_pwc_dd5( -0.75) = 0.47166; exp(x) = 0.472367; dy = -0.000706506
ax_pwc_dd5( -0.25) = 0.777636; exp(x) = 0.778801; dy = -0.00116483
ax_pwc_dd5(  0.25) = 1.2821; exp(x) = 1.28403; dy = -0.00192048
ax_pwc_dd5(  0.75) = 2.11383; exp(x) = 2.117; dy = -0.00316634
ax_pwc_dd5(  1.25) = 3.48512; exp(x) = 3.49034; dy = -0.00522041
ax_pwc_dd5( 1.625) = 5.07937; exp(x) = 5.07842; dy = 0.000952265
ax_pwc_dd5( 1.875) = 6.52022; exp(x) = 6.52082; dy = -0.000600116
ax_pwc_dd5( 2.125) = 8.37213; exp(x) = 8.3729; dy = -0.000770564
ax_pwc_dd5( 2.375) = 10.75; exp(x) = 10.751; dy = -0.000989424
ax_pwc_dd5( 2.625) = 13.8033; exp(x) = 13.8046; dy = -0.00127045
ax_pwc_dd5( 2.875) = 17.7238; exp(x) = 17.7254; dy = -0.00163128
ax_pwc_dd5( 3.125) = 22.7578; exp(x) = 22.7599; dy = -0.00209461
ax_pwc_dd5( 3.375) = 29.2216; exp(x) = 29.2243; dy = -0.00268953
ax_pwc_dd5( 3.625) = 37.5213; exp(x) = 37.5247; dy = -0.00345343
ax_pwc_dd5( 3.875) = 48.1783; exp(x) = 48.1827; dy = -0.00443429
ax_pwc_dd5( 4.125) = 61.8621; exp(x) = 61.8678; dy = -0.00569374
ax_pwc_dd5( 4.375) = 79.4325; exp(x) = 79.4398; dy = -0.00731091
ax_pwc_dd5(4.5625) = 95.8237; exp(x) = 95.8227; dy = 0.00101148
ax_pwc_dd5(4.6875) = 108.581; exp(x) = 108.581; dy = -0.000622117
ax_pwc_dd5(4.8125) = 123.038; exp(x) = 123.039; dy = -0.000704951
ax_pwc_dd5(4.9375) = 139.42; exp(x) = 139.421; dy = -0.000798814
ax_pwc_dd5(5.0625) = 157.984; exp(x) = 157.985; dy = -0.000905175
ax_pwc_dd5(5.1875) = 179.019; exp(x) = 179.02; dy = -0.0010257
ax_pwc_dd5(5.3125) = 202.856; exp(x) = 202.857; dy = -0.00116227
ax_pwc_dd5(5.4375) = 229.865; exp(x) = 229.867; dy = -0.00131702
ax_pwc_dd5(5.5625) = 260.472; exp(x) = 260.473; dy = -0.00149238
ax_pwc_dd5(5.6875) = 295.153; exp(x) = 295.155; dy = -0.00169109
ax_pwc_dd5(5.8125) = 334.452; exp(x) = 334.454; dy = -0.00191626
ax_pwc_dd5(5.9375) = 378.99; exp(x) = 378.986; dy = 0.00327527

Number of function evaluations: 62
	

Notice that these deviations are all well within the chosen tolerance of 0.01.


Copyright © 2001  J. E. Brown   all rights reserved.
write me here
              
Los Alamos, NM  USA