- HOME
- HELP

Tweet

Programs can be created and run interactively. Programming a function is demonstrated in the following example of a function ttwo(x), which multiplies its argument by 2. After the definition it can be used like any other Jasymca function.

>> function y=ttwo(x) y=2*x; end >> ttwo(3.123) ans = 6.246

Following the keyword function is the prototype with a return variable y. This replaces the construct return y of other programming languages. If functions are to be reused later, they should be written to a textfile and saved somewhere in Jasymcas searchpath. The filename must be the function name extended by ?.m?, in the present example ttwo.m. In subsequent sessions the function ttwo can be used without separately loading the file. Several installed functions of Jasymca are provided using this mechanism.

?if x A end? Depending on the condition x one or several statements A are executed. The condition x must be an arbitrary expression, which evaluates to either 0 or 1. The false-case (i.e. x=0) can lead to another branch B: ?if x A else B end? As an example the Heavyside function:

>> function y=H(x) > if (x>=0) > y=1; > else > y=0; > end > end >> H(-2) y = 0 >> H(0) y = 1

Loops with condition x and statement(s) A: while x A end The while-loop is repeated until x becomes false (0).

>> x=1;y=1; >> while(x<10) y=x+y; x++; end >> y y = 46

Loops with counter z and statement(s) A: ?for z = vector A end? In the for-loop the counter is formally initialized by a vector. In each execution of the loop the counter takes on the value of the next element of vector.

>> x=1;y=1; >> for(x=1:0.1:100) y=x^2+y; end >> y y = 3.3383E6

return, continue, break A function may be prt of the loop, and begins another cycle. break permanently leaves the loop.

>> x=1; >> while( 1 ) > if(x>1000) > break; > end > x++; > end >> x x = 1001

Numbers are entered in the usual Computer format: with optional decimal point and decimal exponent following the letter e (or E). The numbers 5364 and $-1.723478265342\; 10^\{12\}$ should be entered like:

>> 5364 ans = 5364 >> -1.723478265342e12 ans = -1.7235E12

Most of the time these will be stored as floating point data (double, IEEE standard 754). They are rounded to 5 significant digits for display, but for calculations the full precision of this format is always preserved (15-16 decimal digits). By switching the format (format long) all significant places are displayed.

>> format long >> -1.723478265342e12 ans = -1.723478265342E12

As an extension Jasymca offers the command format Base Number, which is used to display numbers in a system with arbitrary Base with any Number of significant digits. To display numbers with 15 digits in the binary system we type:

>> format 2 15 >> -1.723478265342e12 ans = -1.1001000101001E40

Using format short returns the display mode to default (short decimal). It should be emphasized, that none of the format commands influences the internal representation and accuracy of floating point numbers. Numbers, which are entered without decimal point and exponent, and which are larger than $10^\{15\}$ are stored as exact rational datatype. These numbers are internally represented as quotient of two variable length integers (java datatype BigInteger), which allows you to perform calculations without any rounding errors. In the first case of the following example a floating point number is generated, in the second case an exact rational:

>> 10000000000000001. ans = 1.0E16 >> 10000000000000001 ans = 10000000000000001

Each floating point number Z can be converted to an exact number using the command rat(Z). The conversion is accomplished by continued fraction expansion with an accuracy determined by the variable ratepsilon (default: $10^\{-8\}$).

>> rat(0.33333333333333333) ans = 1/3

Operations between exact and floating point numbers always lead to the promotion of floating point numbers. Calculations can be performed without rounding errors by rationalizing just the first number.

>> 1/21/525/21/5*7*175*63*15-1 ans = -4.4409E-16 >> rat(1)/21/525/21/5*7*175*63*15-1 ans = 0

Conversely, the command float(Z) converts numbers into floating point format. Both commands also work for composite datatypes, like polynomials and matrices, whose coefficients are transformed in one step. Irrational function values of exact numbers and constants like pi remain unevaluated until the float-command is issued.

>> sqrt(2) ans = 1.4142 >> sqrt(rat(2)) ans = sqrt(2) >> float(ans) ans = 1.4142

The exact datatype is useful especially for unstable problems, like solving systems of linear equations with ill-conditioned matrix. The Hilbert-matrix is an extreme example:

>> det( hilb(20)*invhilb(20) ) ans = 1 % correct >> det( float(hilb(20))*float(invhilb(20)) ) ans = 1.6713E151 % slightly wrong

Imaginary numbers are marked with an immediately following i or j. This will work even if the predefined variables i and j have been overwritten.

>> 2+3i
ans = 2+3i

How to print out the output values if you are working with the script, rather than with the prompt. Just use ?printf? method:

x=log(sqrt(854)); % natural logarithm
printf('Answer=%f\n', x);

In contrast to the examples Octave and Matlab, Jasymca integrates numeric and symbolic datatypes at the core of the program; symbolic math is not treated as an add-on. This means that with few exceptions most operations accept any mixture of numeric and symbolic arguments using the same commands and command syntax. Symbolic variables should not be confused with variables as discussed until now. These latter variables serve as address for an object in memory (the ?environment?), while symbolic variables are algebraic objects on their own. That means if x is a conventional variable, entering x in the text input field makes Jasymca search in the environment for the corresponding object, which then replaces x. If however x is a symbolic variable, the same action will lead to the creation of a first-degree polynomial with variable x and coefficients 1 and 0. In Octave-mode, each symbolic variable x must be declared as symbolic by entering syms x before using it. The command clear x deletes the symbolic (actually any) variable x.

>> x=3; % nonsymbolic variable >> x^2+3-2*sin(x) % placeholder for '3' ans = 11.718 >> syms x % symbolic variable >> x^2+3-2*sin(x) % create function ans = -2*sin(x)+(x^2+3)

Variables are declared by supplying a name and value in the format name=value. The name can be any character sequence. With the exception of the first character it may also contain numbers. The value is any number or expression.

>> x=24+3i x = 24+3i

Some variables are predefined (like pi). The last previous result of a
calculation is stored in the variable ans. All variables are displayed by
the command who. Single variables can be deleted by entering clear variable.
It is possible to define variables whose value is a function. In this case the
function's name must be preceded by the character tosuppressevaluation.
Thesevariablescanbeusedlikethefunctiontheystandfor.Forexample,
whodislikesthebuiltinfunctionrealpart(x)?snamecanshortenittotheMatlab?version:
` matlab> >>real=realpart realpart>>real(24+3i)ans=24`

Hereisanexampleofcalculationoftheskinsurfaceofyourbodyfromheighth$ and weight
W using DuBois' formula:

>> h=182; W=71; >> A=h^0.725*W^0.425*71.84e-4 A = 1.9129

Name | Option | Description |
---|---|---|

format | short | displayformat for floatingpoint numbers: 5 digits, all digits or in BASE with COUNT digits. |

- | long | - |

- | BASE COUNT | - |

syms | | declare |

clear | | delete variables |

who | - | |

path | - | show searchpath |

addpath | path | add to searchpath |

hold | - | lock/unlock graphic window |

hold | on | lock graphic window |

hold | off | unlock graphic window |

Type | Example | O-Mode | Example M-Mode |
---|---|---|---|

Number | -3.214e-12 | -3.214e-12 | |

Vector | [1 2 3 4] | [1,2,3,4] | |

Matrix | [1 2 ; 3 4] | matrix([1,2],[3,4]) | |

Vectorelement | x(2) | x[2] | |

Matrixelement | x(2,1) | x[2,1] | |

Range | 2:5 | - | |

Variable | x = 2.321 | x : 2.321; |

This table shows operators with one example each for Octave-mode (O-mode) and Maxima-mode (M-mode). Operations are executed according to their precedence (prec), and, if these are equal, in order.

Function | O-Mode | M-Mode | Prec | Order |
---|---|---|---|---|

Addition | x + y | x + y | 4 | left-right |

Subtraction | x - y | x - y | 4 | left-right |

Scalar Multiplication | x .* y | x * y | 3 | left-right |

Vector/Matrix Multiplication | x * y | x . y | 3 | left-right |

Scalar Division | x ./ y | x / y | 3 | left-right |

Matrix Division (right) | x / y | - | 3 | left-right |

Matrix Division (left) | x \ y | - | 3 | left-right |

Scalar Exponentiation | x .^ y | x ^ y | 1 | left-right |

Vector/Matrix Exponentiation | x ^ y | x ^^ y | 1 | left-right |

Range | x:y:z | - | 5 | left-right |

Assignment | x = z | x : y | 10 | right-left |

Assignment | x += z | - | 10 | right-left |

Assignment | x -= z | - | 10 | right-left |

Assignment | x /= z | - | 10 | right-left |

Assignment | x *= z | - | 10 | |

Preincrement | ++x | ++x | 10 | left-right |

Predecrement | –x | –x | 10 | left-right |

Postincrement | x++ | x++ | 10 | right-left |

Postdecrement | x– | x– | 10 | right-left |

Adjunct | x' | - | 1 | right-left |

Factorial | x! | - | 1 | left-right |

Function | O-Mode | M-Mode | Prec | Order |
---|---|---|---|---|

Less | x < y | x < y | 6 | left-right |

Less or equal | x ? y | x ? y | 6 | left-right |

Larger | x > y | x > y | 6 | left-right |

Larger or equal | x >= y | x >= y | 6 | left-right |

Equal | x == y | x == y | 6 | left-right |

Not equal | x ~= y | x ~= y | 6 | left-right |

And | x & y | x & y | 7 | left-right |

Or | x | y | x | y | 9 | left-right |

Not | ~x | ~x | 8 | left-right |

diff(function,x) differentiates function with respect to the symbolic variable x. The main variable of function is used if x is not provided. Functions defined by user programs can often be handled as well.

>> syms a,x
>> diff(a*x^3)
ans = 3*a*x^2
>> diff(a*x^3,a)
ans = x^3
>> diff(3*sqrt(exp(x)+2),x)
ans = 1.5*exp(x)/sqrt(exp(x)+2)
>> diff(sin(x)) % no variable specified
ans = 1 % use z=sin(x) as variable
>> diff(sin(x),x) % more reasonable
ans = cos(x)
>> function y=ttwo(x) y=2*x; end
>> diff(ttwo(sin(x)),x)
ans = 2*cos(x)

taylor(function, x, x0, n) calculates the n-th Taylorpolynomial in the symbolic variable x at the point x0.

>> syms x
>> taylor(log(x),x,1,1)
ans = x-1
>> rat( taylor(exp(x),x,0,6))
ans = 1/720*x^6+1/120*x^5+1/24*x^4+1/6*x^3+1/2*x^2+x+1
>> float( taylor(x^3*sin(2*x+pi/4),x,pi/8,2))
ans = 1.057*x^2-0.36751*x+4.1881E-2

integrate(function, x) integrates expression function with respect to the symbolic variable x. Jasymca uses the following strategy: Integrals of builtin-functions and all polynomials are provided:

>> syms x
>> integrate(x^2+x-3,x)
ans = 0.33333*x^3+0.5*x^2-3*x
>> integrate(sin(x),x)
ans = -cos(x)

If function is rational (i.e. quotient of two polynomials, whose coefficients do not depend on x)) we use the standard approach: Separate a polynomial part, then separate a square free part using Horowitz' [11] method, and finally integrate the rest using partial fractions. The final terms are collected to avoid complex expressions.

>> syms x
>> y=(x^3+2*x^2-x+1)/((x+i)*(x-i)*(x+3))
y = (x^3+2*x^2-x+1)/(x^3+3*x^2+x+3)
>> integrate(y,x)
ans = -1/4*log(x^2+1)+(-1/2*log(x+3)+(-1/2*atan(x)+x))
>> diff(ans,x) % control
ans = (x^3+2*x^2-x+1)/(x^3+3*x^2+x+3)

Expressions of type g(f(x))?f?(x) and f?(x)f(x) are detected:

>> syms x
>> integrate(x*exp(-2*x^2),x)
ans = -0.25*exp(-2*x^2)
>> integrate(exp(x)/(3+exp(x)),x)
ans = log(exp(x)+3)

Substitutions of type (a?x+b) are applied:

>> syms x
>> integrate(3*sin(2*x-4),x)
ans = -1.5*cos(2*x-4)

Products polynomial(x)?f(x) are fed through partial integration. This solves all cases where f is one of exp, sin , cos , log , atan.

>> syms x
>> integrate(x^3*exp(-2*x),x)
ans = (-0.5*x^3-0.75*x^2-0.75*x-0.375)*exp(-2*x)
>> integrate(x^2*log(x),x)
ans = 0.33333*x^3*log(x)-0.11111*x^3

All trig and exp-functions are normalized. This solves any expression, which is the product of any number and any type of exponentials and trigonometric functions, and some cases of rational expressions of trig- and exp-functions

>> syms x
>> integrate(sin(x)*cos(3*x)^2,x)
ans = -3.5714E-2*cos(7*x)+(5.0E-2*cos(5*x)-0.5*cos(x))
>> integrate(1/(sin(3*x)+1),x)
ans = -2/3*cos(3/2*x)/(sin(3/2*x)+cos(3/2*x))

The special case ax2+bx+c??????????? is implemented:

>> syms x
>> integrate(sqrt(x^2-1),x)
ans = 0.5*x*sqrt(x^2-1)-0.5*log(2*sqrt(x^2-1)+2*x)

Clever substitutions may be supplied manually through subst(). If all fails, integrate numerically using quad or romberg. The symbolic variable may be omitted if it is the main variable of expression. Integrations can be quickly verified using diff() on the result.

Two routines are supplied for numerical integration, which are both not very sophisticated. quad('expression',ll,ul) is modelled after the Octave/Matlab integration function, but much simpler. Simpson's method is applied with a fixed number of nodes. This function uses the ?eval?-method rather than symbolic variables. The function has to be supplied as quoted string, and must be compatible with vector arguments. Finally, the variable name must be x. romberg uses symbolic function definitions, and a symbolic variable has to be supplied. The maximum number of iterations is set by the variable rombergit (default 11) and accuracy by rombergtol (default: 10?4).

>> quad('exp(-x.^2)',0,5)
s = 0.88623
>> syms x
>> romberg(exp(-x^2),x,0,5)
ans = 0.88623

ode(expression,y,x) solves the linear first-order differential equation y?=f(x)?y+g(x). expression is the complete right-hand-side of the equation, x and y are symbolic variables. Free constants in the solution are marked C.

>> ode(x,y,x)
ans = 0.5*x^2+C
>> syms k
>> ode(-k*y,y,x)
ans = C*exp(-k*x)
>> ode(y*tan(x)+cos(x),y,x)
ans = (0.5*cos(x)*sin(x)+(0.5*x+C))/cos(x)

The basic arithmetic operations are marked with the usual symbols (+ - * / ) . Exponention is performed with the accent character (^). Multiplication and division precede addition and subtraction; any order of evaluation can be forced by parenthesis.

>> 3.23*(14-2^5)/(15-(3^3-2^3)) ans = 14.535 >> 4.5e-23/0.0000013 ans = 3.4615E-17 >> 17.4^((3-2.13^1.2)^0.16) ans = 13.125 >> 17.23e4/(1.12-17.23e4/(1.12-17.23e4/1.12)) ans = 76919

In addition to these arithmetic operators Jasymca provides operators for comparing numbers

< > >= <= == ~=

and for boolean functions

& | ~

. Logical true is the number 1, false is 0.

>> 1+eps>1 ans = 1 >> 1+eps/2>1 % defines eps ans = 0 >> A=1;B=1;C=1; % semicolon suppresses output. >> !(A&B)|(B&C) == (C~=A) ans = 1

The most common implemented functions are the squareroot (sqrt(x)), the trigonometric functions (sin(x), cos(x), tan(x)) and inverses (atan(x), atan2(y,x)), and the hyperbolic functions (exp(x), log(x)). A large number of additional functions are available, see the list in chapter 4. Some functions are specific to integers, and also work with arbitrary large numbers: primes(Z) expands Z into primefactors, factorial(Z) calculates the factorial function. Modular division is provided by divide and treated later in the context of polynomials.

>> log(sqrt(854)) % natural logarithm ans = 3.375 >> 0.5*log(854) ans = 3.375 >> float(sin(pi/2)) % argument in radian ans = 1 >> gammaln(1234) % log( gamma( x ) ) ans = 7547 >> primes(1000000000000000001) ans = [ 101 9901 999999000001 ] >> factorial(35) ans = 1.0333E40 >> factorial(rat(35)) % to make it exact. ans = 10333147966386144929666651337523200000000

Name(Arguments) | Function | Mod |
---|---|---|

float( | | M,O |

rat( | | M,O |

realpart( | realpart of | M,O |

imagpart( | imaginary part of | M,O |

abs( | absolute value of | M,O |

sign( | sign of | M,O |

conj( | | M,O |

angle( | angle of | M,O |

cfs( | continued fraction expansion of | M,O |

primes(VAR) | VAR decomposed into primes | M,O |

Name(Arguments) | Function | Mod |
---|---|---|

sqrt( | squareroot | M,O |

exp( | exponential | M,O |

log( | natural logarithm | M,O |

sinh( | hyperbolic sine | O |

cosh( | hyperbolic cosine | O |

asinh( | hyperbolic areasine | O |

acosh( | hyperbolic areacosine | O |

sech( | hyperbolic secans | O |

csch( | hyperbolic cosecans | O |

asech( | hyperbolic areasecans | O |

acsch( | hyperbolic areacosecans | O |

sin( | sine (radian) | M,O |

cos( | cosine (radian) | M,O |

tan( | tangens (radian) | M,O |

asin( | arcsine (radian) | M,O |

acos( | arccosine (radian) | M,O |

atan( | arctangens (radian) | M,O |

atan2( | arctangens (radian) | M,O |

sec( | secans (radian) | O |

csc( | cosecans (radian) | O |

asec( | arcsecans (radian) | O |

acsc( | arccosecans (radian) | O |

factorial(N) | factorial | M,O |

nchoosek(N,K) | binomial coefficient | O |

gamma( | gammafunction | M,O |

gammaln( | logarithm of gammafunction | M,O |

Parts of an expression may be replaced by other expressions using subst(a,b,c): a is substituted for b in c. This is a powerful function with many uses. First, it may be used to insert numbers for variables, in the example 3 for x in der formula 2x??e?x2.

>> syms x >> a=2*sqrt(x)*exp(-x^2); >> subst(3,x,a) ans = 4.275E-4

Second, one can replace a symbolic variable by a complex term. The expression is automatically updated to the canonical format. In the following example z3+2 is inserted for x in x3+2x2+x+7.

>> syms x,z >> p=x^3+2*x^2+x+7; >> subst(z^3+2,x,p) ans = z^9+8*z^6+21*z^3+25

Finally, the term b itself may be a complex expression (in the example z2+1). Jasymca then tries to identify this expression in c (example: z?x3z2+1?). This is accomplished by solving the equation a=b for the symbolic variable in b (example: z), and inserting the solution in c. This does not always succeed, or there may be several solutions, which are returned as a vector.

>> syms x,y,z >> c=x^3*z/sqrt(z^2+1); >> d=subst(y,z^2+1,c) d = [ x^3*sqrt(y-1)/sqrt(sqrt(y-1)^2+1) -x^3*sqrt(y-1)/sqrt(sqrt(y-1)^2+1) ] >> d=trigrat(d) d = [ x^3*sqrt(y-1)/sqrt(y) -x^3*sqrt(y-1)/sqrt(y) ]

The function trigrat(expression) applies a series of algorithms to expression.

- All numbers are transformed to exact format.
- Trigonometric functions are expanded to complex exponentials.
- Addition theorems for the exponentials are applied.
- Square roots are calculated and collected.
- Complex exponentials are backtransformed to trigonometric functions.

It is often required to apply float(expression) to the final result.

>> syms x >> trigrat(sin(x)^2+cos(x)^2) ans = 1 >> b=sin(x)^2+sin(x+2*pi/3)^2+sin(x+4*pi/3)^2; >> trigrat(b) ans = 3/2 >> trigrat(i/2*log(x+i*pi)) ans = 1/4*i*log(x^2+pi^2)+(1/2*atan(x/pi)-1/4*pi) >> trigrat(sin((x+y)/2)*cos((x-y)/2)) ans = 1/2*sin(y)+1/2*sin(x) >> trigrat(sqrt(4*y^2+4*x*y-4*y+x^2-2*x+1)) ans = y+(1/2*x-1/2)

trigexpand(expression) expands trigonometric expressions to complex exponentials. It is the first step of the function trigrat above.

>> syms x >> trigexp(i*tan(i*x)) ans = (-exp(2*x)+1)/(exp(2*x)+1) >> trigexp(atan(1-x^2)) ans = -1/2*i*log((-x^2+(1-1*i))/(x^2+(-1-1*i)))

These data types are either used for multidimensional objects, or for simultaneous calculations on large numbers of data, e.g. for statistical problems. In this chapter we discuss this latter aspect. Vectors are marked with square brackets. The elements are entered as comma-separated list. The commas may be left if the elements can be distinguished in a unique manner, which however fails in the second example below:

>> x=[1,-2,3,-4] x = [ 1 -2 3 -4 ] >> x=[1 - 2 3 -4] % Caution: 1-2=-1 x = [ -1 3 -4 ]

Colon and the function line space are used to define ranges of numbers as vectors.

>> y=1:10 % 1 to 10, step 1 y = [ 1 2 3 4 5 6 7 8 9 10 ] >> y=1:0.1:1.5 % 1 to 1.5, step 0.1 y = [ 1 1.1 1.2 1.3 1.4 1.5 ] >> y=linspace(0,2,5) % 5 from 0 to 2.5, equidistant. y = [ 0 0.5 1 1.5 2 ]

The number of elements in a vector x is calculated with the function length(x), individual elements are extracted by providing the index k like x(k). This index k must be a number in the range 1 to (including) length(x). The colon operator plays a special role: Used as index, all elements of the vector are returned. Additionally, ranges of numbers can be used as index.

>> y(2) % single element ans = 0.5 >> y(:) % magic colon ans = [ 0 0.5 1 1.5 2 ] >> y(2:3) % index between 2 and 3 ans = [ 0.5 1 ] >> y(2:length(y)) % all from index 2 ans = [ 0.5 1 1.5 2 ] >> y([1,3,4]) % indices 1,3 and 4 ans = [ 0 1 1.5 ] >> y([1,3,4]) = 9 % insert ans = [ 9 0.5 9 9 2 ] >> y([1,3,4]) = [1,2,3] % insert ans = [ 1 0.5 2 3 2 ]

Matrices are handled in a similar way, only with two indices for rownumber (first index) and columnnumber (second index). Rows are separated by either a semicolon or a linefeed during input.

>> M=[1:3 ; 4:6 ; 7:9] M = 1 2 3 4 5 6 7 8 9 >> M([1 3],:) ans = 1 2 3 7 8 9 >> C=M<4 C = 1 1 1 0 0 0 0 0 0

The operators of chapter 2.2 may be applied to vectors and matrices. If scalar, per-element operation is desired, some operators (* / ^) must be preceded by a point to distinguish them from the quite different linear-algebra versions of these operations (see chapter 2.9). Further useful functions are sum(vector) and prod(vector) which return the sum and product of the vectors elements.

Name(Arguments) | Function | Mod | Ref |
---|---|---|---|

linspace( | vector with COUNT numbers ranging from | O | |

length( | number of elements in | M,O | |

zeros(ROWS[,COLUMNS]) | matrix of zeros | M,O | |

ones(ROWS[,COLUMNS]) | matrix of ones | M,O | |

eye(ROWS[,COLUMNS]) | matrix with diagonal one | M,O | |

rand(ROWS[,COLUMNS]) | matrix of random numbers | M,O | |

hilb(RANK) | Hilbertmatrix | M,O | |

invhilb(RANK) | Inverse Hilbertmatrix | O | |

size( | number of rows and columns | M,O | |

sum( | if | M,O | |

find( | indices of nonvanishing elements | M,O | |

max( | largest element in | M,O | |

min( | smallest element in | M,O | |

diag( | if | M,O | |

det( | determinante | M,O | |

eig( | eigenvalues | M,O | |

inv( | inverse | M,O | |

pinv( | pseudoinverse | M,O | |

lu( | LU-decomposition | M,O | |

svd( | singular value decomposition (Lapack) | M,O | |

qr( | QR-decomposition (Lapack) | M,O | |

eigen( | eigenvalues (Lapack) | O |

$\backslash sum\_\{k=1\}^\{k=n\}k$ $\backslash sum\_\{k=1\}^\{k=n\}k^2$$\backslash sum\_\{k=1\}^\{k=n\}k^3$ for n=10, n=100, n=10000. The solution is given below:

>> n=10; k=1:n; >> sum(k), sum(k.*k), sum(k.^3) % point! ans = 55 ans = 385 ans = 3025 >> n=100; k=1:n; >> sum(k), sum(k.*k), sum(k.^3) ans = 5050 ans = 3.3835E5 ans = 2.5503E7 >> n=10000; k=1:n; >> sum(k), sum(k.*k), sum(k.^3) ans = 5.0005E7 ans = 3.3338E11 ans = 2.5005E15

Several standardmatrices are created by means of functions without specifying individual elements: ones(n,m), zeros(n,m), rand(n,m) return matrices with elements 1, 0 or random numbers between 0 and 1. eye(n,m) has diagonalelements 1, else 0, and hilb(n) creates the n-th degree Hilbert-matrix.

>> A=rand(1,3) A = 0.33138 0.94928 0.56824 >> B=hilb(4) B = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 1/4 1/5 1/6 1/7

The following functions are provided for matrix calculations: diag(x) (extracts diagonal elements), det(x) (determinante), eig(x) (eigenvalues), inv(x) (inverse), pinv(x) (pseudoinverse). The adjunct matrix is created using the operator '.

>> det(hilb(4)) ans = 1/6048000 >> M=[2 3 1; 4 4 5; 2 9 3]; >> M' ans = 2 4 2 3 4 9 1 5 3 >> eig(M) ans = [ 11.531 -3.593 1.062 ] >> inv(M) ans = 0.75 0 -0.25 4.5455E-2 -9.0909E-2 0.13636 -0.63636 0.27273 9.0909E-2

The nontrivial functions are all based on the LU-decomposition, which is also accessible as a function call lu(x). It has 2 or 3 return values, therefor the left side of the equation must provide multiple variables, see example below:

>> M=[2 3 1; 4 4 5; 2 9 3] >> [l,u,p]=lu(M) % 2 or 3 return values l = % left triangular matrix (perm.) 0.5 0.14286 1 1 0 0 0.5 1 0 u = % right upper triangular matrix 4 4 5 0 7 0.5 0 0 -1.5714 p = % permutation matrix 0 0 1 1 0 0 0 1 0

Without preceding point the arithmetic operators function as matrix operators, e.g. * corresponds to matrix and vector multiplication.

>> x=[2,1,4]; y=[3,5,6]; >> x.*y % with point ans = [ 6 5 24 ] >> x*y % without point ans = 35

If one of the arguments is a scalar datatype, the operation is repeated for each element of the other argument:

>> x=[2,1,4]; >> x+3 ans = [ 5 4 7 ]

Matrix division corresponds to multiplication by the pseudoinverse. Using the operator \ leads to left-division, which can be used to solve systems of linear equations:

>> M=[2 3 1; 4 4 5; 2 9 3]; >> b=[0;3;1]; >> x=M\b % solution of M*x = b x = -0.25 -0.13636 0.90909 >> M*x % control ans = 0 3 1

Systems of linear equations can (and should) be solved directly with the function linsolve(A,b).

The application jasymca (not the applet or midlet) contain JLAPACK [9], the Java-port of the LAPACK [10]-routines with extended and better algorithms for matrix calculations. However, these are limited to matrices with real coefficients in floating point format. The LAPACK routines are accessed by the following functions:

%Singular value decomposition of A (1 or 3 returnvalues).

>> A=[2 3 1; 4 4 5; 2 9 3]; >> svd(A) ans = [ 12.263 3.697 0.9705 ]

%QR-decomposition of A (2 returnvalues).

>> A=[2 3 1; 4 4 5; 2 9 3]; >> [q,r]=qr(A) q = -0.40825 -5.3149E-2 -0.91132 -0.8165 -0.4252 0.39057 -0.40825 0.90354 0.13019 r = -4.899 -8.165 -5.7155 0 6.2716 0.53149 0 0 1.4321

Solves $A\cdot x=b$ (1 returnvalue). Example in chapter 2.13.1.

Solves $A\cdot x=b$, overdetermined (1 return value). For an example see insert ``Comparison of LAPACK and Jasymca Routines''.

% Eigenvalues of A (1 returnvalue).

>> A=[2 3 1; 4 4 5; 2 9 3]; >> eigen(A) ans = [ 11.531 1.062 -3.593 ]

Comparison of LAPACK and Jasymca Routines We calculate the 4-th degree regression polynomial for the following x,y-data:

>> x=[1:6],y=x+1 x = [ 1 2 3 4 5 6 ] y = [ 2 3 4 5 6 7 ] >> polyfit(x,y,4) p = 5.1958E-14 -9.6634E-13 -2.4727E-12 1 1

The coefficients p(1),p(2),p(3) should vanish since x and y represent a perfect straight line. This is an unstable problem, and it can be easily extended to make Jasymca completely fail. In our second attempt we use the Lapack-routine linlstsq:

>> x=[1:6],y=x+1; >> l=length(x);n=4; >> X=(x'*ones(1,n+1)).^(ones(l,1)*(n:-1:0)) >> linlstsq(X,y') ans = -1.6288E-18 -7.0249E-17 1.0653E-15 1 1

The coefficients p(1),p(2),p(3) are now significantly smaller. This particular problem can be solved exactly using Jasymca-routines and exact numbers, which avoids any rounding errors:

>> x=rat([1:6]);y=x+1; >> polyfit(x,y,4) p = 0 0 0 1 1

Solving systems of linear equations is accomplished by either the function linsolve(A,b) (all versions of Jasymca), or linsolve2(A,b) (LAPACK, not in applet and midlet). In both cases A is the quadratic matrix of the system of equations, and b a (row or column) vector representing the right-hand-side of the equations. The equations may be written as A?z=b and we solve for z.

>> A=[2 3 1; 4 4 5; 2 9 3]; >> b=[0;3;1]; >> linsolve(A,b) ans = -0.25 -0.13636 0.90909 >> linsolve2(A,b) % not Applet or Midlet ans = -0.25 -0.13636 0.90909

For large numeric matrices one should use the LAPACK-version if available. The Jasymca version can also handle matrices containing exact or symbolic elements. To avoid rounding errors in these cases it is advisable to work with exact numbers if possible:

>> syms x,y >> A=[x,1,-2,-2,0;1 2 3*y 4 5;1 2 2 0 1;9 1 6 0 -1;0 0 1 0] A = x 1 -2 -2 0 % symbolic element 1 2 3*y 4 5 % symbolic element 1 2 2 0 1 9 1 6 0 -1 0 0 1 0 0 >> b = [1 -2 3 2 4 ]; >> trigrat( linsolve( rat(A), b) ) ans = (-6*y-13/2)/(x+8) (20*y+(-9*x-151/3))/(x+8) 4 ((-3*x+10)*y+(-49/4*x-367/6))/(x+8) (-34*y+(13*x+403/6))/(x+8)

?Equation? in the following means the equation expression = 0. Equations are solved for a symbolic variable x by the function solve(expression, x). If expression is a quotient, then nominator = 0 is solved. Jasymca uses the following strategy to solve equations: Unordered List ItemFirst, all occurrences of the variable x in expression are counted, both as free variable and embedded inside functions. Example: In x3?sin(x)+2x2?x?1????? x occurs three times: as free variable, in sin(x) and in x?1?????. Unordered List ItemIf this count is one, then we are dealing with a polynomic equation, which is solved for the polynomial's main variable, e.g. z. This works always, if the polynomial's degree is 2 or of it is biquadratic, otherwise only, if the coefficients are constant. In the next step the solution is solved for the desired variable x. As an example: Jasymca has to solve sin2(x)?2sin(x)+1=0 for x. It first solves z2?2z+1=0 for z and then sin(x)=z for x. Examples with free variables:

>> syms x,b >> solve(x^2-1,x) ans = [ 1 -1 ] >> solve(x^2-2*x*b+b^2,x) ans = b

An example with function variable (exp(j?x)):

>> syms x >> float( solve(sin(x)^2+2*cos(x)-0.5,x) ) ans = [ 1.438i -1.438i -1.7975 1.7975 ]

If count is 2, only one case is further considered: The variable occurs free and inside squareroot. This squareroot is then isolated, the equation squared and solved. This case leads to additional false solutions, which have to be sorted out manually.

>> syms x >> y=x^2+3*x-17*sqrt(3*x^2+12); >> solve(y,x) ans = [ -32.501 26.528 -1.3931E-2-2.0055i -1.3931E-2+2.0055i ]

In all other cases Jasymca gives up.

Coupled systems of equations can be solved by the function algsys([expressions],[symbolic variables]). First, all linear equations are solved using the Gauss-method, then each equation is fed through solve() and the solution used to eliminate one variable in all other expressions. The equations are treated in the order they are supplied. This method only works for simple systems. The solution is provided as vector of solution vectors, each individual solution in as linear factor: In the first example below there is one solution with xs=-2/3, a2=3/4, a0=2, a1=0, the second example has two solutions.

>> syms xs,a0,a1,a2 >> algsys([2-a0,a1-0,a2*xs^2+a1*xs+a0-3-xs, 2*a2*xs+a1+1],[a2,a1,a0,xs]) ans = [ [ xs+2/3 a2-3/4 a0-2 a1 ] ] >> syms a,xs >> algsys([a*xs+3*a-(3-xs^2),a+2*xs],[a,xs]) ans = [ [ -sqrt(6)+(xs+3) 2*sqrt(6)+(a-6) ] [ sqrt(6)+(xs+3) -2*sqrt(6)+(a-6) ] ] >> float(ans) ans = [[ xs+0.55051 a-1.101 ] [ xs+5.4495 a-10.899 ]]

Jasymca can handle polynomials with symbolic variables. In this chapter, however, we work with the Matlab/Octave/Scilab-approach of using vectors as list of polynomial coefficients: A polynomial of degree n is represented by a vector having n+1 elements, the element with index 1 being the coefficient of the highest exponent in the polynomial. With poly(x) a normal polynomial is created, whose zeros are the elements of x, polyval(a,x) returns function values of the polynomial with coefficients a in the point x, roots(a) calculates the zeros, and polyfit(x,y,n) calculates the coefficients of the polynomial of degree n, whose graph passes through the points x and y. If their number is larger than n+1 a least square estimate is performed. The regression analysis in exercise 8 was performed using this method. The roots of a 4th-degree polynomial are -4,-2,2,4 and it intersects the y-axis at (y=-64). Calculate its coefficients:

>> a=poly([-4,-2,2,4]) a = [ 1 0 -20 0 64 ] >> a = -64/polyval(a,0) * a a = -1 0 20 0 -64 >> a = polyfit([-4,-2,2,4,0],[0,0,0,0,-64],4) a = -1 0 20 0 -64

Example 2: On a grid with x,y-Koordinaten we have the following greyvalues of a digital image:

y\x 1 2 3 4 1 98 110 122 136 2 91 112 131 141 3 73 118 145 190 4 43 129 170 230

Which greyvalue do you expect at position x=2,35;y=2,74? Calculate the bicubic interpolation.

Solution:

>> Z=[98,110,122,136; > 91,112,131,141; > 73,118,145,190; > 43,129,170,230]; >> i=1; p=polyfit(1:4,Z(i,:),3); z(i)=polyval(p,2.35); >> i=2; p=polyfit(1:4,Z(i,:),3); z(i)=polyval(p,2.35); >> i=3; p=polyfit(1:4,Z(i,:),3); z(i)=polyval(p,2.35); >> i=4; p=polyfit(1:4,Z(i,:),3); z(i)=polyval(p,2.35); >> p=polyfit(1:4,z,3); zp=polyval(p,2.74) zp = 124.82

We have learnt that polynomials may be represented by the vector of their coefficient. Using a symbolic variable x we will now create a symbolic polynomial p. Conversely, we can extract the coefficients from a symbolic polynomial using the function coeff(p, x, exponent). The command allroots(p) returns the zeros.

>> a=[3 2 5 7 4]; % coefficients >> syms x >> y=polyval(a,x) % symbolic polynomial y = 3*x^4+2*x^3+5*x^2+7*x+4 >> coeff(y,x,3) % get one coefficient ans = 2 >> b=coeff(y,x,4:-1:0) % or all at once b = [ 3 2 5 7 4 ] >> allroots(y) % same as roots(a) ans = [ 0.363-1.374i 0.363+1.374i -0.697-0.418i -0.697+0.418i ]

Up to this point there is little advantage of using symbolic calculations, it is just another way of specifying a problem. The main benefit of symbolic calculations emerges when we are dealing with more than one symbolic variable, or, meaning essentially the same, when our polynomial has nonconstant coefficients. This case can be treated efficiently only with symbolic variables. Notice in the example below how the polynomial y is automatically multiplied through, and brought into a canonical form. In this form the symbolic variables are sorted alphabetically, i.e. z is main variable compared to x. The coefficients can be calculated for each variable separately.

>> syms x,z >> y=(x-3)*(x-1)*(z-2)*(z+1) y = (x^2-4*x+3)*z^2+(-x^2+4*x-3)*z+(-2*x^2+8*x-6) >> coeff(y,x,2) ans = z^2-z-2 >> coeff(y,z,2) ans = x^2-4*x+3

The command allroots functions with variable coefficients also, but only, if the polynomials degree in the main variable is smaller than 3, or it is biquadratic. If roots of other variables x are searched, one should use the more general solve(p,x), which will be discussed in more detail later.

>> syms x,z >> y = x*z^2-3*x*z+(2*x+1); >> allroots(y) ans = [ sqrt((1/4*x-1)/x)+3/2 -sqrt((1/4*x-1)/x)+3/2 ] >> solve(y,x) ans = -1/(z^2-3*z+2)

The decomposition of p in linear, quadratic, cubic etc factors is accomplished by sqfr(p). Returned is a vector of factors sorted in ascending order of the exponents.

>> syms x >> y=(x-1)^3*(x-2)^2*(x-3)*(x-4) y = x^7-14*x^6+80*x^5-242*x^4+419*x^3-416*x^2+220*x-48 >> z=sqfr(y) z = [ x^2-7*x+12 x-2 x-1 ]

The division of two polynomials p and q in one polynomial and remainder is calculated using divide(p,q). If the polynomials have more than one variable, an optional variable can be specified, which will be used for division. gcd(p,q) returns the greatest common denominator of two expressions. Both functions also work with numbers as arguments.

>> divide(122344,7623) ans = [ 16 376 ] >> divide(2+i,3+2*i) ans = [ 1/2 1/2 ] >> syms x,z >> divide(x^3*z-1,x*z-x,x) ans = [ x^2*z/(z-1) -1 ] >> divide(x^3*z-1,x*z-x,z) ans = [ x^2 x^3-1 ] >> gcd(32897397,24552502) ans = 377 >> gcd(z*x^5-z,x^2-2*x+1) ans = x-1

realpart(expression) and imagpart(expression) is used to decompose complex expressions. Symbolic variables in these expressions are assumed to be real-valued.

>> syms x >> y=(3+i*x)/(2-i*x) y = (-x+3i)/(x+2i) >> realpart(y) ans = (-x^2+6)/(x^2+4) >> imagpart(y) ans = 5*x/(x^2+4)

When you are working in the prompt, the answers are given automatically. If you have chosen to print output by executing a macro file (with the extension ?m?), use the ?printf()? method as in this example:

x=log(sqrt(854)); % natural logarithm printf('Answer=%f\n', x);

You can output any symbolic variable as well. Here is a more complicated example where we print a value, symbolic equation and a matrix:

syms x,y; a=2*sqrt(x)*exp(-x^2); ff=subst(3,x,a); printf("%f",ff); ss=trigrat(sqrt(4*y^2+4*x*y-4*y+x^2-2*x+1)); printf("%f\n",ss); M=[1:3 ; 4:6 ; 7:9]; C=M<4; printf("%f\n",C);

Data may be graphed using the plot(x,y)-function, x and y being equalsized vectors, which denote the coordinates of the data points to be plotted. A third optional argument plot(x,y,option) specifies plot options like colors and symbols, see exercise 7ff. The graphic may then be decorated (axis, title) and exported as encapsulated-postscript file for inclusion in text documents. With hold (or hold on) the graphic gets locked so that subsequent plot commands use the same window. Repeating hold (or hold off) deletes the graphic. The method

plot(x,y [,option])

plots x versus y using option x and y are equalsized vectors. option is a string, specifying color (one of r,g,b,y,m,c,w,k) and symbol (one of +,*,o,x). Default: blue lines. Logarithmic and semilogarithmic plots are provided with the functions loglog, linlog and loglin. Here are the major plot methods:

Name(Arguments) | Function |
---|---|

plot(x,y [,option]) | Plot x versus y |

loglog(x,y[,option]) | logarithmic plot |

linlog(x,y[,option]) | semilogarithmic plot |

loglin(x,y[,option]) | semilogarithmic plot |

print(name ) | write graphic in eps-file |

Example 1 let us plot the function: $y=\backslash frac\{1\}\{1+2x^2\}$ in the range x=0.01?100 linear and logarithmic. The code is below:

>> x=0.01:0.01:100; y=1./(1+0.5*x.*x); plot(x,y) >> x=0.01:0.01:100; y=1./(1+0.5*x.*x); loglog(x,y)

Example 1a let us plot the function cos(x):

>> x=-10:0.1:10; >> plot(x,cos(x));

Example 2: Display of Lissajous-figures: From the vector t=0:0.1:4*pi; create the trigonometric expressions x=sin(0.5*t+1); and y=cos(1.5*t);. The plot x vs. y is called Lissajous-figure. Create different figures by variating the constants 0.5,1,1.5 in the definition. Partial solution:

>> t=0:0.1:4*pi; >> x=sin(0.5*t+1); >> y=cos(1.5*t); >> plot(x,y)

Example 3 Calculate the first 100 elements of the sequences $x\_n\; =\; \backslash frac\{n-1\}\{n\};\; x\_n=\backslash frac\{n+1\}\{n\};\; x\_n=\backslash frac\{n+(-1)^n\}\{n\}$ Plot xn versus n using the command plot(n, xn). Variate the plotoptions (colors, symbols).

>> n=1:100; >> x1=(n-1)./n; x2=(n+1)./n; x3=(n+(-1).^n)./n; >> plot(n,x1) % Standard: blue, lines >> hold Current plot held. >> plot(n,x2,'r') % Color: r,g,b,c,m,y,w. >> plot(n,x3,'g')

Example 4 Plot the datapoints of the following table using plot and colored symbols. Calculate the linear regression using polyfit, and plot the regression line in the same graph. Add title and labels, and export the graphic to a file suitable for inclusion in a text document.

x 0 1 2 3 4 5 6 7 8 9 y -3.1 -0.7 1.8 4.1 6.2 8.9 11.3 13.5 16 18.3

>> x=0:9; >> y=[-3.1,-0.7,1.8,4.1,6.2,8.9,11.3,13.5,16,18.3]; >> plot(x,y,"+r") % Symbol: o,x,+,* >> hold Current plot held. >> plot(x,polyval(polyfit(x,y,1),x)) % Regression >> xlabel("x-Achse") >> ylabel("y-Achse") >> title("Beispielgraph") >> print("graph.eps") % not Applet or Midlet!