Numerical Calculations

>	nextprime(1007898743453453465767945830458254303);

>	evalf(Pi,100);

Latex command

Latex snippets are easily generated with the latex command in Maple.

>	Int(1/(x^2+1), x) = int(1/(x^2+1), x);

>	latex(Int(1/(x^2+1), x) = int(1/(x^2+1), x));

\int \! \left( {x}^{2}+1 \right) ^{-1}{dx}=\arctan \left( x \right)

>

Introduction to the CurveFitting Package

>

restart;

The following command allows you to use the short form of the function names in the CurveFitting package.  Note that the PolynomialInterpolation , RationalInterpolation , Spline , and ThieleInterpolation functions replace the obsolete functions interp, ratinterp, spline, and thiele, respectively.

>	with(CurveFitting);

A feature of the CurveFitting package is that the routines allow you to specify the data points in two ways.  The points can be provided as two lists, the first containing the independent values and the second containing the dependent values.  The following command computes a polynomial that interpolates the points {(0,1), (1,3), (2,-1), (3,2)}.

>	PolynomialInterpolation([0,1,2,3], [1,3,-1,2], v);

The points can also be provided as a list of lists.

>	PolynomialInterpolation([[0,1],[1,3],[2,-1],[3,2]], v);

The ability to specify the data as a list of lists facilitates plotting the original data with the result.  Example plots are provided in the next section.  The following command allows you to use the short form of the function names in the plots package.

>	with(plots):

Warning, the name changecoords has been redefined

Fitting a Curve Through a Set of Points

Below, we define a list of 8 points, then plot these points by using the plots[pointplot] function.  (Note that the short form of the function name is used below.)

>	points1 := [[0,1],[1,2.5],[3,2.3],[4.2,5],[5,3.5],[5.8,4.2],[7,7],[8,10]]: pntplot1 := pointplot(points1,symbol=BOX):

The next command uses the PolynomialInterpolation function to compute the unique polynomial that interpolates the data points.

>	polycurve := PolynomialInterpolation(points1, v):

The plot command shown below generates a plot of the result. The plots[display] function displays both this plot and the previously computed plot of the data points.  It can be seen that the degree 7 polynomial passes through each of the 8 data points.

>	polyplot := plot(polycurve,v=0..8): display([pntplot1,polyplot]);

If n data points are supplied, then the interpolating polynomial has degree less than or equal to n-1.  One well-known problem with fitting an interpolating polynomial through a large set of data points is the undesirable oscillations produced by a high-degree polynomial.  In the plot above, the function value for v=2 is far from the value for v=1 or v=3.

A smoother curve can be obtained by using the Spline function.  This function returns a piecewise polynomial of default degree 3 that passes through the 8 data points.   The second derivative is set to zero at the endpoints, and this results in a "natural" spline.  The degree of the spline can be controlled by the degree option, described on the Spline help page.

>	splcurve := Spline(points1, v);

The plot below shows that a smoother curve, without the unwanted oscillations, is produced.

>	splplot := plot(splcurve,v=0..8): display([pntplot1,splplot]);

Computing a Least Squares Fit

The LeastSquares function can be used to find a curve that best fits the data in a least-squares sense, that is, minimizes the sum of the squares of the differences between the estimated values and the actual data.  Unlike the curves described in the previous section, the least-squares curve may not necessarily pass through the given points.  The least-squares method is often used to fit models to experimental data.

The commands below show how the LeastSquares function is used to compute the best linear fit through the points defined at the beginning of the previous section.  In the resulting plot, it can be seen that the curve passes near but not necessarily through all the points.  

>	lscurve := LeastSquares(points1,v); lsplot := plot(lscurve,v=0..8): display([pntplot1,lsplot]);

In the previous example, the default curve used to fit the points is a linear polynomial av+b with parameters a and b.  The LeastSquares function allows you to provide a different type of curve and to specify the parameters to optimize.   In addition, weights associated with the data points can be defined.  See the LeastSquares help page for more details about these options.

The LeastSquares function requires that the model curve be linear in the parameters.  In some cases, a nonlinear model can be transformed to allow a least-squares fit.  For example, if you wish to use the model

w=a*exp(b*v),

you can take the logarithm of both sides of the equation and apply the transformation

{y=log(w), c=log(a)}

to obtain the model

y=c+b*v,

which is linear in the parameters c and b.

>	newpoints1 := map(x->[x[1],evalf(log(x[2]))], points1); newlscurve := LeastSquares(newpoints1,v);

This result gives values for c and b.  The following commands show the resulting plot, after transforming back to the original model.

>	a := exp(coeff(newlscurve,v,0)); b := coeff(newlscurve,v,1); newlsplot := plot(a*exp(b*v),v=0..8): display([pntplot1,newlsplot]);

Introduction to the Optimization Package

>

restart;

The following command allows you to use the short form of the command names in the Optimization package.  

>	with(Optimization); with(plots):

Warning, the name changecoords has been redefined

The commands in the Optimization package allow you to specify the objective function and the constraints in several different ways.  This worksheet deals with the algebraic form of the objective function and constraints.  The other forms are Matrix form and procedure form where Matrix form is covered in a separate example worksheet.

In general, an optimization problem is of the form

The Optimization package handles several types of objective and constraint functions.  Each type of problem is a subtype of the next one.

Objective Function/Constraint/Command

Linear/Linear/LPSolve

Quadratic/Linear/QPSolve

Nonlinear(General)/Nonlinear(General)/NLPSolve

There is also a minimize command that will attempt to automatically select the best algorithm for a given problem.

In addition, there is a command, LSSolve, which minimizes functions of the form

Linear Programming Example 1

For a first example, we have a simplex two dimensional linear programming problem.

subject to

>	obj := -2*x-y;

>	cnsts := [y<=4*x+1/2,y<=-5*x+2,x>=0,y>=0];

The plot below shows the feasible region in yellow, the contours of the objective function, and the optimal solution as a green circle.

>	p1 := inequal(cnsts, x=-0.5..2, y=-0.5..2, optionsexcluded=(colour=white), optionsfeasible=(colour=yellow)): p2 := contourplot(obj, x=-0.5..2, y=-0.5..2): p3 := pointplot({[0.1666,1.166]}, symbolsize=13, colour=green): display(p1, p2, p3);

The  command returns the optimal function values, as well as the point at which the optimal value occurs.

>	LPSolve(obj, cnsts);

Alternatively, we could use the first two constraints and the nonnegative option.

>	LPSolve(obj, cnsts[1..2], assume=nonnegative);

The first element of the solution is the minimum value that the objective function obtains, while satisfying the constraints.  The second element indicates a point where the minimum is reached.  This point is not necessarily unique.

Linear Programming example 2

We can also include equality constraint as the next example shows.  This example also demonstrates the maximize option.

>	obj := 2*x+y+z;

>	cnsts := [y<=4*x+1/2,y<=-5*x+2,z=-2/3*x-2/3*y+20/3];

>	LPSolve(obj, cnsts, assume=nonnegative, maximize);

Quadratic Programming Example

For a first example of a quadratic program, we will use:

>	obj := 4*x^2-y+x+5;

>	cnsts := [x+y-7*x>=-11, 11/2*x+y<=0,x>=-4];

The plot below shows the feasible region in yellow and the contours of the objective function.  The green circle indicates the optimal point.

>	p1 := contourplot(obj, x=-4..4, y=-10..10, contours=30): p2 := inequal(cnsts, x=-4..4, y=-10..10, optionsexcluded=(colour=white), optionsfeasible=(colour=yellow)): p3 := pointplot({[-0.8125,4.486]}, symbolsize=13, colour=green): display(p1,p2,p3);

The QPSolve command returns the optimal function values, as well as the point at which the optimal value occurs.

>	QPSolve(4*x^2-y+x+5, {x+y-7*x>=-11, 11/2*x+y<=0,x>=-10});

Nonlinear Programming

The NLPSolve command allows us to solve problems of a fairly general nature.  As long as the derivative of the objective function and constraint function is continuous, we can use NLPSolve.

Our first basic example of nonlinear programming is a simple univariate minimization of a polynomial.

>	p := 0.4126984123e-2*x^7-.1091666666*x^6+7.078095236*x+1.136388888*x^5-18.19500000*x^2-5.845833331*x^4+15.23138888*x^3+2;

This plot shows the function and the minimum found by NLPSolve.

>	p1 := plot(p, x=4.5..5.5): p2 := pointplot({[4.8877,0.6883]}, symbolsize=13, colour=green): display(p1, p2);

We can minimize this function in the range [0..6].

>	NLPSolve(p, x=0..6);

We can also maximize the function.

>	NLPSolve(p, x=0..6, maximize);

We can minimize more complex functions such as:

>	f := x*erf(x)*cos(x);

>	s := NLPSolve(f,x=1..20);

Below we plot the function f and the minimum point.  We can see that NLPSolve finds a local minimum, not necessarily a global one.

>	p1 := plot(f,x=1..20): p2 := pointplot([[rhs(s[2][1]),s[1]]],color=green,symbolsize=14): display(p1,p2);

We can also minimize multivariate functions such as

>	f := 100*(y-x^2)^2+(1-x)^2;

with

>	NLPSolve(f, x=-10..10, y=-10..10, initialpoint=[x=-1.2, y=1]);

Curve Fitting with Least Squares

>	restart: with(Optimization):

We can use the least squares command to do curve fitting.  For example, if we have the following data points, ( ):

>	data := [[1,1.5],[2,3.5],[2.5,1.9],[3.1,4.5],[4.3,1.9],[4.7,2.4],[5.8,3],[6,5.7],[6.6,3.4],[6.7,0.5],[7.1,4],[7.4,1.7]]:

and we want to fit the data with a quintic polynomial.  Let our quintic be,

.

Now we need to find a, b, and c such that

is a minimum.

LSSolve takes as input a list of the residues .

>	p := a*x^5+b*x^4+c*x^3+d*x^2+e*x+f;

>	residues := map((d) -> eval(p, x=d[1])-d[2], data);

Now we can call LSSolve with the residues.

>	sol := LSSolve(residues);

Our polynomial becomes:

>	poly := eval(p, sol[2]);

We can now plot the data and the curve

>	with(plots): p1 := pointplot(data): p2 := plot(poly, x=0.9..8): display(p1, p2);

Warning, the name changecoords has been redefined

>