Reference.mws

PHYS2030 Computational Methods for Physicists and Engineers

A Maple Primer: this introduction to Maple 7/8/9  provides a tutorial and reference to the question:

how does one do Calculus and Linear Algebra in Maple?

While at some point you should know most of the commands explained by heart, you can use this worksheet during homework and tests.

Note that there are tutorials in Maple's help system, and that each command can be queried by ?commandname, e.g., ?diff

1.1 Basics:

Commands are entered into Maple at the prompt (> symbol), and are terminated by a semicolon or colon. A colon suppresses the display of the result:

>    3+9;

12

To assign the result to a variable:

>    a:=3+9;

a := 12

Want to know what a  is?

>    a;

12

>    b:=4+10:

>    b;

14

>    b-a;

2

>    b*a;

168

>    b^a;

56693912375296

>    b**a;

56693912375296

Both versions work, but we will avoid the Fortran-style **  operator!

At this point we know how to do arithmetic. We can save our session for future reference. Maple stores it as a *.mws  file. Note that when we re-open the file, it will display its content, but it will not be active unless we execute (enter) the commands. This can also be done for an existing worksheet from the Edit  menu in one go.

Let us continue to the level of graphing calculator  proficiency: use the trig and exponential functions and graph them.

>    f1:=sin(2*x);

f1 := sin(2*x)

>    f2:=exp(-3*x);

f2 := exp(-3*x)

f1  and f2  are not functions, but EXPRESSIONS in Maple. We will come to functions/maps later!

The mathematical symbol 'pi' requires to use an upper-case Pi :

>    plot(f1,x=-Pi..Pi);

[Maple Plot]

Now suppose we want to show both functions on one graph, and distinguish them. We will need to restrict the range on the y-axis, or the graph becomes unwieldy.

>    plot([f1,f2],x=-Pi..Pi,-1..2,color=[red,blue]);

[Maple Plot]

Here we used the concept of a list, which is an ordered set: we placed the two expressions to be graphed in a definite order, and the colors were associated correspondingly.

The asymptotic behavior of functions is important to us, mostly because we will need to perform integrals, often over an infinite x-range in upper-year physics.

Exercise 1.1.1:

Combine functions such as exp(sin(x)) and graph them. Observe their properties (asymptotic behavior, periodicity, etc.). Graph functions and their inverses (exp, ln;  sin, arcsin; etc.), and observe the relationships between their graphs.

1.2 Simple Calculus:

To understand the behavior of functions we look at derivatives. Expressions are differentiated using the diff  command.

>    f:=exp(sin(s));

f := exp(sin(s))

>    f1:=diff(f,s);

f1 := cos(s)*exp(sin(s))

>    f2:=diff(f,s$2);

f2 := -sin(s)*exp(sin(s))+cos(s)^2*exp(sin(s))

>    plot([f,f1,f2],s=-2*Pi..2*Pi,color=[red,blue,green]);

[Maple Plot]

Exercise 1.2.1:

Observe the behavior of the function and its derivatives at the critical points (vanishing derivative=stationary point, vanishing 2nd derivative=turning point), by inspecting the graph.

How do we find the critical points? This is a question of computer algebra. The solve  and fsolve  commands will be explored to find the zeroes.

>    sol:=solve(f1,s);

sol := 1/2*Pi

Note that solve found only one root! The function is periodic, and this root should be repeated infinitely many times. In addition there are other roots.

The message: be careful with answers from Maple (or any other computer algebra system).

Let us start with the one root Maple found. We want to verify it:

>    subs(s=sol,f1);

cos(1/2*Pi)*exp(sin(1/2*Pi))

>    simplify(%);

0

We learned the substitute command, and that it does not automatically simplify the result. Alternatively:

>    eval(f1,s=sol);

0

Do we have a minimum or a maximum here?

>    eval(f2,s=sol);

-exp(1)

>    evalf(%);

-2.718281828

We have a negative second derivative, and thus a maximum. (use the parabola x^2/-x^2 to remember the condition!)

Now let us learn a numerical method to arrive at the roots of the derivative function. We use a black-box method. fsolve requires that you bracket the desired root, we use 1 and 2 as brackets to find the above solution in numerical form.

>    sol1:=fsolve(f1,s=1..2);

sol1 := 1.570796327

>    eval(f1,s=sol1);

-.5575287929e-9

Note that the function evaluation does not yield zero. This is caused by the limited precision, which is set by default to 10. Digits is an internal variable which allows us to increase the precision.

>    Digits:=14;

Digits := 14

>    sol1:=fsolve(f1,s=1..2);

sol1 := 1.5707963267949

>    eval(f1,s=sol1);

-.91898820644691e-14

We find that we can push the accuracy closer towards the desired result. In contrast to floating point number system based computing languages we are not limited to precision restrictions such as 8 or 16 or 32. We can set Digits  to almost anything we like (or need).

Nevertheless, the problem remains: when do we recognize a number to be close enough to zero that it may actually be zero?

>    evalf(10^(-Digits));

.10000000000000e-13

Exercise 1.2.2:

Use the graph as a guide to find the other critical points numerically, and verify the minimum/maximum property. Observe what happens when no root can be found (e.g., by trying to solve for a root of the function f  itself).

Now we calculate some anti-derivatives:

>    int(f,s);

int(exp(sin(s)),s)

This is the display when Maple can't do it. It just shows the integral sign. In contrast to derivatives, which can always be obtained mechanically, integrals of elementary functions often cannot be found in terms of elementary functions.

Exercise 1.2.3:

Use known integrals of elementary functions to test Maple's int  engine. It is an engine, because it does not use look-up, but algorithmic approaches, something you learned in 1st year Calculus.

An important fact, however, is that definite integrals can be approximated numerically in these cases where the anti-derivative is not known as a simple function. We use a black-box approach, and an inert version of the integration command:

>    evalf(Int(f,s=0..Pi));

6.2087580357111

This number represents the area under the curve between the vertical lines at s=0  and s=Pi .

Exercise 1.2.4:

Come up with examples where the function crosses the x-axis in the range of integration, and verify that the sign of the function matters, i.e., calculate the areas between a  and x0  and x0  and b  as well as the area from a  to b  and explain! Here a,b  are the boundaries of integration, and x0  is the location of the zero. You can use the sin(x)  function, and bracket only one root.

What else can we do?

>    sum(n^2,n=1..N);

1/3*(N+1)^3-1/2*(N+1)^2+1/6*N+1/6

Comparable to integration: Maple can calculate closed-form sums.

>    a:='a':

>    limit(sin(a*x)/x,x=0);

a

Maple can calculate limits. The statement a:='a':  serves to unassign the variable a , i.e., to turn it back into symbol a . Maple also does one-sided limits

>    f:=abs(x);

f := abs(x)

>    f1:=diff(f,x);

f1 := abs(1,x)

>    limit(f1,x=0);

undefined

>    limit(f1,x=0,left);

-1

>    limit(f1,x=0,right);

1

>    plot([f,f1],x=-2..2,color=[red,blue]);

[Maple Plot]

Exercise 1.2.5:

Explore other functions with cusps, e.g., |cos(x)|, and explore the behavior of the derivative function at the cusp.

Improper integrals

Many integrals in physics involve an infinite range, e.g., we could have an expression for the radiated power of a light source at a given wavelength (frequency). Suppose we want to know the entire radiated intensity. This requires an integral over all wavelengths (or frequencies). Problems with the integrand at the endpoints lead to potentially infinite results. Modern physics (quantum mechanics) was discovered as a result of the demand to avoid problems with a classical calculation (Planck's hypothesis).

>    int(x^2*exp(-x),x=0..infinity);

2

>    plot(x^2*exp(-x),x=0..10);

[Maple Plot]

What about the following generalization:

>    int(x^2*exp(-a*x),x=0..infinity);

limit(-(a^2*x^2*exp(-a*x)+2*a*x*exp(-a*x)+2*exp(-a*x)-2)/a^3,x = infinity)

Why can't Maple do the limit?

>    int(x^2*exp(-a*x),x=0..infinity) assuming a>0;

2/a^3

Exercise 1.2.6:

Choose bell-shaped curves [such as 1/(a^2+x^2) ] and find improper integrals over the entire x-axis. Graph the integrand for some choices of a  and convince yourself why the area under the curve depends on a  as given by the symbolic answer.

1.3 Linear Algebra:

Step 1: call the LinearAlgebra package. Before calling it we reset the Maple session to clear all variables. Call the LinearAlgebra package with a semicolon terminator so that the names of defined procedures come up.

>    restart; with(LinearAlgebra):

The two most important objects defined in the package are Matrix and Vector. We define variables by entering content using lists and lists of lists:

>    V:=Vector([1,2,3]);

V := Vector(%id = 524836)

>    A:=Matrix([[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]]);

A := Matrix(%id = 674404)

>    a11:=1: a12:=5: a13:=7: a21:=0: a22:=1: a23:=2: a31:=9: a32:=6: a33:=1:

>    Determinant(A);

16

>    X:=LinearSolve(A,V);

X := Vector(%id = 747036)

>    whattype(X);

Vector[column]

>    A.X;

Vector(%id = 744168)

>    A;

Matrix(%id = 674404)

>    evalm(A.X - V);

vector([9/2*a11-7*a12+9/2*a13-1, 9/2*a21-7*a22+9/2*a23-2, 9/2*a31-7*a32+9/2*a33-3])

>    evalm(A &* X - V);

vector([0, 0, 0])

We see that Maple9 is not free of problems (at the level of version 9.01). The above two lines should work in the same fashion.

Exercise 1.3.1:

Set up your own inhomogeneous systems of equations of dimensionalities 2,3,4, etc.. Check that the coefficient matrix is non-singular by calculating the determinant, then find the solution, and verify it.

When one needs to solve many systems with the same coefficient matrix, but different right-hand sides, it may be opportune to calculate the inverse of the matrix:

>    Ai:=MatrixInverse(A);

Ai := Matrix(%id = 780352)

>    Ai . A;

Matrix(%id = 853144)

yikes, not again!

>    evalm(Ai &* A);

matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

Let us explore a little the depths of symbolic computation scoping rules. The matrix A was defined from symbolic entries, and then the entries were assigned to be numbers. We can unassign them and watch what happens.

>    a11:='a11': a12:='a12': a13:='a13': a21:='a21': a22:='a22': a23:='a23': a31:='a31': a32:='a32': a33:='a33':

>    A;

Matrix(%id = 674404)

We are able to turn back to the original general result. Not so for the matrix inverse, because it was computed based on the knowledge of the matrix entries:

>    Ai;

Matrix(%id = 780352)

The matrix inverse can be calculated symbolically:

>    MatrixInverse(A);

Matrix(%id = 898944)
Matrix(%id = 898944)
Matrix(%id = 898944)
Matrix(%id = 898944)
Matrix(%id = 898944)
Matrix(%id = 898944)

An important message:

The symbolic matrix inverse requires a substantial amount of storage, which does not arise when the entries are numbers. The number of terms grows to the tune of n!  for an n-by-n matrix.

Suppose that we want to try out how far Maple can go - it would be clumsy to define the matrix by hand. We will learn a sequence command that allows us to automate things:

>    seq(i^2,i=1..10);

1, 4, 9, 16, 25, 36, 49, 64, 81, 100

To assemble a matrix with symbolic entries we also need to know how to concatenate symbols:

>    m:=8;

m := 8

>    a||m;

a8

Let us pick the matrix size:

>    N:=5;

N := 5

>    [seq(a1||j,j=1..N)];

[a11, a12, a13, a14, a15]

That is neat. Now let us do the same with the other index and put it all together:

>    A:=Matrix([seq([seq(a||i||j,j=1..N)],i=1..N)]);

A := Matrix(%id = 1046088)

Wow! The line below writes screens full when terminated with a semicolon:

>    Ai:=MatrixInverse(A):

You can exhaust your computational resources easily by increasing N . While the numerical inversion of a non-singular 20-by-20 matrix is trivial (within Maple or otherwise), the symbolic inversion of a general 9-by-9 matrix will make even the most recent computers chew up their memory.

Inhomogeneous systems of linear equations are relatively straightforward to deal with. Let us proceed with homogeneous systems of equations that arise most often in the eigenvalue-eigenvector problem context. Let us go back to a numerical matrix:

>    A:=Matrix([[1,2,3],[1,4,9],[0,1,0]]);

A := Matrix(%id = 1208852)

>    Determinant(A);

-6

>    V:=Vector([1,1,1]);

V := Vector(%id = 1233796)

We chose a vector which points along the diagonal in Euclidean 3-space R3. When we multiply this vector with A, we obtain a new vector with different length (Euclidean vector norm), and which points in a different direction:

>    V1:=A . V;

V1 := Vector(%id = 20505596)

>    VectorNorm(V);

1

>    VectorNorm(V1);

14

This is clearly the maximum norm, and not the Euclidean norm!

>    VectorNorm(V,Euclidean),VectorNorm(V1,Euclidean);

3^(1/2), 233^(1/2)

This is clumsy. Let us define our own procedure to compute the norm. First we need a way to measure the number of components in our vectors, then we use a simple add function that works in analogy to the sequence construct. Finally we need to know that procedures always return the last computed entry:

>    nops(V);

3

>    VNorm:=proc(W::Vector): sqrt(add(W[i]^2,i=1..nops(W))); end:

>    VNorm(V);

3^(1/2)

>    VNorm([1,2,3]);

Error, invalid input: VNorm expects its 1st argument, W, to be of type Vector, but received [1, 2, 3]

This is the result of our restriction for the dummy argument W  to be of type Vector , it won't work on a list.

Interestingly: we have been able to use loop constructs without defining loops. This is accomplished using the seq, add, mul  functions. We will get to explicit loops later.

Now the problem of eigenvectors: for each N-by-N matrix A  we can find at most N  distinct vector directions such that the result of A.V  will point in the same direction as V . Such vectors are called eigenvectors of A , and the problem is how to find them.

The problem is usually formulated as

A.V = lambda V    ,  where lambda  is an element of the number field; it is a placeholder for a scale factor.

Read this equation carefully: it says that for a given matrix A we are asking to find vectors V and scale factors lambda, such that the action of A onto these vectors corresponds to a stretch with factor lambda.

We can manipulate the right-hand side of the equation by squeezing in an indentity matrix, and moving it over to the left:

>    C:=IdentityMatrix(3);

C := Matrix(%id = 20562260)

>    B:=A-lambda*C;

B := Matrix(%id = 20626396)

Now we have the problem

B.V = 0   (where 0 stands for a zero-vector)

and we are left with the problem to solve a homogeneous system of equations.

We use the placeholder lambda  to our advantage. We recall the result of an important theorem from linear algebra.

A system of equations has a unique solution if the determinant of the coefficient matrix does not vanish, i.e., if the coefficient matrix is non-singular.

In the context of the homogeneous system of equations this means that the only solution acceptable in this case is the trivial solution where all components V_i=0 .

Interesting solutions can be found only if the determinant of the coefficient matrix B  vanishes, i.e., when B  is singular. We use the freedom of the undetermined lambda  to gurantee such solutions,

>    cp:=Determinant(B);

cp := 7*lambda+5*lambda^2-6-lambda^3

The values of lambda  that allow interesting solution vectors V  are obtained from the roots of this polynomial (called the characteristic p.).

An analysis of the cubic (for N=3 ) shows that there can be three real roots, or one real root and a complex conjugate pair. A limiting case is when the minimum in the parabola touches down on the lambda -axis, and the pair becomes a repeated root.

>    plot(cp,lambda=-4..7);

[Maple Plot]

In our example we have three real roots, but note that the entries in the original matrix A  could have been shifted a little to push the minimum onto the lambda -axis, and above.

The roots are called the eigenvalues. For each eigenvalue a solution vector V  can be found, and it is undetermined in its length. First find the roots:

>    Evals:=solve(cp);

Evals := 6, 1/2*5^(1/2)-1/2, -1/2*5^(1/2)-1/2

>    evalf(Evals);

6., .6180339880, -1.618033988

>    Eigenvalues(A);

Vector(%id = 20671772)

Now calculate the eigenvectors one by one:

>    B1:=eval(B,lambda=Evals[1]);

B1 := Matrix(%id = 20723444)

>    V1:=LinearSolve(B1,Vector([0,0,0]));

V1 := Vector(%id = 20790364)

The eigenvector associated with the first eigenvalue is determined up to an arbitrary scale parameter _t0[3].

>    V1:=simplify(1/VNorm(V1) * V1) assuming _t0[3]>0;

V1 := Vector(%id = 538488)

>    VNorm(V1);

1

We have eliminated the freedom by normalizing the eigenvector to unit Euclidean norm, i.e., unit length in real 3-space.

Now calculate the second eigenvector:

>    B2:=eval(B,lambda=Evals[2]);

B2 := Matrix(%id = 1248120)

>    V2:=LinearSolve(B2,Vector([0,0,0]));

V2 := Vector(%id = 666828)

>    V2:=simplify(1/VNorm(V2) * V2) assuming _t0[3]>0;

V2 := Vector(%id = 659104)

This looks messy, but is exact. To get an idea of the direction evaluate to floating point:

>    evalf([V1,V2]);

[Vector(%id = 652164), Vector(%id = 743520)]

Now

>    simplify(DotProduct(V1,V2));

-1/92*46^(1/2)*(9*5^(1/2)+37)/(64+27*5^(1/2))^(1/2)

>    evalf(%);

-.3776155210

>    evalf(add(V1[i]*V2[i],i=1..3));

-.3776155208

The two vectors are not collinear, i.e., they are linearly independent. They are not orthogonal.

Exercise 1.3.2:

Find the third eigenvector of the matrix A  and show that it is linearly independent from the other two.

If the matrix is real, symmetric, and has distinct eigenvalues, then the eigenvectors will turn out to be perpendicular. They are always guaranteed to be linearly independent, and that is the most important feature. One can always apply an orthonormalization procedure (such as Gram-Schmidt) to turn them into an orthonormal basis for the vector space in question.

This latter property makes the eigenvector-eigenvalue problem an important cornerstone for function space as used in quantum mechanics.

The problem is so ubiquitous that we need a shortcut to obtain the complete solution. In symbolic mode the eigenvectors are returned normalized with the maximum norm:

>    Eigenvectors(A);

Vector(%id = 737968), Matrix(%id = 995180)

If we evaluate the matrix to floating-point entries first, the Eigenvector  function applies a numerical diagonalization procedure which returns eigenvectors normalized with the Euclidean norm.

>    Evecs:=Eigenvectors(evalf(A));

Evecs := Vector(%id = 702916), Matrix(%id = 1111004)

The appearance of the placeholder for imaginary parts appears annoying. Let us experiment a little with complex numbers, before we proceed.

>    eq:=z^2+1=0;

eq := z^2+1 = 0

>    solve(eq,z);

I, -I

The imaginary unit is displayed and entered with the upper-case i = I , unlike in textbook math notation. This is done so that we can use i  also for an index, or other occasions.

>    z1:=2*a+3*b*I;

z1 := 2*a+3*I*b

>    Re(z1),Im(z1),conjugate(z1);

Re(2*a+3*I*b), Im(2*a+3*I*b), conjugate(2*a+3*I*b)

>    evalc([%]);

[2*a, 3*b, 2*a-3*I*b]

Note that evalc  cannot operate on the sequence of expressions, but it can operate on the list formed from these. This is so because according to the syntax, appending objects by commas (a sequence) is tantamount to passing several arguments. Putting the brackets around the expression sequence turns it into a single object, a list which is an ordered set. Also note that the functions Re, Im, conjugate  do not evaluate without a call to evalc .

Now let us demonstrate the following. We extract the eigenvectors from Evecs  in the form of a matrix made up of the eigenvectors as columns.

>    Evecs[2];

Matrix(%id = 1111004)

>    whattype(%);

Matrix

>    S:=Re(Evecs[2]);

Error, invalid input: simpl/Re expects its 1st argument, x, to be of type algebraic, but received Matrix(3, 3, [[...],[...],[...]], datatype = complex[8], storage = rectangular, order = Fortran_order, shape = [])

yikes, this could have worked; evalf worked in this way after all!

>    S:=map(Re,Evecs[2]);

S := Matrix(%id = 1181332)

The matrix is not orthogonal. The columns (or rows) are normalized, but not mutually orthogonal. Therefore, the inverse is NOT obtained by transposing it.

>    Si:=Transpose(S);

Error, attempting to assign to `Si` which is protected

OK, there are predefined names that we shouldn't use.

>    Str:=Transpose(S);

Str := Matrix(%id = 1237112)

>    S . Str;

Matrix(%id = 20402396)

>    Str . S;

Matrix(%id = 20475212)

>    Sinv:=MatrixInverse(S);

Sinv := Matrix(%id = 20625312)

>    Sinv . S;

Matrix(%id = 20646344)

This is as close as we can get to a unit matrix when working to 15 digits.

>    S . Sinv;

Matrix(%id = 786312)

The point of the exercise is to show that the original matrix A  is transformed to diagonal form by the transformation matrix made up from the eigenvectors. The order matters, namely:

>    S . A . Sinv;

Matrix(%id = 840352)

>    Sinv . A . S;

Matrix(%id = 1269224)

Therefore, the problem of finding the eigenvectors and eigenvalues is called also the diagonalization problem.

The physics message:

Suppose the matrix A  represents some physics variable in an arbitrary reference frame. This could be the moment of inertia tensor of some 3d object (the matrix A  would be real and symmetric based on the definitions of the principal and deviation moments).

Finding the eigenvectors is tantamount to finding coordinate axes such that when these axes are chosen the tensor is represented by a very simple matrix, namely the matrix of eigenvalues on the diagonal. In classical mechanics this proves that for any shape (a potato!) three axes can be found such that the object can be spun without wobbling. When one chooses an eigenvector of the moments of inertia matrix as rotation axis, one can spin the object smoothly.

Your car mechanic knows about this in practical terms when balancing the wheels of your car: a small weight is added (of a mass and to a location as specified by a sophisticated machine) to ensure that the axle direction becomes prinicipal axis for the wheel. With this technique imperfections in the rubber, as well as dents to the rim can be compensated for.

Exercise 1.3.3:

Use simple matrices A  that are real and symmetric and find their eigenvectors and eigenvalues. Work with the symbolic form only if the eigenvalues and eigenvectors are very simple. Verify that the matrix inverse is given by matrix transpose in this case, and verify the diagonalization calculation.

Quadratic forms: (in 2-space, can be extended to 3-space)

>    A:=Matrix([[1,2],[2,1]]);

A := Matrix(%id = 455680)

>    QF:=Transpose(Vector([x,y])) . A . Vector([x,y]);

QF := x*(x+2*y)+y*(2*x+y)

Question: why did we need the transpose on the left?

By default Vector  generates a column vector!

We choose a scale for the geometric shape defined by the quadratic form:

>    eq:=simplify(QF)=1;

eq := x^2+4*x*y+y^2 = 1

To graph the equation one can proceed in various ways. In the important plots  package there is a convenient function for this case:

>    with(plots):

Warning, the name changecoords has been redefined

>    implicitplot(eq,x=-2..2,y=-2..2,scaling=constrained);

[Maple Plot]

>    Eigenvalues(A);

Vector(%id = 1092436)

A hyperbola is characterized by having one positive and one negative eigenvalue. Let us perform the transformation to normal form:

>    S:=Eigenvectors(A)[2];

S := Matrix(%id = 975644)

>    Sinv:=MatrixInverse(S);

Sinv := Matrix(%id = 511864)

>    Adiag:=Sinv . A . S;

Adiag := Matrix(%id = 725320)

>    QFdiag:=Transpose(Vector([x,y])) . Adiag . Vector([x,y]);

QFdiag := -x^2+3*y^2

>    eqdiag:=simplify(QFdiag)=1;

eqdiag := -x^2+3*y^2 = 1

>    implicitplot(eqdiag,x=-2..2,y=-2..2,scaling=constrained);

[Maple Plot]

Clearly, when we choose to transform the object to principal form first, no cross terms appear. It is of interest to superimpose the eigenvectors onto the original graph:

>    P1:=implicitplot(eq,x=-2..2,y=-2..2,scaling=constrained,color=red):

>    line1:=S[1,1]*x+S[2,1]*y=0;

line1 := -x+y = 0

>    line2:=S[1,2]*x+S[2,2]*y=0;

line2 := x+y = 0

>    P2:=implicitplot([line1,line2],x=-2..2,y=-2..2,scaling=constrained,color=[blue,green]):

>    display(P1,P2);

[Maple Plot]

Exercise 1.3.4:

Repeat the above procedure for a matrix with two positive eigenvalues. When calling the implicitplot  it may be needed to adjust the range for the x and y variables so that the entire shape becomes visible. The interpretation of the eigenvalues (their inverses) should become quite obvious in this case.

1.4 Differential equations:

>    restart;

Let us explore a few simple equations: relationships between a function and some of its derivatives define an equation that determines the function. For example: the exponential function is defined as the function which is equal to its derivative.

>    DE:=diff(y(x),x)=y(x);

DE := diff(y(x),x) = y(x)

>    sol:=dsolve(DE);

sol := y(x) = _C1*exp(x)

>    rhs(sol);

_C1*exp(x)

A special solution is obtained by fixing an initial condition. Let us generalize slightly: we formulate a problem where the rate of decay is proportional to the amount of material present. A positive material constant a  controls the proportionality.

>    DE:=diff(N(t),t)=-a*N(t);

DE := diff(N(t),t) = -a*N(t)

Suppose that at time zero we have 100 atoms.

>    IC:=N(0)=100;

IC := N(0) = 100

Equation and initial condition are placed into a set:

>    sol:=dsolve({DE,IC});

sol := N(t) = 100*exp(-a*t)

The constant a  has dimensions of inverse time. Note that functions can only be evaluated for dimensionless argument! The constant t0=1/a  thus has something to do with the so-called half-life. It describes the time scale at which a population shrinks to 1/e  of its initial value. We simply measure time in units of this scale for the following graph:

>    plot(subs(a=1,rhs(sol)),t=0..3,N=0..100);

[Maple Plot]

Exercise 1.4.1:

Define a pair of first-order decay equations, where one population N1(t)  decays with a given rate 1/t0 , and feeds a second population N2(t)  with this same rate which itself decays at a rate of 0.5/t0 .  Solve the pair of equations and graph their solutions. Explore what happens when the daughter population N2  decays with twice the rate of the parent.

Define the harmonic oscillator equation for a frictionless oscillator.

>    restart;

>    DE:=m*diff(x(t),t$2)=-k*x(t);

DE := m*diff(x(t),`$`(t,2)) = -k*x(t)

A second-order ordinary differential equation requires us to specify two initial conditions. These correspond to specifying intial displacement from equilibrium, as well as intial speed.

We choose a displacement of one unit (would be 1 meter in SI units), and no initial speed.

The formulation of this initial condition is somewhat tricky:

>    IC:=x(0)=1,D(x)(0)=0;

IC := x(0) = 1, D(x)(0) = 0

>    sol:=dsolve({DE,IC});

sol := x(t) = cos(1/m^(1/2)*k^(1/2)*t)

We recognize the circular frequency omega=sqrt(k/m) . The period is given as T=2*Pi/omega , and we measure time in units where T=1 .

>    sol1:=subs(k=1,m=1,rhs(sol));

sol1 := cos(t)

>    plot([sol1,diff(sol1,t)],t=0..3*Pi,color=[red,blue],title="position(red) and velocity(blue) in a harmonic oscillator");

[Maple Plot]

Exercise 1.4.2:

Modify the initial condition to consider the following cases:  (a) the mass starts at equilibrium with a positive velocity; (b) the mass starts as above, but with a negative velocity.

Finally, we show how to use dsolve  in numeric mode. A nonlinear equation is one where powers or other functions of the dependent variable (or its derivatives) appear in the equation. The full-blown mathematical pendulum is such an example. The variable is the rotation angle in this case.

>    DE:=diff(theta(t),t$2)=-g/l*sin(theta(t));

DE := diff(theta(t),`$`(t,2)) = -g/l*sin(theta(t))

Suppose we are in SI units, and approximately g=10 m/sec^2, and the pendulum arm is of 1 m length:

>    g:=10: l:=1:

>    DE;

diff(theta(t),`$`(t,2)) = -10*sin(theta(t))

Start the pendulum near the top from rest and keep in mind that the angle is measured in radians:

>    IC:=theta(0)=3.13,D(theta)(0)=0;

IC := theta(0) = 3.13, D(theta)(0) = 0

>    sol:=dsolve({DE,IC},numeric,output=listprocedure):

>    Theta:=eval(theta(t),sol):

>    Theta(1);

3.00447005752296592

>    plot(Theta(t),t=0..10);

[Maple Plot]

We can also define the derivative of the solution, given that we have a second-order equation:

>    Omega:=eval(diff(theta(t),t),sol):

>    plot([Theta(t),Omega(t)],t=0..10,color=[red,blue]);

[Maple Plot]

Exercise 1.4.3:

Explore the mathematical pendulum solution for different initial conditions: (a) a smaller intial displacement: observe how the period of oscillation depends on the intial angle, except when it is small enough such that the linearized regime is reached when sin(x(t)) can be replaced by x(t); (b) the same initial displacement as above, but with additional initial angular speed.

Note that inaccuracies in the numerical solution may occur for intial angles rather close to Pi .

Some explanations are in order to understand how the initial condition on the velocity is specified. We have worked so far only with expressions, and differentiated them. The subtle point about the initial condition is that first we have to differentiate the map x(t), and then evaluate it at t=0, the initial point.

Maps (functions, procedures) in Maple:

If we would like the convenience of true maps that we enjoy when calling sin(s), exp(t), etc., from built-in functions, we need to define our own maps. This is done in a way that reminds us of the formal mathematical definition of a map. The arrow is formed using the hyphen and 'greater than' signs.

>    f:=x->x^2*sin(x);

f := proc (x) options operator, arrow; x^2*sin(x) end proc

>    f(s);

s^2*sin(s)

>    f(x);

x^2*sin(x)

>    f(2);

4*sin(2)

Maps are differentiated with the D  operator:

>    g:=D(f);

g := proc (x) options operator, arrow; 2*x*sin(x)+x^2*cos(x) end proc

>    g(s);

2*s*sin(s)+s^2*cos(s)

>    diff(f(s),s);

2*s*sin(s)+s^2*cos(s)

The conversion of a map into an expression has been done above by calling the map with argument x  or some other independent variable name. f(x), g(x)  have become expressions in x .

To convert the other way around one uses the unapply  procedure. For example:

>    h_expr:=diff(f(t),t$2);

h_expr := 2*sin(t)+4*t*cos(t)-t^2*sin(t)

>    h:=unapply(h_expr,t);

h := proc (t) options operator, arrow; 2*sin(t)+4*t*cos(t)-t^2*sin(t) end proc

Maps can be graphed directly, but one is not allowed to enter a name for the independent variable.

>    plot([f,g,h],-2..3,color=[red,blue,green]);

[Maple Plot]

1.5 Programming in Maple

We complete our reference worksheet/tutorial by providing an introduction to the most important loop elements and decision-making syntax elements.

We have introduced simple sequencing elements which allowed us to avoid writing explicit do loops. Now we introduce multi-line do loops, and ultimately also multi-line procedures.

To have an example to work on we consider the Taylor series expansion of a known function. We begin by using the Maple built-in procedure, and then construct the Taylor polynomail ourselvcs.

>    restart;

>    taylor(sin(x),x=0,8);

series(1*x-1/6*x^3+1/120*x^5-1/5040*x^7+O(x^8),x,8)

In order to work with the result we have to subtract the order term:

>    sinT8:=convert(taylor(sin(x),x=0,8),polynom);

sinT8 := x-1/6*x^3+1/120*x^5-1/5040*x^7

>    plot([sin(x),sinT8],x=0..2*Pi,-1.5..1.5,color=[red,blue]);

[Maple Plot]

How are the terms generated? We choose a maximum order N ; we need the derivatives of the function and evaluate it at the expansion point, namely the origin. For simplicity we work with expressions at first, and only later with maps.

>    N:=7;

N := 7

>    f:=sin(x);

f := sin(x)

>    fT:=eval(f,x=0): for n from 1 to N do: fd:=diff(f,x$n); fT:= fT + eval(fd,x=0)*x^n/n!; od;

fd := cos(x)

fT := x

fd := -sin(x)

fT := x

fd := -cos(x)

fT := x-1/6*x^3

fd := sin(x)

fT := x-1/6*x^3

fd := cos(x)

fT := x-1/6*x^3+1/120*x^5

fd := -sin(x)

fT := x-1/6*x^3+1/120*x^5

fd := -cos(x)

fT := x-1/6*x^3+1/120*x^5-1/5040*x^7

One can suppress the output during the do-loop by terminating the loop with od:  - that is the usually preferred mode, and then one can put print  statements into the loop for debugging.

The loop can also be terminated with the end do:  statement. We will stick, however, with the historical terminator that represents an inverted do .

This loop is simple enough that we could have programmed a one-liner:

>    eval(f,x=0)+add(eval(diff(f,x$n),x=0)*x^n/n!,n=1..N);

x-1/6*x^3+1/120*x^5-1/5040*x^7

Keep in mind, however, that the explicit do loop is easier to read and debug.

Exercise 1.5.1:

Test the loop with different functions and compare your results with Maple's result from the taylor  procedure. Change the order of expansion. Change the expansion point.

Note that we have separated out the first term in the Taylor polynomial, because:

>    add(eval(diff(f,x$n),x=0)*x^n/n!,n=0..N);

Error, wrong number (or type) of parameters in function diff

We can incorporate the first term by avoiding the call to diff  using an if-then-else-fi  branch construct:

>    for n from 0 to N do: if n=0 then fT:=eval(f,x=0); else fd:=diff(f,x$n); fT:= fT + eval(fd,x=0)*x^n/n!; fi; od;

>    fT;

x-1/6*x^3+1/120*x^5-1/5040*x^7

The branching causes the output to be suppressed (that may not have been intentional by the software developers).

An important comment concerns the recursive updating  of the expression fT :

The only reason why a statement like fT:= fT + blah;   works is that fT  as been initialized to some value (expression) before the first appearance of the statement. Otherwise an infinite loop might be generated, although Maple usually catches attempts of this sort:

>    y1:=y1+5;

Error, recursive assignment

The above statement has no meaning if y1  has not been previously declared. Undeclared symbols are not initialized (Fortran would initialize a previously unused variable to zero). We can't add anything to an unitialized symbolic variable.

There are more loop programming elements that control the flow, such as while  - it is well worth looking at the examples in the help pages for this type of control.

Let us now implement a procedure that uses the differentiation operator, i.e., it works on maps to achieve the same goal. It will be an example that shows how to pass arguments into a procedure. Note that one cannot pass answers back from the procedure, except via the final executed statement. If a procedure needs to update several variables, it is done through a global  declaration of these.

The differentiation operator D  is applied repeatedly using the D@@n  syntax. Consult the help page ?D  and look at the examples!

>    MyTaylor:=proc(fmap,point,maxorder::nonnegint) local order,res,x;

>    res:=fmap(point);

>    for order from 1 to maxorder do:

>    res:=res + (D@@order)(fmap)(point)*(x-point)^order/order!;

>    od;

>    unapply(res,x);

>    end:

A note on writing (typing) procedures:

as a beginner keep to writing one statement per command prompt, unless they are very straightforward;

write the statements at isolated prompts (execution groups), and do not Enter the statements;

after finishsing with the end statement review the program flow; check for uncompleted do -loops, and if-then-else-fi  branching constructs.

join the execution statements (prompts) into a single execution group for the entire procedure (from proc  to end: );

now submit the procedure to the Maple engine; if you terminate it with end;  the syntax will be displayed as understood by Maple.

>    MyTaylor(exp,2,5);

proc (x) options operator, arrow; exp(2)+exp(2)*(x-2)+1/2*exp(2)*(x-2)^2+1/6*exp(2)*(x-2)^3+1/24*exp(2)*(x-2)^4+1/120*exp(2)*(x-2)^5 end proc

>    unapply(simplify(taylor(exp(x),x=2,6)),x);

proc (x) options operator, arrow; series(exp(2)+exp(2)*(x-2)+1/2*exp(2)*(x-2)^2+1/6*exp(2)*(x-2)^3+1/24*exp(2)*(x-2)^4+1/120*exp(2)*(x-2)^5+O((x-2)^6),x-2,6) end proc

>    cosT6:=MyTaylor(cos,0,6);

cosT6 := proc (x) options operator, arrow; 1-1/2*x^2+1/24*x^4-1/720*x^6 end proc

>    plot([cos(x),cosT6(x)],x=0..2*Pi,-1.5..1.5,color=[red,blue]);

[Maple Plot]

Keep in mind that MyTaylor  returns a procedure (a map), which means that it will be plotted in the same way as the Maple built-in functions. To plot without converting to expressions change the specification of the independent variable range:

>    plot([cos,MyTaylor(cos,Pi,6)],0..2*Pi,-1.5..1.5,color=[red,blue]);

[Maple Plot]

Exercise 1.5.2:

Explore local Taylor expansions by making use of the defined procedure MyTaylor . In particular, consider effects from a finite radius of convergence by looking at the function 1/(1+x^2)  or similar expressions while expanding about x=0 . Generate expansions about other points and graph your results.

Exercise 1.5.3:

Design a procedure of your own which uses Maple's taylor  procedure (or our own MyTaylor ) to investigate the convergence behavior of the Taylor series about x=0  (McLaurin series) of a function provided by the user up to some fixed order N . The last statement executed by the procedure should be a plot which displays the function together with the polynomials of order n=1..N  in different colors.

>