# Calculus

This set of exercises involves using Maxima to analyze differentiation,
integration, and the evaluation of limits.

## Derivatives

First we review a small set of the many functional forms for which 
Maxima can provide derivatives. The next four cells relate to first
derivatives of functions of a single variable.

### Simple Derivatives

First, take the derivative of a polynomial. Write the original 
expression something like y:a*x^n + b*x^(n-2) + c;. You can 
specify the exponents as you like and you may add as many 
terms as you like, as long as the number of terms with nonzero
exponents exceeds 2.

/*Input. Execute (ctrl-enter) for results */
y: a*x^n + b*x^(n-2) + c; diff(y,x);

Differentiate sin(%theta) with respect to the angle, %theta.

/*Input. Execute (ctrl-enter) for results */
diff(sin(%theta),%theta);

Define a growth function as a*b*exp(r*t) or, equivalently, as 
a*b*%e^(r*t). Note the % before e, to distinguish Euler's constant
(2.718...) from the letter e. Confirm that the derivative with respect
to time is itself a growth function, such that positive values of 
a and r  are positive.

/*Input. Execute (ctrl-enter) for results */
a*b*exp(r*t); 
/*
The use of a list is not required. In the command
diff(%,t), the % refers to the expression before 
the list of commands. You cannot use % inside a
list.
*/
[a*b*%e^(r*t), diff(%,t)];

State the function f(x):=a*log(x) and instruct Maxima to take its 
derivative with respect to x. 

/*Input. Execute (ctrl-enter) for results */
f(x):=a*log(x); diff(f(x),x);

Create the named expression, y: 5000 -20*x + x^2 - 0.05*x^3;. 
How does observation suggest that y should be have as x 
increases, for positive values of x? Take the derivative with 
respect to x. Does this expression look like one that would behave
as you expect?

/*Input. Execute (ctrl-enter) for results */
y: 5000 -20*x + 2*x^2 - 0.01*x^3;
dydx : diff(y,x);

The graph below shows how the slope of y changes as x increases.
As before, just use the graph for interpretation without concern 
regarding its commands.

Note that generally it is not a good idea to put the function and 
its derviative(s) on the same graph; their units are different. 
The vertical axis on the y expression is units of y, whatever 
that for the purpose at hand. The vertical axis on the second
axis is in units of x per unit of x, whatever x's units might be.

yGraph: gr2d(yrange=[0,15000],ytics=5000,
  explicit(y,x,0,150))$
dydxGraph:gr2d(xaxis=true,explicit(dydx,x,0,150))$
wxdraw(yGraph,dydxGraph)$

Determine the expression for the second derivative. Set that 
expression equal to zero and confirm that the first derivative,
graphed above, reaches its maximum value when x = 200/3.

/*Input. Execute (ctrl-enter) for results */
diff(y,x,2); solve(%,x);

### Partial and Cross-Partial Derivatives

Functional expressions can involve any number of variables and 
functional forms can be combined. Suppose that z = a*x^b*y^c*exp(r*t).
Bind the expression on the right to the name z, and then take 
the three derivatives with respect to x, y, and z.  Use kill(all) to 
remove assignments made above.

Assume that a, b, c, and r are positive. In which direction does an 
increased value of y cause the partial derivative of z with respect
to x to move? What about an increase in t? 

In which direction does an increase in x or y cause the trend 
growth rate (the derivative with respect to time) to change?

/*Input. Execute (ctrl-enter) for results */
kill(all)$
z : a*x^b*y^c*exp(r*t);
[diff(z,x),diff(z,y),diff(z,t)];

Define a new variable z1 by substituting these values into z:
[a=5,b=0.5,c=1.4,r=0.1]. Take the three partial derivatives of this
function. 

Does the derivative of z1 with respect to x increase or decrease
as x increases? Is it ever negative, for positive x values? How does
it react to changes in y? To changes in t?

Does the derivative of z1 with respect to y increase or decrease
as y increases?  How does it react to changes in y? To changes in t?

Answer the same three questions for the third partial derivative
in the list.

/*Input. Execute (ctrl-enter) for results */
z1 :subst([a=5,b=0.5,c=1.4,r=0.1],z);
[diff(z1,x),diff(z1,y),diff(z1,t)];

Have Maxima produce the three second partial derivatives and 
interpret them. 

/*Input. Execute (ctrl-enter) for results */
[diff(z1,x,2),diff(z1,y,2),diff(z1,t,2)];

Have Maxima produce the six cross-partial derivatives and interpret 
them. Are all signs (all positive) as you expected?

/*Input. Execute (ctrl-enter) for results */
[diff(z1,x,1,y,1),diff(z1,x,1,t,1),diff(z1,y,1,x,1),
diff(z1,y,1,t,1),diff(z1,t,1,x,1),diff(z1,t,1,y,1)];

## Integration

Often we begin with a function that is the derivative of another
function, and we wish to determine the nature of the original 
function. Integration provides the requisite information. Integration
can involve a function of a single variable or of many variables. Often,
it is possible to determine the desired integral using analytical means;
sometimes we must settle for numerical approximations.

### Integration with a Single Variable

Suppose that the known relationship is this: a*x^b - (c+d*x).
Use the integrate function to determine the original function's form.
Use this command set: assume(b>0)$ integ1: integrate(a*x^b-(c+d*x),x);.
(Try it without assume(b>)$ first, in order to see why providing this
information to Maxima is a good idea.

Notice that Maxima does not return the obligatory constant of integration.
It's up to you, the analyst, to remember this aspect.

/*Input. Execute (ctrl-enter) for results */
kill(all)$ assume(b>0)$ 
integ1: integrate(a*x^b - (c + d*x),x);

Substitute these values into the expression above: a=1000, b=-2.5,
c=50, and d=15.

/*Input. Execute (ctrl-enter) for results */
subst([a=1000,b=-2.5,c=50,d=15],integ1);

An important use of integration is to determine the area under a 
curve (or between two curves) over some range of values. Use the
command integ2: integrate(a*x^b - (c + d*x),x,x1,x2); to determine
the expression for the definite integral for this expression. 

/*Input. Execute (ctrl-enter) for results */
integ2: integrate(a*x^b - (c + d*x),x,x1,x2);

Determine the area given the parameter values used above, and with
x1 = 1 and x2 = 5.

/*Input. Execute (ctrl-enter) for results */
subst([a=1000,b=-2.5,c=50,d=15,x1=1,x2=5],integ2);

Determine the values of the integral at the upper and lower ends
of the specfied integral and confirm that the value of the 
definite integral above is the difference between the two values.

/*Input. Execute (ctrl-enter) for results */
lower: subst([a=1000,b=-2.5,c=50,d=15,x=1],integ1);
upper: subst([a=1000,b=-2.5,c=50,d=15,x=5],integ1);
upper-lower;

The graph below shows the expression that is being evaluated
and identifies the area as the area above the x axis between 
1 and approximately 2.5 less the area below the x axis over 
the rest of the range of x values.

wxdraw2d(yrange=[-200,2000],xaxis=true,
explicit(1000*x^-2.5 - (50 + 15*x),x,0.1,6),
line_type=dots,points_joined=true, 
points([[1,-200],[1,2000]]),points([[5,-200],[5,2000]]) 
)$

The graph below decomposes the expression, so that the 
area is the area between the thin blue line and the 
thicker black line between x = 1 and x = 2.5, 
approximately, less the area between the black line
and the blue line over the rest of the x range.

wxdraw2d(yrange=[0,2000],xaxis=true,
line_width=2,
explicit(1000*x^-2.5,x,0.1,6), color=black,
explicit(50 + 15*x,x,0,6),line_type=dots,
points_joined=true, points([[1,0],[1,2000]]),
 points([[5,0],[5,2000]])   )$

The graph above emphasizes that the area defined above is the 
difference between the areas under each of two expressions.
Separate the function into its two component parts, determine
the area under each between x = 1 and x = 5, and then calculate
the difference between the two. 

/*Input. Execute (ctrl-enter) for results */
area1: integ1: integrate(a*x^b,x,1,5);
area2: integ1: integrate(c + d*x,x,1,5);
netarea: area1-area2;
subst([a=1000,b=-2.5,c=50,d=15],[area1,area2,netarea]);

Maxima offers a number of numerical options that can approximate
areas when integrate cannot determine the areas analytically. 
Consider the Romberg option, which is most like integrate in 
syntax. Replace integrate with romberg in the commands above.
Rename the areas to indicate that a different technique is being used.

Because romberg is numerical, the numbers must be substituted before 
the computation begins. You will discover that the Romberg approximation
is quite close to the analytical solution. This is often, but not always, 
true.

/*Input. Execute (ctrl-enter) for results */
[expr1,expr2]:subst([a=1000,b=-2.5,c=50,d=15],[a*x^b,c+d*x]);
[areaR1: romberg(expr1,1,5), areaR2: romberg(expr2,x,1,5),
netareaR: areaR1-areaR2];

### Integration with More than One Variable

Create the expression z:a*t*x^b - (c*t + d*x );

/*Input. Execute (ctrl-enter) for results */
kill(all)$ z:a*t*x^b - (c*t + d*x );

Integrate z with respect to x and then integrate the 
result with respect to t. Name the result integr1.
Assume that b < -1.

/*Input. Execute (ctrl-enter) for results */
assume(b<-1)$
integrate(z,x); integr1: integrate(%,t);

Substitute the following into z: a=150,b=-0.5,c=50,d=15.

/*Input. Execute (ctrl-enter) for results */
params:[a=150,b=-0.5,c=50,d=15];
subst(params,integr1);
/*or*/
subst([a=150,b=-0.5,c=50,d=15],integr1);

Given the parameters above, create an expression for z. Assign
the name z1 to this expression.

/*Input. Execute (ctrl-enter) for results */
z1: subst(params,z);

The graph below shows the expression for values of t ranging from 0 to 5
and values of x ranging from 1 to 5. It shows that higher values of x 
reduce the expression's value, while higher values of t increase the
value. Over most of the range, values are positive. The number that
we seek is area between the green surface and the gray floor (where
z = 0).

wxdraw3d(xlabel="x",ylabel="t",zlabel="z",
key = "The floor: x = y =0",explicit(0,x,0,5,t,0,5), 
color=green, key="z = -15*x + 150*t/x^0.5 - 50*t",
explicit(z1,x,1,5,t,0,5),view=[80,320],
dimensions=[640,640]);

Integrate the expression z1 over the ranges shown above, t from 0 to 5 and 
x from 1 to 5. Preface the integrate command(s) with ratprint:false to avoid
a list of substitutions that Maxima makes in determining this value.

/*Input. Execute (ctrl-enter) for results */
ratprint:false$ 
/*
The float(%); command is needed for comparison with
the numerical approximation below.
*/
integrate(z1,t,0,5); integrate(%,x,1,5); float(%);
/*
or the order can be reversed
*/
integrate(z1,x,1,5); integrate(%,t,0,5); float(%);
/*
or the commands can be nested
*/
integrate(integrate(z1,x,1,5),t,0,5); float(%);

Replace integrate with romberg to confirm that the numerical approximation
yields virtually the same value as the analytical solution.

/*Input. Execute (ctrl-enter) for results */
romberg(romberg(z1,x,1,5),t,0,5);

## Limits

Maxima's limit function can be used to determine numeric values
that sepcific functions tend toward as a variable approaches some
value that interests the analyst. Use limit(1/x,x,0); to determine
the limit of 1/x as x approaches zero.

/*Input. Execute (ctrl-enter) for results */
limit(1/x,x,0);

The result above is not as informative as it seems. The text reply
infinity denotes complex infinity. Among other possibilities, this 
could indicate an arbitrarily large positive number ("positive
infinity") or an arbitrarily large negative number ("negative infinity").
Which of these occurs depends on whether x approaches 0 from above
or from below.

Enter the commands [limit(1/x,x,0,plus),limit(1/x,x,0,minus)]; to 
retrieve the limits as x approach zero from above and from below.

Placing the two commands in a list is not required.

/*Input. Execute (ctrl-enter) for results */
[limit(1/x,x,0,plus),limit(1/x,x,0,minus)];

The graph illustrates what's going on in the example above.

wxdraw2d(yrange=[-10,10],xaxis=true,
yaxis=true,key="1/x",explicit(1/x,x,-5,5));

An important application of the limit concept is the derivation of
rules for taking derivatives. Define the function f(x):= a + b*x^2 + c*x^n;.
Then define the expression (f(x+h) - f(x))/h; where h is an increment
to x. Now use the command limit(%,h,0); to evaluate this limit.
Take the first derviative of f(x) to confirm that the two approaches 
are the same for this function. (The two results are not presented 
in quite the same way; confirm that they are equivalent.)

/*Input. Execute (ctrl-enter) for results */
f(x):= a + b*x^2 + c*x^n;
(f(x+h) - f(x))/h;
limit(%,h,0);
diff(f(x),x);

Confirm that using the limit approach and taking the 
derivative yields the same result for g(x)=a*exp(B*x).
You will have to answer this question: Is x B an integer?
A list of possible answers is provided. Select no, n, or N 
and ctrl-enter or shift-enter. (You might have to hit
ctrl-enter or shift-enter twice.)

/*Input. Execute (ctrl-enter) for results */
g(x):= A*exp(B*x);
(g(x+h)-g(x))/h;
limit(%,h,0);
diff(g(x),x);

Once you have told Maxima that x*B is not an integer, this information
stays in its memory until you restart Maxima. Apply the steps above 
to g2(x) = x^A*exp(B*x).

Confirm that the third and fourth output lines are equivalent.

/*Input. Execute (ctrl-enter) for results */
g2(x):= x^A*exp(B*x);
(g2(x+h)-g2(x))/h;
limit(%,h,0); expand(%);
diff(g2(x),x);

A Taylor series converts any differentiable function into a polynomial
of some specified order. Use the command taylor(sqrt(x), x, 9, 3); to
approximate the square root of 7 with a third-order polynomial.
The result appears as a ploynomial +..., where the +... indicates 
that a residual exists, so that the result is contains an error.

Substitute 7.0 for x. Use the command float(sqrt(7.0)); to compare the 
value that Maxima generates with the third-order Taylor approximation.

This expansion around x = 9 can provide good results for values of x that
are also near to 9. It should not be used for x values that are far from 9.

/*Input. Execute (ctrl-enter) for results */
taylor(sqrt(x),x,9,3);
subst(7.0,x,%);
sqrt(7.0);

Copy the cell from above and paste it below. Then place
parentheses around the first line and convert it to the 
command taytorat(taylor(sqrt(x),x,9,3)); which expresses
the polynomial in canonical form and removes the +... term.

/*Input. Execute (ctrl-enter) for results */
taytorat(taylor(sqrt(x),x,9,3));
subst(7.0,x,%);
sqrt(7.0);

The graph below provides a sense of how close the Taylor
approximation is, given the expansion is about x = 9.

wxdraw2d(
point_size=2,points_joined=true, line_type=dots,
points([[0,3],[9,3],[9,0]]), line_type=solid,
explicit(sqrt(x),x,0,25),key="Taylor",
line_width=2,
explicit(taytorat(taylor(sqrt(x),x,9,3)),x,0,25)
)$

Copy the command cell above the graph cell and place it below this cell.
Replace 3 with 5. Confirm that the result is that taylor provides a 
more precise approximation of sqrt(5).

/*Input. Execute (ctrl-enter) for results */
taylor(sqrt(x),x,9,5);  taytorat(%);
float(subst(7,x,%));
float(sqrt(7));

Confirm that the series in the first line of output is this:
f_0(9)/0! + f_1(9)/1! + f_2(9)/2! + ... + f_5(9)/5!, where f_i (9) is the i-th derivative
of f(x), evaluated at x = 9. The zero-th derivative is the original function. See
notes at: http://www.sosmath.com/calculus/tayser/tayser01/tayser01.html.

Use these commands:
 makelist(diff(sqrt(x),x,i), i,0,5)$ to determine the six derivatives indicated above.
 diffList: subst(x0,x,%); to name the list and to evaluate it at x = x0.
 x0xList: makelist((x-x0)^i,i,0,5); to create the list that raises the difference 
   between x0 and x to the n-th power.
 facList: makelist(i!,i,0,5); to make a list of factorials for i = 0, 1, ..., 5.
 x0xList*diffList/facList; to create the series of expressions, stated in terms of x0.
subst(9,x0,%); to evaluate the preceding expression when x0 = 9.

/*Input. Execute (ctrl-enter) for results */
makelist(diff(sqrt(x),x,i), i,0,5)$ diffList: subst(x0,x,%);
x0xList: makelist((x-x0)^i,i,0,5);
facList: makelist(i!,i,0,5); 
x0xList*diffList/facList;
subst(9,x0,%);

## Differential Equations

Define an expression for a derivative named growth as follows:
growth: 'diff(y(t),t)= r*y(t);

Be sure to precede diff( ) with the apostrophe that inhibits 
computation.

kill(all)$
/*Input. Execute (ctrl-enter) for results */
growth: 'diff(y(t),t)= r*y(t);

Generate a solution for this differential equation using desolve( ).
Name the solution soln1.

/*Input. Execute (ctrl-enter) for results */
soln1: desolve(growth,y(t));

Confirm that the solution satisfies the expression for growth.

/*Input. Execute (ctrl-enter) for results */
subst(soln1,growth);
/*Check this by using diff( ). Denote y(0) as y0, a constant.*/
diff(y0*exp(r*t),t); /* or */ diff(y0*%e^(r*t),t);

diff(y0*%e^(r*t),t);

Suppose that y(0) is not known, but y(50) = 1000. Use atvalue( )
to express y(75) as a multiple of y(0).

/*Input. Execute (ctrl-enter) for results */
atvalue( y(t), t = 50, 1000)$
soln2 : desolve(growth, y(t));
subst([r=0.05,t=75],soln2);

Given y(50) and a value of r, y(0) can be determined. Do this
by setting the expression for y(t) equal to 1000.

/*Input. Execute (ctrl-enter) for results */
soln3: solve(1000 = y(0)*%e^(0.05*50),y(0));

Create a functional expression f(t, r) for this growth function.
Evaluate this function for r = 0.05 and for t = 50 and t = 75.

/*Input. Execute (ctrl-enter) for results */
y0: float(rhs(soln3[1]));
f(t,r) := ''y0*exp(r*t);
[f(50,0.05), f(75,0.05)];

wxdraw2d(yrange = [0, 1.05*f(75,0.05)],
explicit(f(t,0.05),t,0,75))$

One important application of the example above is the Malthusian
prediction. Malthus expressed population growth as b*y - m*y, where
b is the birth rate and m is the mortality rate. Malthus's system contained
another simple differential equation diff(x(t),t) = a, where a is the rate 
at which the means to feed and maintain a population grows. Think of 
x as the number of people that can be maintained and y as the number
of people. 

We can solve the two equations as a pair. Name the subsistence 
differential equation food with this command: 
food: 'diff(z(t),t) = a;. This says that food (or subsistence goods
more generally) grows by a unit each time period, while the 
population is growing in an exponential fashion.

/*Input. Execute (ctrl-enter) for results */
food: 'diff(z(t),t) = a;

Solve the system of equations as follows:
solnyz: desolve([growth,food],[y(t),z(t)]);, where
solnyz is 

/*Input. Execute (ctrl-enter) for results */
solnyz: desolve([growth,food],[y(t),z(t)]);

Confirm that if r > 0 (b > m, in Malthusian terms) and a > 0 
but finite, then at some finite time period y grows sufficiently 
to use up all available product. That is, at some point y(t)=z(t).
kill(y0)$
solve(subst(y0,y(0),rhs(solnyz[1]))
    - subst(z0,z(0),rhs(solnyz[2])),t);

/*Input. Execute (ctrl-enter) for results */
kill(y0)$
/*Input. Execute (ctrl-enter) for results */
solve(subst(y0,y(0),rhs(solnyz[1]))
    - subst(z0,z(0),rhs(solnyz[2])),t);

The graph shows a case in which population growth 
catches up with food growth at t = 914.23. To draw the 
two graphs, the expressions growthcurve and foodcurve.
The initial value of y is 10 and the growth rate is 1%
per period. The initial value of z is 2000 and the 
growth rate is 100 units per period. (We assume that 
one unit of z is required per person.)

[growthcurve: subst([y(0)=10,r=0.01],rhs(solnyz[1])),
foodcurve: subst([z(0)=2000,a=100],rhs(solnyz[2])),
find_root(growthcurve-foodcurve,0.1,10000)];
wxdraw2d(
explicit(growthcurve,t,0,1000),
explicit(foodcurve,t,0,1000))$

An alternative to modeling a population growth, one that
is consistent with Malthusian concerns is to posit that the
population inhabits an envrionment that imposes a limit
on that population. This limit is often called the carrying
capacity and is denoted by K below. 

One model of population growth is the logistic curve.
The next few cells look briefly at this model. The material
cited here is for your information and goes beyond the 
material in the text.
http://people.uncw.edu/hermanr/mat367/FCABook/FCA_Main2013.pdf
http://people.uncw.edu/hermanr/mat463/ODEBook/Book/ODE_Main.pdf

The first cell provides the expression for the population 
growth function and binds that expression to the name
dPdt.

/*Input. Execute (ctrl-enter) for results */
kill(all)$          
dPdt : r*P*(1-P/K);

We do not use Maxima's analytical approach to analyzing this
expression. Rather we use an alternative, which can be used
where analytical solutions cannot be achieved. We use the
Runge-Kutta (rk) method of tracing the population size.
Doing so requires that r, P0 (the inital value of P), and K
be specified. The first command below makes the 
requisite substitutions and assigns the result the name
dPdt0.

Next a list of data is produced using the rk( ) command. 
The arguments in this command are these: the expression,
the variable whose values is being determined, that 
variable's initial value, and a list. The list defines the
independent variable (t in this case), its inital value,
its final value, and the increment size. The specification
below considers P for t ranging from 0 to 150 at 15 period
intervals. 

The transpose(matrix()) command is used only to align
the results.

/*Input. Execute (ctrl-enter) for results */
dPdt0 : subst([r=0.05,K=200],dPdt)$
data : rk(dPdt0,P,50,[t,0,150,15])$
transpose(matrix(data));

The graph shows the 16 points defined above and 
confirms that the population size approaches K = 200
asymptotically.

wxdraw2d(xlabel="t",ylabel="P",
point_size=2,points_joined=true,
points(data))$

This growth function can be solved analytically. See
http://math.usu.edu/~powell/biomath/mlab3-02/node3.html
for a derivation.

P(t,r,P0,K):= K*P0/(P0 + (K-P0)*%e^(-r*t));

The graph below shows a continuous representation of the 
expression.

wxdraw2d(xlabel="t",ylabel="P",
explicit(P(t,0.05,50,200),t,0,150))$