Sign up

Jasymca

Login

Do math, code and ?chat on your phone tablet or laptop.

Contact us: info@holoklip.com

Quick Guide

Numbers

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);

Symbolic Variables

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

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    

Commands

Name Option Description
format short displayformat for floatingpoint numbers: 5 digits, all digits or in BASE with COUNT digits.
- long -
- BASE COUNT -
syms sym1 [,sym2,?] declare sym1,? as symbolic variable.
clear var1 [,var2,?] delete variables var1,?
who -
path - show searchpath
addpath path add to searchpath
hold - lock/unlock graphic window
hold on lock graphic window
hold off unlock graphic window
Define
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;
Calculate

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
Comparison and Logical Operators
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

Calculus

Differentiation

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)
Taylorpolynomial

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
Indefinite Integral

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.

Numerical Integration

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
Differential Equations

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)

Functions

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    
Scalar
Name(Arguments) Function Mod
float(var) var as floating point number M,O
rat(var) var as exact number M,O
realpart(var) realpart of var M,O
imagpart(var) imaginary part of var M,O
abs(var) absolute value of var M,O
sign(var) sign of var M,O
conj(var) var conjugate complex M,O
angle(var) angle of var M,O
cfs(var) [varT]) continued fraction expansion of var with accuracy varT M,O
primes(VAR) VAR decomposed into primes M,O
Scalar functions
Name(Arguments) Function Mod
sqrt(var) squareroot M,O
exp(var) exponential M,O
log(var) natural logarithm M,O
sinh(var) hyperbolic sine O
cosh(var) hyperbolic cosine O
asinh(var) hyperbolic areasine O
acosh(var) hyperbolic areacosine O
sech(var) hyperbolic secans O
csch(var) hyperbolic cosecans O
asech(var) hyperbolic areasecans O
acsch(var) hyperbolic areacosecans O
sin(var) sine (radian) M,O
cos(var) cosine (radian) M,O
tan(var) tangens (radian) M,O
asin(var) arcsine (radian) M,O
acos(var) arccosine (radian) M,O
atan(var) arctangens (radian) M,O
atan2(var1, var2) arctangens (radian) M,O
sec(var) secans (radian) O
csc(var) cosecans (radian) O
asec(var) arcsecans (radian) O
acsc(var) arccosecans (radian) O
factorial(N) factorial N! M,O
nchoosek(N,K) binomial coefficient (NK) O
gamma(var) gammafunction M,O
gammaln(var) logarithm of gammafunction M,O

Programming

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.

Branches

?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

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
Jumps

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

Symbolic Transformations

Substitution

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) ]
Simplifying and Collecting Expressions

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)))

Vectors and Matrices

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.

Matrix and Vector Operations
Name(Arguments) Function Mod Ref
linspace(var1,var2,COUNT) vector with COUNT numbers ranging from var1 to var2 O
length(vector) number of elements in vector 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(matrix) number of rows and columns M,O
sum(var) if var is a vector: sum of elements, if var is a matrix: sum of columns M,O
find(var) indices of nonvanishing elements M,O
max(var) largest element in var M,O
min(var) smallest element in var M,O
diag(var,[OFFSET]) if var is a vector: matrix with var as diagonale, if var is matrix: diagonale as vector M,O
det(matrix) determinante M,O
eig(matrix) eigenvalues M,O
inv(matrix) inverse M,O
pinv(matrix) pseudoinverse M,O
lu(matrix) LU-decomposition M,O
svd(matrix) singular value decomposition (Lapack) M,O
qr(matrix) QR-decomposition (Lapack) M,O
eigen(matrix) eigenvalues (Lapack) O
Working with vectors

\sum_{k=1}^{k=n}k \sum_{k=1}^{k=n}k^2 \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
Functions

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).

Lapack

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:

svd(A)
%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(A)
 %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
linsolve2( A, b)
Solves $A\cdot x=b$ (1 returnvalue). Example in chapter 2.13.1. 
linlstsq(A, b)
Solves $A\cdot x=b$, overdetermined (1 return value). For an example see insert ``Comparison of LAPACK and Jasymca Routines''. 
eigen(A)
% 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

Equations

Systems of Linear Equations

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)
Nonlinear Equations

?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.

Systems of Nonlinear Equations

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 ]]

Polynomials

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
Polynomial functions. Summary table

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
Roots

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)
Squarefree Decomposition

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 ]
Division and Common Denominator

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
Real- and Imaginary Part

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)

Printf output

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);
Plot

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=\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 = \frac{n-1}{n}; x_n=\frac{n+1}{n}; x_n=\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!
For More Help