Chapter 2

Introduction to Maple

2.1 Introduction

In the course of the following chapters, we will derive and solve some of the most fundamental equations in microfluidics and fluid physics. It seems reasonable to use a suitable Algebra package, e.g., for visualizing the results and solutions to the equations found. In many ways, this allows us to get a better feel for the effects that these equations describe.

In this section, we will briefly look at some of the most important functions that Maple provides for us. If the reader already has some understanding of how these Algebra packages work, it may be worth skipping this section. For readers entirely inexperienced with these tools, it may be useful to refer to one of the many great introductory courses and textbooks on Maple. There are also many useful online resources, e.g.,www.mapleprimes.com. A very good textbook is the Maple Advanced Programming Guide [1] which covers the fundamentals as well as some of the very fine details of Maple.

Version. All examples used in this book have been made with Maple 18 using an academic license. Most of the examples will work fine on previous versions of Maple although some details, e.g., the plotting interface, have seen some major changes in recent Maple versions and may therefore, not work directly in older version. Please refer to the online documentation for a check on the feature compatibility.

2.2 Elementary Maple Commands

In this section, we will discuss some of the most fundamental operations in Maple. Feel free to consult the online help files and documentations for more information on the respective functions.

Command Termination. In Maple all commands are terminated with a semicolon “;”. The following code will assign the value 4 to the variable a

a:=4;

Omitting the semicolon results in an error

Warning, inserted missing semicolon at end of statement

Maple will usually output a line of code indicating that the command has been successfully executed. These lines usually only repeat the command entered. If a command is terminated by a colon “:” this output is omitted. This is often referred to as the silent mode. Please note that the silent mode will output an error if the command inserted is wrong.

Assignment. We have already performed an assignment without being aware of it. Assignments in Maple are made by usage of “:=”. Once we made the following assignment

a:=4;

we will obtain the assigned value 4 whenever we query a. This means that by inserting

a;

Maple will output

4

It is important to note that the mere equal sign “=” indicates an equation in Maple. Assignments must be made by using “:=” and not just “=”. It is common to miss this especially when working with Maple for the first time.

Equations. Equations in Maple are simply written as

x^2+2*x=4;

The left-hand side of this equation can be extracted using the function lhs(...), whereas the right-hand side of the equation can be extracted using rhs(...). It is possible to assign equations to variables, e.g.,

equl:=x^2+2*x=4;

solve. The function solve is one of the most important functions in Maple. As the name suggests, this function allows us to solve an equation. As an example, by typing

solve(4*x^2-2=0,x);

Maple would return

(1/2)*sqrt(2), -(1/2)*sqrt(2)

which are obviously the two correct results.

eval and evalf. The Maple function eval is used to force the evaluation of an expression using, e.g., a given value for a variable. The function evalf forces the evaluation to a floating point value. As an example

poly:=x^3+3*x+2;
eval(poly,x=2);

will return the value 16. Using evalf we can force the evaluation of π by typing

evalf(Pi);

will return the value 3.141592654.

Functions. Functions are defined much the same way as assignments are made. If we want to create a function f (x) = x2 we type

f:=x->x^2;

This creates a function f which is dependent on the independent variable x. We can now type, e.g.,

f(3);

whereupon we would obtain the correct result, i.e., the value 9. Obviously, we can reassign the result of a function call to a different variable, e.g.,

c:=f(3);

where the value 9 would be assigned to the variable c.

A function can also be defined using unapply(expression,variable). This will turn expression into a function of variable. As an example, we could have defined f by typing

f:=unapply(x^2,x);

which would have yielded the same function as defined previously.

Built-In Functions.Maple has wide range of built-in functions which we will not discuss in detail here. These functions include trigonometric, exponential, Bessel functions, and so forth. Obviously, Maple understands fundamental calculus such as addition and multiplication.

Vectors and Matrices. In Maple, vectors are created using the function Vector. This means thatb:=Vector([4,4,4,4,4,4,4,4,4,4]);

will create a vector called bsi1_e with 10 entries, each of which carries the value 4. In Maple, the square brackets [...] are used to denote an ordered set of values.

Likewise, matrices are constructed using the function Matrix. This means that

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

will create a tridiagonal matrix A. As you can see, a matrix is nothing other than an ordered set of vectors, i.e., an ordered set of ordered sets of values.

Values from vectors and matrices are accessed using indices. The index is one-based in Maple, which means that b[1] would return the first entry in the vector b. The first entry in the second line of A is queried using A[2,1]. Here, the first index is the line index, whereas the second index is the column index. This is due to the fact that the matrix A is essentially a vector filled with vectors.

The number of elements within a vector can be obtained using the function numelems. Please refer to listing 26.4 for an example of this function.

Comments. In Maple, comments are denoted with the hashtag “#”. Every line starting with a hashtag will be ignored by Maple.

The Function restart. We have already seen that we can assign values to variables. Variables can be reassigned at any given time. In longer worksheets, one may be confused as to which variables have already been assigned values and which ones have not. In order to start a calculation from scratch it is not necessary to close and reload the worksheet. The command restart; will result in a restart of the Maple server. At this point, all prior definitions are erased and all variables are cleared. It is a good idea to call this function at the top of every worksheet. This ensures that whenever the worksheet is restarted, there are no leftover assignments that may potentially cause trouble.

The Function read. We will use the function read(...) quite frequently. read takes one argument which is a file name. It will simply paste the content of this file at the respective position in the Maple worksheet. The file is supposed to contain Maple code. If the file contains definitions, these definitions will become available within the Maple worksheet. Throughout this book, we will use the file Core.txt which is explained in detail in section 2.3. This file contains a number of useful function definitions which allow us to make the Maple code within the actual worksheet significantly shorter and more readable.

Plotting.Maple has a very advanced plotting interface. It is very easy to quickly plot a function, and thus get a good idea about an obtained solution. Take the function f (x) we just defined, we could plot this function using the following code

plot(x->f(x),0..2);

The command plot invokes the plotting interfaces. The expression x->f(x) indicates that the independent variable is supposed to be x and the values of f (x) are obtained. The plotting interval is indicated by the range 0..2. The output of the plotting command is shown in Fig. 2.1.

f02-01-9781455731411
Fig. 2.1 The function f (x) = x2as plotted via the Maple command plot.

We will rarely use the plotting interface directly, as we will define a wrapper function called quickplot which is explained in detail in section 2.3.

The Function typeset. The function typeset in Maple is used to create a string that will be formatted mathematically correct. This function is really only used when, e.g., axis labels are supposed to be provided in string form (as is the norm in the Maple plotting interface), but should be typeset mathematically correct. An example would be the expression a2si6_e which could be converted into a mathematically correctly printed string by typing

s:=typeset(a/2);

If the string s is passed as a parameter to a plotting function, the resulting fraction would be correctly typeset.

The Function with.Maple has a wide set of packages and tools that can be used, depending on the application at hand. However, not all of these packages are loaded at startup. This is mainly due to the fact that some functions are different depending on the context in which they are used. Distribution functions may be required for statistical tasks, but they may also exist in physics packages and could be implemented with slightly different syntax.

In order to access functions defined in a given module, the function with(...) is used. The package to load is simply indicated in the brackets.

The Function display. The function display is used to combine two plots into one. Extending the example, we could combine the plot of f (x) and 2 · f (x) into one plot by writing

display(plot(x->f(x),0..2),plot(x->2*f(x),0..2));

Procedures. In Maple, procedures can be thought of as subroutines which take a given number of arguments, process these arguments, and (potentially) return a value. As an example, take the function AdamsMoulton2 which is defined in listing 27.5

AdamsMoulton2:=proc(df,_from,_to,F0,points) local ret,delta,i,Fpred:
  ret:=Matrix(points,2):
  delta:=(_to-_from)/(points-1):
  ret[1,1]:=_from:ret[1,2]:=F0:
  for i from 2 to points do:
 ret[i,1]:=evalf(ret[i-1,1]+delta):
 Fpred:=ret[i-1,2]+delta*df(ret[i-1,1],ret[i-1,2]):
 ret[i,2]:=evalf(ret[i-1,2]+delta/2*(df(ret[i,1],Fpred)+df(ret[i-1,1],ret[i-1,2]))):
  end do:
 return ret;
end proc:

Here we define a procedure indicated by the Maple primitive proc. The procedure is expected to have the arguments df, _from, _to, F0, and points. It would be possible to indicate the type of values expected. As an example, if points is expected to be an integer, we could write points::integer. We could even indicate a default value in the form of points::integer:=100 where Maple would assume a value of 100 if the user fails to supply a value.

After the procedure definition we can list a number of local variables which are denoted by the primitive local. These are variables which only exist inside of the procedure and cannot be accessed outside of it. It is possible not to declare these variables, but to simply use them in the code. However, Maple will issue a warning of the type

Warning, `ret` is implicitly declared local to procedure `AdamsMoulton2`

if we fail to declare the variable ret as being local. The definition of the procedure ends with the primitive end proc. If the procedure should return a value, the primitive return followed by the value must be indicated. In the above example, the procedure returns the value of the local variable ret.

The Functionsprintf. Especially inside of procedures, it may sometimes be necessary to output values or a message to the user. This can be done using the command sprintf(...). This function has a very powerful interface for formatted output. Please check the online documentation for more details.

if/then/elif/end if.Maple handles branching using if/then/elif/end if constructions. As an example, the variable a can be tested using

if a < 1 then
 #code segment 1
elif a < 5 then
 #code segment 2
else
 #code segment 3
end if

If a < 1 the first code segment is executed. If a < 5 the second code segment is executed. In all other cases, the third code segment is executed. These types of branching are very useful, especially inside of procedures.

for/next.Maple also handles for/next constructs which are known from other programming languages. As an example, the following code will simply add up all values from 1 to 10

sum:=0:
for i from 1 to 10 do
  sum:=sum+i:
end do:

Here i is the index variable which counts from 1 to 10. If no value is given as from, the default value of 1 is assumed. The default value by which the counter is incremented is 1. However, this value can also be explicitly set.

try/catch/end try. For all readers familiar with programming, the concept of try/catch/end try constructs is well known. Maple will try to execute code enclosed in a try block. If it fails to execute the code, e.g., because one of the functions returns an error, the program is not aborted, but continued at the catch section. Using such constructs, it is often possible to test a potential solution, but avoiding program crashes in case of failure.

An example of a try/catch block can be found in listing 27.3 where we find the following expression

try
 tosolve:=Fi=ret[i-1,2]+(df(ret[i-1,1],ret[i-1,2])+df(ret[i,1],Fi))*delta/2:
 ret[i,2]:=rhs(solve({tosolve},{Fi})[1]):
catch:
 print(sprintf(“Cannot solve %a”,tosolve));
 return;
end try:

Here, Maple will try to solve the equation assigned to tosolve. Due to the nature of the equation, it may be impossible for Maple to come up with a solution (for more details on this code segment, check listing 27.3). Usually, Maple would simply abort in such a case. However, as we have enclosed the code in a try/catch block, the error will be caught and the code execution continues by outputting information to the user.

Working With Units.Maple allows calculating with units. As an example, we can assign a unit to a value by typing

p:=1*Unit(mbar/mm);

Maple allows very convenient operations on units, including unit conversions. Please note that we need to load specific packages in order to allow unit calculations in Maple. We will see plenty of examples for calculations with units throughout this book.

Summary. This section summarized the most important concepts and functions we require when working with Maple. As stated, it is not intended to be a complete reference to Maple, but should merely serve as a quick reference on the most important concepts. Feel free to explore the program for yourself and become familiar with it by going through the very good tutorials.

2.3 The File Core.txt

Throughout this book we will use the file Core.txt in which we store several useful Maple procedures. The complete listing of the file is shown in listing 2.1. In the following, we will briefly explain the individual procedures without going too much into detail as to their implementation.

Listing 2.1

[Maple] The file Core.txt which contains a set of useful functions we will require regularly. The file can be found in the supplementary material.

1 with(plots): with(Units): with(Statistics): with(PolynomialTools): with(ColorTools): with
     (plottools): with(FileTools): with(DEtools): with(Optimization): with(LinearAlgebra):
    with (Student[NumericalAnalysis]):with(ArrayTools):with(ScientificConstants):
2
3 #this function generates a 2D plot
4 quickplot := proc(content::list,rangeX::`..`,rangeY::`..`:=NULL,labelx:=“x”,labely:=“y”,legends::
    list:=NULL,IsPointPlot::boolean:=false)
5    return plot(content,rangeX,`if`(rangeY=NULL,NULL,rangeY),labels=[labelx,labely],legend=legends,`
   if`(IsPointPlot,style=point,style=patch),_rest):
6 end proc:
7
8 #this function generates a 3D plot
9 quickplot3d := proc(content::list,rangeX,rangeY,rangeZ::`..`:=NULL,labelx:=“x”,labely:=“y”,labelz :=“z”)
10 return plot3d(content,rangeX,rangeY,`if`(rangeZ=NULL,NULL,view=rangeZ),labels=[labelx,labely, labelz],plotlist=true,_rest):
11 end proc:
12
13 #this function extracts two columns from an array for plotting
14 extracttoplot:=proc(data,extract) local _points,_show,i;
15   #check if extract contains exactly two entries (x and y)
16    if numelems(extract)<>2 then
17      printf(“ERROR: extract must contain two entries [x,y]”);
18      return;
19    end if;
20   _points:=Size(data,1):
21   _show:=Matrix(_points,2);
22   for i to _points do
23     _show[i,1]:=data[i,extract[1]]:
24     _show[i,2]:=data[i,extract[2]]:
25   end do:
26   return _show;
27 end proc:
28
29 #this function builds an estimation table
30 buildestimation:=proc(_from,_to,functionnames,functions,showestimated::boolean:=true,points::integer:=100,inlineresults::boolean:=true) local i,j,_estimated,_evaluatedat,_output,_toplot, _plots;
31    #check if the number of function names matches the given functions
32   if numelems(functionnames)<>(numelems(functions)+1) then
33     printf(“ERROR: Number of function names (%d) must be number of functions (%d) + 1 (we don’t need a function for the function values at which we want to evaluate)”,numelems(functionnames),numelems(functions));
34   return;
35    end if;
36    #this is the return array
37    _evaluatedat:=Array(1..points);
38    #we need an array with the entries - the names of the functions and their values
39    _estimated:=Array(1..2,1..1);
40    #these are the function names
41    _estimated[1,1]:=Array(functionnames);
42    #and this is where we will plot the values
43    _estimated[2,1]:=Matrix(points,numelems(functionnames));
44    #now we assemble the estimations
45    for i to points do
46   #the first value needs to be the lower range value exactly otherwise we will get an error
47   if i=1 then
48     _evaluatedat[i]=evalf(_from);
49   else
50     if i=points then
51       #the last value needs to be the upper range exactly otherwise we will get an error
52       _evaluatedat[i]:=evalf(_to);
53     else
54       #we interpolate the values in between
55       _evaluatedat[i]:=evalf(((_to-_from)/(points-1))*(i-1)+_from);
56     end if:
57    end if:
58    _estimated[2,1][i,1]:=_evaluatedat[i];
59    #and calculate the other functions
60    for j to numelems(functions) do
61       _estimated[2,1][i,j+1]:=evalf(functions[j](_evaluatedat[i]));
62   end do:
63    end do:
64    #now we reassemble this into an array
65    _output:=Array(1..2):
66    _output[2]:=_evaluatedat:
67    _output[1]:=_estimated:
68    #do we want to display
69    if showestimated then
70       _toplot:=Array(1..numelems(functionnames));
71       for j to numelems(functionnames) do
72         _toplot[j]:=Array(1..points):
73       end do:
74     for i to points do
75      for j to numelems(functionnames) do
76      _toplot[j][i]:=_estimated[2,1][i,j];
77     end do:
78    end do:
79    #prepare the lists for plotting
80    _plots:=Array(1..numelems(functionnames));
81    for j to numelems(functionnames) do
82    _plots[j]:=plot(zip((X,Y)->[X,Y],_evaluatedat,_toplot[j]),axes=boxed,style=point,title =functionnames[j]);
83    end do:
84    print(display(_plots)):
85   end if:
86   #which return do we want
87   if inlineresults then
88       return approxsoln=_output[1],output=_output[2],range=_from..._evaluatedat[points];
89   else
90     return _output:
91   end if:
92  end proc:
93
94  #this function plots the results of a numerically solved ODE
95  plotode:=proc(results, toplot) local i,j,_plots,_toplot,_points,_evaluatedat;
96   _points:=Size(results[2,1],1):
97   _toplot:=Array(1..numelems(toplot));
98    _evaluatedat:=Array(1.._points);
99    for j to numelems(toplot) do
100     _toplot[j]:=Array(1.._points):
101    end do:
102    for i to _points do
103      _evaluatedat[i]:=results[2,1][i,1];
104      for j to numelems(toplot) do
105        _toplot[j][i]:=results[2,1][i,toplot[j]];
106     end do:
107    end do:
108    #prepare the lists for plotting
109      _plots:=Array(1..numelems(toplot));
110      for j to numelems(toplot) do
111         _plots[j]:=plot(zip((X,Y)->[X,Y],_evaluatedat,_toplot[j]),axes=boxed,style=point,title= results[1,1][toplot[j]]);
112    end do:
113    print(display(_plots)):
114  end proc:

2.3.1 Includes

As can be seen in line 1, the file opens with a number of Maple packages being loaded. These are all the packages which are required throughout this book. Loading them from the file Core.txt makes the actual worksheets significantly shorter, as we do not have to choose the packages to include ourselves. Obviously, not all worksheets will require all included packages. However, the only downside of loading too many packages is a slight loss in time which is hardly noticeable.

2.3.2 quickplot

One function which we will require over and over again is the function quickplot. This function is defined in line 4. The function expects a list of functions to plot as the argument content as well as the ranges in x- (rangeX) and y-direction (rangeY). The latter argument is NULL by default and can be omitted if the range in y-direction should not be restricted. Next, the function expects the labels in x- and y-direction labelx and labely which are “x” and “y” per default. The next argument is a list of strings which would be used as a legend in case more than one function must be plotted. In this case, the user could provide, e.g.,[“function 1”,”function 2”]. The last argument is a boolean value which decides if the function is to be plotted as a solid line or as points, i.e., as discreet values.

In listing 27.5 we use this function to plot the function e-t2si5_e in the interval between 0 and 4 using the line

quickplot([t->exp(-t/2)],0..4,"x","f(x)");

As you can see, the definition of quickplot passes the Maple primitive _rest when calling plot. This variable stores all potential addition arguments which may be passed onto quickplot. Examples include, e.g., the plot symbols if a point plot is to be created. Please refer to listing 27.9, line 9 for an example of such an additional argument.

2.3.3 quickplot3d

Similar to the function quickplot, we define a second function which is termed quickplot3D (see line 9). As the name suggests, this function will create a three-dimensional plot. It takes a list of functions to plot as the argument content as well as the ranges in x- and y-direction as the arguments rangeX and rangeY. Both of these arguments are mandatory. The range along the z-axis is optional and will be chosen by Maple if not explicitly set. The arguments labelx, labely, and labelz can also be passed as arguments to the function. Omitting them will result in the axes simply being labeled “x”, “y”, and “z”.

2.3.4 extracttoplot

The function extracttoplot defined in line 14 is used primarily for displaying the result of numerical calculations. It expects a matrix as parameter data and an array of two indices as parameter extract. It will then extract the two columns denoted by the indices into a two-column matrix. The first index will indicate the matrix column for the x-values, whereas the second index indicates the column for the y-values. The created array is then returned. Usually, this function is required in order to plot array data. For an example, please refer to listing 27.7.

2.3.5 buildestimation

The function buildestimation is a rather lengthy piece of code which is used for creating estimations required for numerical solutions of differential equations. We use this function extensively in section 24, as an example see listing 24.1. The function expects the range of the independent variable given as lower boundary _from and upper boundary _to. It also expects the name of the functions that should be solved for as parameter functionnames. These must be given in form of an array. The actual functions are given as arrays in the parameter functions. The parameter showestimated is a boolean value. If this value is true, the function will output plots of the estimation functions. Otherwise the estimations are calculated, but not displayed to the user. This mainly serves as a means of checking whether or not the given functions look the way they are supposed to look. The next parameter points indicates the total number of points to create. The default value is 100. The last parameter inlineresults is mainly for testing purposes and can be omitted.

The function first checks if the given number of function names matches the number of functions given (see line 31). It then creates the array _estimated which holds the data. This array has two entries: The first entry is again an array with the name of the functions (see line 41). The second entry is a matrix which holds the values of the respective estimated function at a given value of the independent variable (see line 43). The function then fills this matrix with the values of the estimated functions with values (see line 45). It first calculates the current value of the independent variable and stores it in _evaluatedat[i]. It then calculates the value of the respective estimation function for this value of the independent variable and stores it at the correct location of the array _estimated (see line 61).

Finally, the procedure creates a single array termed _output into which the arrays _estimated and _evaluatedat are grouped (see line 65). If the user requested to see the functions, the procedure will create plots and output them for the user to see line 69.

2.3.6 plotode

The last function that is defined in Core.txt is the function plotode (see line 95). This function essentially extracts the numerical data obtained from the numerical solver in Maple and copies the data into array format. This is required in order to be able to plot the result using the function quickplot. We will see an example of this function in listing 24.1.

2.4 The File Corefunctions.txt

The second file which we will require is the file CoreFunctions.txt which is shown in listing 2.2. This file is a collection of useful functions which we use frequently when studying numerical solutions to equations. The file defines the following functions

Listing 2.2

[Maple] The file CoreFunctions.txt which contains a set of useful numerical functions. The file can be found in the supplementary material.

1   #fourth-order Runge-Kutta method for solving ODE systems
2   RungeKutta4ODEs:=proc(df,dependent,_from,_to,initialvalues,points) local ret,numdep,delta,i,j,c:
3  numdep:=numelems(dependent):
4  ret:=Matrix(points,1+numdep):
5  c:=Matrix(numdep,4):
6  delta:=(_to-_from)/(points-1):
7  ret[1,1]:=_from:
8  for j to numdep do:
9   ret[1,1+j]:=initialvalues[j]:
10  end do:
11  for i from 2 to points do:
12   ret[i,1]:=evalf(ret[i-1,1]+delta):
13   for j to numdep do:
14     c[j,1]:=evalf(delta*df[j](ret[i-1,1],seq(ret[i-1,1+k],k=1..numdep))):
15   end do:
16   for j to numdep do:
17     c[j,2]:=evalf(delta*df[j](ret[i-1,1]+delta/2,seq(ret[i-1,1+k]+c[k,1]/2,k=1..numdep))):
18   end do:
19   for j to numdep do:
20     c[j,3]:=evalf(delta*df[j](ret[i-1,1]+delta/2,seq(ret[i-1,1+k]+c[k,2]/2,k=1..numdep))):
21   end do:
22   for j to numdep do:
23     c[j,4]:=evalf(delta*df[j](ret[i-1,1]+delta,seq(ret[i-1,1+k]+c[k,3],k=1..numdep))):
24   end do:
25   for j to numdep do:
26     ret[i,1+j]:=evalf(ret[i-1,1+j]+1/6*(c[j,1]+2*c[j,2]+2*c[j,3]+c[j,4])):
27   end do:
28   end do:
29   return ret;
30   end proc:
31
32  #successive over-relaxation
33  DoSOR:=proc(m::Matrix,b::Vector,x0::Vector,omega,kmax:=10) local i,j,k,xcurrent,summed,n:
34    n:=Size(m,1):
35    xcurrent:=Vector(1..n):
36    for i to n do
37     xcurrent[i]:=x0[i]:
38   end do:
39   for k to kmax do
40     for i to n do
41      summed:=0:
42     for j to n do
43      if i<>j then
44        summed:=summed+m[i,j]*xcurrent[j]
45     end if:
46   end do:
47   xcurrent[i]:=(1-omega)*xcurrent[i]+omega*evalf((b[i]-summed)/(m[i,i])):
48  end do:
49  end do:
50  return xcurrent:
51 end proc:
52
53 #sequential sense LU-decomposition with pivoting
54 DoLUDecompositionSequentialDensePivot:=proc(m::Matrix) local n,i,j,k,l,piv,store:
55    n:=Size(m,1):
56    piv:=Vector(1..n):
57    for i from 1 to n do
58     piv[i]:=i:
59   end do:
60   for i from 1 to n do
61    for l from i to n do
62     if abs(m[piv[l],i])>abs(m[piv[i],i]) then
63    store:=piv[l]:
64    piv[l]:=piv[i]:
65    piv[i]:=store:
66    end if:
67    end do:
68    for j from i+1 to n do
69      m[piv[j],i]:=m[piv[j],i]/m[piv[i],i]:
70    end do:
71    for k from i+1 to n do
72     for j from i+1 to n do
73     m[piv[k],j]:=m[piv[k],j]-m[piv[k],i]*m[piv[i],j]:
74     end do:
75    end do:
76   end do:
77   return piv,m;
78   end proc:
79
80 #sequential dense LU-substitution: solving linear systems of equations
81 DoLUSubstitutionDensePivot:=proc(piv::Vector,m::Matrix,b::Vector) local n,i,k,x,y,summed:
82   n:=Size(b,1):
83   y:=Vector(1..n):
84   x:=Vector(1..n):
85   #forward-substitution - find y[i]
86   for i from 1 to n do
87    summed:=0:
88    for k from 1 to i-1 do
89      summed:=summed+m[piv[i],k]*y[k]:
90    end do:
91    y[i]:=evalf(b[piv[i]]-summed):
92   end do:
93   #back-substitution - find x[i]
94   for i from n by -1 to 1 do
95    summed:=0:
96    for k from i+1 to n do
97      summed:=summed+m[piv[i],k]*x[k]:
98    end do:
99    x[i]:=evalf(1/(m[piv[i],i])*(y[i]-summed)):
100  end do:
101  return x:
102 end proc:
103
104  #Thomas algorithm for solving tridiagonal linear systems of equations
105  DoThomas:=proc(m::Matrix,b) local n,x,i:
106   n:=Size(m,1):
107   x:=Vector(1..n):
108   m[1,2]:=m[1,2]/m[1,1]:
109   b[1]:=b[1]/m[1,1]:
110   for i from 2 to n-1 do
111    m[i,i+1]:=m[i,i+1]/(m[i,i]-m[i,i-1]*m[i-1,i]):
112    b[i]:=(b[i]-m[i,i-1]*b[i-1])/(m[i,i]-m[i,i-1]*m[i-1,i]):
113   end do:
114   b[n]:=(b[n]-m[n,n-1]*b[n-1])/(m[n,n]-m[n,n-1]*m[n-1,i]):
115   x[n]:=evalf(b[n]):
116   for i from n-1 to 1 by -1 do
117    x[i]:=evalf(b[i]-m[i,i+1]*x[i+1]):
118   end do:
119   return x;
120  end proc:

 RungeKutta4ODEs: the fourth-order Runge-Kutta method which is discussed in section 27.2.8; the implementation of the function is explained in listing 27.8

 DoSOR: the successive over-relaxation which is discussed in section 25.3.4; the implementation of the function is explained in listing 25.3

 DoLUDecompositionSequentialDensePivot: the sequential dense LU-decomposition with pivoting discussed in section 25.3.6; the implementation of the function is explained in listing 25.10

 DoLUSubstitutionDensePivot: the LU-substitution discussed in section 25.3.7; the implementation of the function is explained in listing 25.11

 DoThomas: the Thomas-Algorithm discussed in section 25.3.8; the implementation of the function is explained in listing 25.12

As stated, we will not detail the specifics and the implementation of these functions at this point. This will be done in later chapters.

2.5 The Neptunlib

Throughout this book, we will implement a number of functions in C. For this we will build a custom dynamically linked library (DLL) termed NeptunLib.dll. In order to build this DLL, the files NeptunLib.c (which contains the code) and NeptunLib.h (which is the header file) are required. Both can be found in the supplementary material to this book.

Compiler Options. The files were created, compiled, and tested on a Windows 7 (64-bit) using the integrated development environment (IDE) Eclipse with gcc as the compiler of choice. The compiler options used were

-std=c99 -O3 -c -fmessage-length=0

These options set the C dialect to c99 and the optimization level to highest (-O3).

Linker Options. The linker options used were

fvisibility=hidden -Wl,--add-stdcall-alias -shared

which sets the mangling such that Maple can read the exported functions of the DLL without problems. The type of library created was a shared library (shared),i.e., a DLL.

You can find a precompiled version of the NeptunLib.dll in the supplementary material if you do not wish to compile and link it yourself.

2.6 Summary

In this section we have briefly introduced the Algebra package Maple and highlighted some of its most important features. This should serve as a sufficient introduction for everyone interested in following and potentially extending the worksheets and the listings which are described in this book. As stated, there are many great online resources for readers more interested in the specifics and the capabilities of Maple. We have also introduced the two files Core.txt and CoreFunctions.txt, which contain a number of useful functions which we will use throughout this book.

References

[1] Geddes K., Labahn G., Monagan M. Maple 12 advanced programming guide. 2008 (cit. on p. 9).

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset