calculation
Insert and edit math into text (use window)
Introductory Scientific Computing
A Maple Tutorial*
copyright, January 2001 , RH Landau
Physics Department, Oregon State University
Maple Building Blocks
From Unix Prompt (= %)
% xmaple & Start X Windows version of Maple in background (background = &)
% maple & Start text version of Maple in background (background = &)
Select Modes
T "text" input Useful for writing notes or writing up homework on worksheet
[> Maple input To enter mathematics
calculation
Insert and edit math into text (use window)
Execute Commands
Place mouse anyplace after > and press [Enter] to execute.
> restart; De-assign all variables (like a new session); good to do for new group
> 3+7; Enter and Return may/may not be same key
>
Basic Operations & Representations
> (4/2)*3-7^8+7**8; After many operations (2 forms of power 8) this should be 6
> 2485/3479; Maple reduces to 5/7, no loss of precisio n
> 2485./3479; Decimal point => float; default to 10 digits, loss of precision
> Digits := 40; Change number of digits variable to 40 (NB := not =)
> 2485./3479; Try again with 40 digits now
> evalf( 2485/3479, 20 ); Evaluate rational number as float, (optional: 20 digits)
> 1.0*10^(-6); Floating point number in scientific notation
> Float( 1, -6 ); Alternative form of above
> 1E-6; Alternative form of above, not `` e''
> 1.*E-6; Variable E to -6 power
> log(1.*E-6); test
> exp(1.); exp(1); Two versions of e (base of natural logs)
> log( exp(1) ); test
> 1.E-6; Syntax error without the *, or with the "."
> sqrt(10!); Maple even simplifies for you
> Pi; Name for 3.1415...
> ln(x); log10(x); abs(x); cos(x); Functions needing argument x
> arccos(x); cosh(x); P(l,x); BesselK(k,x); "
> sin(Pi); sin(pi); The first one is what you want unless you define pi
> (2+3*I)/(5-4*I); Complex numbers
> convert( %, polar ); Polar form of complex number
> ?help; Use ? in front of anything for help on that command
>
Sums & Products
> sum( (-1)^i * x^(2*i)/(2*i)!, i=0..2 ); Evaluate first 3 terms in cos series
> sum( (-1)^i * x^(2*i)/(2*i)!, i=0..infinity ); Evaluate all terms
> Sum( (-1)^i * x^(2*i)/(2*i)!, i=0..infinity ); Symbolic sum
> value(%); Evaluate previous expression (% replaces old ")
> Product( (i^2+3*i-11)/(i+3), i=0..10 ); Form product of symbols
> value(%); Evaluate previous expression
>
Symbolic Manipulations
> restart; Clean the slate
> Sum( i, i=1..n ); This should be n(n+1)/2 when evaluated
> sum( i,i = 1 .. n ); The active form. lowercase not upper
> value(%); Not form of expected result
> expand(%); factor(%); expand(%); so be manipulative!
> expand( (x+y)^15 );
> simplify( cos(x)^2 + sin(x)^2 ); Simplify using math identities
> convert( sin(x), exp ); Convert to exponential notation (many variants)
>
Statements, Expressions, Functions & Procedures
> restart; De-assign all variables (like new workshhet)
> 1+4; exp(x)-1; Expressions (evaluated or stored symbol s)
> x :=3; y := 3* z^2 + 1; Assign value or expressions to x & y; has side effects
> eval( y, z=2 ); Evaluate y for z=2
> f :=(x) -> x* sin(x); Defines function f(x) as x sin(x)
> R :=(x,y,z) -> sqrt( x*x+y*y+z*z ); Define function of 3 variables
> g := unapply( 3*z^2+1, z ); g(2); Convert expression into function
> x := 'x'; De-assign x
> z := x * sin(x); A procedure setting variable z equal to expression
> zfun := unapply( z,x ); Define function zfun(x) from procedure z
Solving Equations
Single Equation (subscripted variable)
> solve( a* x^2 +b*x +c , x ); Search for x satisfing a*x^2+b*x+c = 0
> solve( a*x^2 +b*x +c = 0 , x ); Alternate form of above
> y := a* x^2 +b*x +c = 0; Alternate form using variable
> solve( y, x ); "
> eval( y, x = a/2 ); Verify answer
> roots( x^2-1 ); Gives real, interger roots & multiplicity of polynomials
> roots( x^2 +1, I ); Gives complex, interger roots & multiplicity
> ans := solve( a* x^2 +b*x +c =0, x); Form sequence (subscripted variables)
> ans[2];
> fsolve( x*sin(x) -1 , x ); Solve for one floating point (numeric) values of x
> fsolve( x*sin(x) -1 , x=0..6 ); Solve for 1 numeric values of 0 < x <6
>
Simultaneous Equations
> eqn1 := a +3*b +4*c = 41; 2 equations, 3 unknowns
> eqn2 := 5*a +6*b +7*c = 20; "
> solve( {eqn1, eqn2}, {a, b} ); Solve for a,b wrt c
> x := Linsolve(A,b) mod 5; Solve matrix equation [A] x = b (see too linear alg)
>
Plotting Along
> restart; Clean the slate
> with( plots ); Load extended plotting package
> with( plottools ); Load extended plotting tools
> plot( tan(x), x=-7...3 ); Basic plot
> plot(tan(x), x =-7..7, y =-5..5, labels =[`x`,`tan(x)`], title =`y=tan(x)`); With embellishments, strings in back quotes
> plot( { cos(x), sin(x), x^2 }, x = -Pi..Pi ); Multiple functions on one plot
> plot( max(0, cos(x)), x = -2*Pi..2*Pi ); Limit range to positive [also min(.), abs(.) ]
> plot([sin(x), cos(x), x =0..2*Pi], scaling = constrained); Phase-space, same x,y scale
> implicitplot({x^2+y^2 = 1, y=exp(x)},x =-2..2, y =-2..2); Implicit 2 variable eqtn
> plot3d( x*(x^2-3*y^2), x=-1..1, y=-1..1,title='saddle' ); 3D plot of saddle function
> animate( cos(t*x) *sin(t*x), x = 0..Pi, t=1..20 ); 2D animation, click graph for controls
> animate3d( cos(t*x) * sin(t*y), x =0..Pi,y=0..Pi, t=1..2 ); 3D animate; click for controls
>
Plotting of Data
> with(statplots);with(stats); Need statistics packages (not math!)
> Xdata := [1,2,3,6,10]; A list of all x values
> Ydata := [1,4,9,36,100]; A list of corresponding y values
> scatterplot( Xdata, Ydata, color=RED); Can have several plots
>
Calculus
Differentiation
> restart; Clean the slate
> f := (x) -> ln(x)/(1-x); Define function to differentiate
> diff( f(x), x ); Derivative of expression wrt x
> diff( f(x), x$2 ); Second derivative (nb $ not ^ )
> diff( (sin(x))^y, x ); Partial derivative of f(x,y) wrt x
> diff( (sin(x))^y, x,y ); Second partial derivative of f(x,y) wrt x and y
> Diff( f(x), x ); Inert form of above
> f_prime := value(%); Evaluate derivative and call it f_prime
> D(f)(x); Derivative of function f wrt x . ( D = operator, Diff = command). (See below)
> Diff( f(x,y,z), x,x,y,z,z,z ); Differentiate multivariate functions
> g := (x) -> D(f)(x); Define new function as derivative of another function
> h := (x,y,z) -> sin(x*y*z); Define multivariate function
> g := D [1] (h); Differentiate h wrt its first argument
> Limit( (1+x/n)^n, n =infinity ); Limit (then value(%) or use limit )
>
D - Differential operator
> D (cos); Derivative of cos (transforms one function into another)
> D (cos) (x); Explicit dependence on x shown
> D (f) (x); Derivative
> g := (x) -> D(f)(x); Define new function, derivative of other function
> h := (x,y,z) -> sin(x*y*z); Define multivariate function
> g := D[1](h); Differentiate h wrt its first argument
>
Integration
> restart; Clean the slate
> Int( f_prime, x); Antiderivative (integral)
> value (%); "
> simplify (%); "
> Int ( x/sqrt(1+x), x ); value(%); Indefinite integrate wrt x (inert); evaluat e
> int ( x/sqrt(1+x), x ); Indefinite integral wrt x (implied constant of integration)
> diff ( %, x ); Check answer by differntiation
> Int ( x/sqrt(1+x), x = 0..1 ); Evaluate definite integration (inert)
> int ( x/sqrt(1+x), x = 0..1 ); evalf(%); Evaluate definite integration (active); floats
> g := unapply ( int( f(x), x ),x ); Integrate to produce a function
> g := (y) -> int( f(x,y), x = a..b ); Defines function of y after integrate out x
>
Series Expansions
> ans := series ( exp(x)/x, x=0 ); Power series expansion about x=0
> truncate := convert ( ans, polynom ); Truncate series to polynomial for plottin g
> plot ( { exp(x)/x, truncate }, x = 0..4 ); Plot exact vs truncated series
> Order := 12; Trunacate as 12th order polynomial
> truncate := convert(ans, polynom); "
>
Differential Equations
First Order ODE
> restart; Clean the slate
> with( DEtools ); Load differential equation tools
>
diff_eq := D( D(x) ) (t) = -w*w *x;
Simple harmonic motion,
> diff_eq := ( D@@2 ) (x) (t) = -w*w *x; Alternate form of above
>
Second Order ODE with Plot
> restart;
> with( DEtools );
> DEplot( [diff ( theta(t), t ) - omega = 0, diff( omega(t), t) + sin(theta) = 0], [t, theta, omega], 0..12, { [0,0,1.0] }, scene = [t,omega] );
>
Second Order PDE
> restart;
> with( DEtools );
>
diff_eq1 := D ( D(y) ) (x) + 5*D(y)(x) + 6*y(x) = 0;
;
operator acts on expressions, not functions
> init_con := y(0)=0, D(y)(0)=1; Initial conditions
> dsolve( {diff_eq1, init_con} , {y(x)} ); Solve the equation
>
System of ODE's
>
sys := (D@@2)(y)(x) = z(x), (D@@2)(z)(x) = y(x);
> dsolve( {sys}, {y(x), z(x)}); Maple will generate initial condition parameters
>
Linear Algebra
Defining Matrices &Vectors
> with(LinearAlgebra); Need for more than basics; good way to see what's available
> M := matrix(3,3, [4,-2,1,3,6,-4,2,1,8]); Define 3 x 3 matrix with 9 elements
> N := matrix(2,2, [2,x,1, y^3]); Define 2 x 2 matrix with symbol s
> A := matrix(3,3, (i,j)-> i*j); Matrix with function as elements
> B := array(1..3, 1..3, [[-2,2,-3],[2,1,-6],[-1,-2,0]]); Another way to do matrix M above
> ID := IdentityMatrix(3);
> v := vector([1,2,3]); Enter column vector with list; displayed as row with commas
> ucol := matrix(3, 1,[1,2,3]); An explict column vector
> urow := matrix(1,3,[1,2,3]); An explicit row vector, NB no commas!
Matrix Operations
> evalm(A &* B); Matrix multiply of too square matrices
> evalm(A &* ID); Multiply by identity matrix
> evalm(A^2); Square of matrrix
> evalm(1/B); inverse of nonsingular matrix
> det(A); Find determinant (0 => singular)
> inverse(N); Find inverse another way, even symbolic.
> P := multiply(A,B); Matrix multiplication another way
> evalm(A &* u); Multiply square matrix by column vector
> multiply(A, u); Or in other words,
> multiply (u, A);
Eigenvalues, Eigenvectors, Linear Equations
> eigenvals(B); Eigenvalues of matrix B, note double root
> eigenvects(B); Eigenvectors of 3x3 matrix (eigenvalue,multiplicity,vectors; ...)
> A := matrix(3,3, [4,-2,1,3,6,-4,2,1,8]); ) Define matrix
> b := vector([12, -25,32]; Define vector
> x := linsolve (A,b); Solve A x = b for unknown vector x
Generating Fortran & C Code (Maple V)
> with (linalg);
> A := matrix ( 3, 3 ); Define 3 x 3 matrix
> print (A); View matrix
> A_inv := inverse( A ); Define A_inv to contain symbolic inverse
> codegen[ fortran ] (A_inv); Produce Fortran Code
> codegen[C] (A_inv); generate C code
>
Statistics & Fitting
Initialize (load)
> restart;
> with ( stats );
>
Least Square Fits
> with ( fit ); Load fitting package
> X := [2,4,6,8]; Define independent variable
> Y := [0,20,28,44]; Define dependent variable
> leastsquare [ [x,y], y = a*x^2+b*x+c, { a, b, c } ] ( [X,Y] ); Fit quadratic
>
Mean & Variance
> with ( describe ); Load package for means et al
> mean_X := mean ( X );
> median ( X );
> variance ( X );
> dev := standarddeviation ( X );
>
Statisitcal Plots (Maple V)
> with(statplots); Load statistical plotting package
> X := [4.5,4.0,5.4,1.6,5.7,3.5,7.8,8.1,6.0,13.4,2.8,5.1,3.7,8.0,3.6, 7.2,3.8,2.1,2.6,4.6]; X data.
> Y := [7.4,4.4,2.8,5.4,9.9,-1.4,1.0,1.1,4.8,0.45,-0.78,9.9,4.8, -3.1,4.4,5.3,1.91,-0.79,1.4,1.6]; Y data.
> scatterplot( X, Y, color=black); Make scatterplot .
>
Maple Procedures (see too Generating Fortran & C Code)
> restart;
> sub := proc (x, y)
x^y;
end; Procedure sub to raise x to the y power
> sub (2,5); This should be 2^5=32
> sub (i,j); This should be algebraic, i^j
>
With loops, ifs, and local variables
>
Chebyshev := proc( n )
local p, k;
p[0] := 1;
p[1] := x;
if n <= 1 then
RETURN ( eval(p) )
fi;
for k from 2 to n do
p[k] := expand( 2*x*p[k-1] - p[k-2] )
od;
RETURN ( eval(p) )
end;
> a := Chebyshev(5); a[3];
>
* Development supported in part by the National Science Foundation Curriculum Development Grant NSF?? and the Education, Outreach and Training Thrust Area of the National Partnership for Advanced Computational Infractrure (NPACI).