This workbook contains the wxMaxima commands that appear in 
Sections 7 and 8 of Chapter 1.

The cell below is optional. The purpose of each of the commands
is indicated in the final section of Chapter 1 in the text.

load(draw)$ 
set_draw_defaults(  fill_color=gray40,
draw_realpart=false,font_size=10)$
fpprintprec:5$ ratprint:false$kill(all)$

# Chapter 1. Introduction

## Calculus

This section introduces Maxima symbolic approach to 
differential and integral calculus. 

### Differentiation

The first command specifies that dy is required, so Maxima
returns dy = (dy/dx)*dx, where Maxima's notation for dx is del(x).  
Maxima} interprets the second command as requesting dy/dx.
We place these two input items in a list and list inside matrix( ).

matrix(["dy requested", "dy/dx requested"],["- - -",  "- - -"],
[diff(log(x)/x),           diff(log(x)/x, x)]);

Multiple Variables
For z = x^a y^b, the diff command, without reference to a variable, tells 
Maxima to fully differentiate this function with respect to x, y, a, and b. 
Maxima} does not know that a and b are constants. 
The second command Maxima to produce the partial derivative of z 
with respect to x. 

z : x^a*y^b$    diff(z);         diff(z, x);

The next cell uses declare( ) to state that a and b are constants. 

declare([a, b], constant)$     diff( z );

Maxima can extract higher-order derivatives. With multiple variables,
it can extract cross-partial derivates.

Define f(x, y) = x^4/y. The example below shows 
the first and second derivatives with respect to $x$, denoted $f_x$ 
and $f_{xx}  $. It also shows the first derivative with respect to $y$
and the cross-partial derivative, $f_{yx}$.

kill(all)$    f(x,y) := x^4/y$  print( "f(kx,y) = ", f(x,y))$

We use subscripts to indicate the partial derivatives. One use of brackets 
is to generate subscripts. The output below shows the first and second
partial derivatives of f(x,y) with respect to x.

print(f[x],"=",diff(f(x,y),x),  ", and",f[xx],"=",diff(f(x,y),x,2))$

Next we generate the first partial derivative of f(x,y) with respect to y and
the cross-partial derivative. Revise these commands to generate the 
cross-partial derivative f[xy] and confirm that it equals f[yx].

print(f[y], "=", diff(f(x, y), y ), ",  
      and    ", f[yx], "=", diff(f(x,y), y,1, x, 1)  )$

### Integration

The syntax of indefinite integration resembles that of 
differentiation, with the integrate( ) command replacing 
diff( ), except that multiple integration requires multiple 
calls to integrate( ). 

The first item in the table below shows integral of f(x,y) 
with respect to (wrt) x.  The second item in the matrix below 
is a list of the results of two equivalent ways of conducting
mulitple integration.

Exercise: Confirm that reversing the order of in which x and y
are entered into integrate results in the same final expression. That is,
integrate f(x,y) wrt y and then integrate the result wrt x.

matrix(["Integral of f(x,y) wrt x", 
"Two ways to do multiple integration"],
 [Ix: integrate(f(x, y),x), 
  [integrate(Ix,y), integrate(integrate(f(x,y),x),y)] ] );

Indefinite and Definite Integrals
The integral can be evaluated over a specified range. We continue
to use f(x,y), so to determine a value requires integrating over
a range for one variable. This yields an expression. Then integrate
that expression over a range of the other variable. Sequence does
not matter.

print("x first, 0<=x<=5 yields an expression ",Ix05:integrate(f(x,y),x,0,5),
"and then y, 1<=y<=2 yields a value " ,integrate(Ix05,y,1,2))$
print(Iy02: integrate(f(x,y),y,1,2), integrate(Iy02,x,0,5))$

The cell above contains exact solutions. If we represent the 
limits of integration as floating point values, then the result
is a floating point representation of the solution.

print("x first, 0<=x<=5 yields an expression ",
 Ix05:integrate(f(x,y),x,0.0,5.0),
"and then y, 1<=y<=2 yields a value " ,
integrate(Ix05,y,1.0,2.0))$
print("The counterparts when integration order is reversed are ",
Iy02: integrate(f(x,y),y,1.0,2.0), "and ", 
integrate(Iy02,x,0.0,5.0))$

The limits of integration may include infinity and -infinity.
This time, we embed the integrals within float( ) in order to 

float([f(x):=(1/sqrt((2*%pi)))*exp(-x^2/2),
integrate(f(x),x,0,1.96 ),integrate(f(x),x,0,inf),
integrate(f(x),x,-inf,inf)]);

The text says that failure to inform Maxima that a parameter has a specific
characteristic can result in a query before it proceeds. Execute the command
below. In response to the question, type p; and ctrl-enter. You could enter
p, pos, or positive with the same result.

integrate( exp(-a*t^2),t,-inf,inf);p;

To avoid the query, provide Maxima the requisite information, using assume( ).

assume(a>0)$ integrate( exp(-a*t^2),t,-inf,inf);

Here is an example for which Maxima cannot provide a definite integral. A visual
check, using a linear approximation, says that the area should be just over
(1/2)*(1*1).

kill(all)$   expr:exp(-x^3)$ integrate(expr,x,0,1);
wxdraw2d(explicit(expr,x,0,1))$

Numerical Integration

The Romberg approximation, command romberg( ), and the quad_quags( ) method are
two alternatives for numerical integration.  The output shows that the two differ 
slightly. The quad_quags( ) output provdes a measure of accuracy that romber( ) 
does not offer.

fpprintprec:12$ 
matrix(["Romberg output","quad_quags output"],
       [romberg(expr,x,0,1),quad_qags(expr,x,0,1)]);

Like integrate( ), romberg( ) can be extended to integrals of more than one variable.
Nesting achieves the desired result. The commands must be nested. That is, 
romberg(f(x,y),y,1,2) yields no result because the output is not a number but an
expression. As always, when Maxima cannot evaluate a command, it just returns 
the command.

kill(all)$  fpprintprec:16$  f(x,y):=x^4/y;
/*invalid command */   romberg(f(x,y),y,1,2); 
romberg( romberg(f(x,y),y, 1, 2), x, 0, 5);

### Limits

Maxima can compute limits for many functions.  For, f(x) = sin(x)/x, the limit is the 
same when x approaches 0 from either above or below. 

kill(all)$   [f(x):=sin(x)/x,  limit(f(x),x,0)];
wxdraw2d(explicit(f(x),x,-1,1))$

For many functions, however, the value of the limit depends on the direction 
from which the limiting value of x is approached. Consider g(x) = |x|/x.

kill(g)$    [g(x):=abs(x)/x, limit(g(x), x, 0),
limit(g(x),x,0,minus),  limit(g(x),x,0,plus)];
wxdraw2d(yrange=[-2,2],color=blue,
 key = "|x| / x",explicit(g(x),x,-2,2));

If a function grows without bound as x approaches the specified value, 
wxMaxima returns the word "infinity" if Maxima does not know whether
the limit is a very large positive value or a very large negative value.

With enough information, wxMaxima prints the symbol for infinity or 
-infinity and below.

[limit(1/x,x,0),limit(1/x, x, 0, minus), 
limit(1/x, x, 0, plus)];

### Taylor and Maclaurin Series

Maxima's  taylor( ) command provides a Taylor series representation 
of an expression. The output below shows the results of expansion 
around x = 0 (Maclaurin series) and x = pi of sin(x)/x.

taylor(sin(x)/x,x,0,4);
T4: taylor(sin(x)/x, x,%pi,4);

These series must be truncated if they are to be used in further 
analysis. taytorat( ) truncates a series. We use the fourth-degree 
polynomial approximation with expansion around pi, T4.

fpprintprec:5$ T4rat:taytorat(T4); 
expand(float((T4rat)));

The graph below shows that the approximation used above tracks the 
actual function well in the neighborhood of pi (3.14...). 

For x > 2*pi or so, however, the approximation becomes 
increasingly bad.

wxdraw2d(yrange=[-2,2],
explicit(sin(x)/x, x,0,4*%pi),color=blue,
key="T4rat",explicit(T4rat,x,0,4*%pi)    )$

## Differential Equations

As with integration, Maxima can analyze differential equations either symbolically 
or numerically. The example below analyzes the behavior of a function that
produces the differential equation specified below. The single quote
before diff( ) instructs Maxima not to implement differentiation. The expression 
named \verb!expr! is equivalent to dx/dt = alpha*x(t).

kill(all)$  expr :'diff(x(t),t)+%alpha*x(t)=0;

The functional notation,  x(t),  must be used when desolve( ) 
is used. We solve this ODE twice. In the first case, the initial 
value of x, x(0), is not specified. In the second case, the value is 
specified. 

sol1 : desolve(expr,x(t)); 
atvalue(x(t),t=0,%pi)$  sol2:desolve(expr,x(t));

Apply the solutions to the expression dx/dt  - alpha*x = 0 
in order to check the solution. Values of 0 are returned, 
confirming the solution.

[diff(rhs(sol1),t ) + %alpha*rhs(sol1),
diff(rhs(sol2),t ) + %alpha*rhs(sol2)];

### Numerical Analysis

Numerical methods can be used when analytical solutioins are not
available (or when they are, but all that is required is a numerical
value). The Runge-Kutta numerical  method, command: rk( )), 
can be applied. That method is applied below to the example above, 
but with a specified value of alpha.
rk( )produces a list of five pairs of output. Each pair identifies
a value of t and the value x that is implied by that value of t.

kill(all)$   dxdt:-0.05*x$  rk(dxdt,x,1,[t,0,4,1]);

Exercise: Confirm that these five values are approximately those 
implied by sol1 given x(0)=1 and alpha=0.05.

The first three commands provide one way to carry out this 
exercise. The matrix( ) command produces a table. The cons( )
command, inside matrix( ) constructs a new list but adding text
to the original lists.

tList:makelist(i,i,0,4); 
x(t):=1*exp(-.05*t); 
xList: x(tList);
matrix(cons("t",tList),cons("x",xList));

### Second-Order ODEs

Maxima can treat higher-order ODEs. The example below shows 
how to analyze a second-order ODE. The assume( ) command 
avoids a dialog with Maxima regarding omega's nature.

kill(all)$                assume(%omega > 0)$
expr : 'diff( x(t), t, 2) + %omega^2*x(t) = 0;
atvalue(  diff(x(t), t, 1), t= 0, v0  )$
sol : desolve (expr, x(t));
ratsimp(diff( rhs(sol), t, 2) + %omega^2*rhs(sol));

### Extensions

No input or out for this section.