Matlab Breitbach Ss08

G. Breitbach Math. Tools & Simulation SS2008

MATLAB / SIMULINK

Introduction into the MATLAB/SIMULINK Software Package 1. The MATLAB Environment 2. Numbers, Functions, Operations 3. Symbolic Calculations 4. Graphical Tools, Plotting of Functions 5. Vectors, Matrices and Linear Equations 6. Script Files 7. Differential Equations

Simulation Example A A Simple Dynamic Absorber

Simulation Example B The Motion of a Structure Consisting of Two Hinged Bars Carrying Two Masses m, M

Simulation Example C Energy Balance of a House

Breitbach, Math. Tools and Simulation, SS2008

1

Introduction into the M A T L A B/S I M U L I N K Software Package What is MATLAB? MATLAB is a program for computation and visualisation. It is available on computer systems ranging from personal to supercomputers. MATLAB ist controlled by commands, and it is programmable. MATLAB can be used together with other programs (e.g. FORTRAN, C etc.) Range of areas using MATLAB: Research and development in the industrial field Teaching mathematics, especially linear algebra and differential equations Research in numerical analysis and scientific computing Electronics Control theory Physics Economics Chemistry Biology The central object in MATLAB is the matrix (The name MATLAB is derived from MATrix LABoratory) and the operations with matrices. A matrix is an array of numbers arranged in rows and columns. A vector is a special matrix consisting of one row or one column only. A single number is a 1x1-matrix. There are Arithmetic operations Logical operations Mathematical functions Graphical tools I/O operations (data exchange) In principle all operations are related to matrices and vectors SIMULINK is based on MATLAB. It is designed with special graphical features for model based and system level analysis.

1. The MATLAB Environment Start MATLAB by clicking the MATLAB ikon on the desktop. The Command window is opened. The prompt >> is waiting for commands. Beside this window MATLAB has two other windows: The EDIT/DEBUGGER window and the figure window MOST IMPORTANT COMMAND FOR THE BEGINNER HAVING PRODUCED NONSENSE AND RUNNING IN TROUBLES: Breitbach, Math. Tools and Simulation, SS2008

2

Type simply exit informations are lost.

at the prompt. The program is closed and all useless

Instructions&Problems 1.1: It is possible to use typical DOS-commands like dir or cd (change directory). The directories connected with the MATLAB-program can be found with the dir and the cd command at the prompt. If we are working in a special directory, all files of that directory are listed with dir. Using the command cd ..(please observe that the blanc between the cd and the two dots are important!) changes to the next level directory. In the bin directory You will find a subdirectory work. Insert now cd work to go into the work directory. Again the dir command delivers eventually some files saved in work. It is a good practice to specify a working directory and to be there directly after start up. On the desktop clicking the MATLAB-ikon with the right mouse button and in „Eigenschaften“ insert the directory C:\matlab\work. Then calling MATLAB the current directory will always be C:\matlab\work. Instructions&Problems 1.2: MATLAB can be used like an ordinary calculator is used. Type 3+4 and you get the answer ans = 7 In case you create a statement too long for a single line use the e l l i p s i s (...) e.g. 1+2+3+4+5+6 ... (first line) +7+8+9 (second line) Try the other operations beside addition: - Subtraction; * Multiplication; / Division (and a “ \ Division“!); ^Exponentiation. Dealing with normal numbers the meaning of the operators is quite clear: 4+5 or 4*5 or 2^2 or 4/5. The operation 4\5 is more strange. Type it and see, what the result is: 1.25. 4\5 gives the same result as 5/4 ! We have in MATLAB two types of division: the right (that is the ordinary division) and the left division. It will be explained later more when dealing with matrix operations.

A statement like v = 15 creates a variable v stores the value 15 in it and saves it in a part of the computer memory (workspace). Define at the prompt: v=3 and after pressing enter insert z=4 (then enter again). Type then at the prompt v+z to obtain the result ans = 7. The predefined variable ans is assigned the value of the latest computed expression which is not given a name.

Breitbach, Math. Tools and Simulation, SS2008

3

Insert: v=21;z=17 in one line. z=17 is echoed. Calculate v+z to obtain 38. Both variables are overwritten by the new specified values. Now insert: v=50;z=70; there is no echo. The addition gives 130. Closing a command using the semicolon ; produces no echo. Insert: s=v*z (with or without semicolon at the end of the statement) and create the variable s with the value 3500. Hint: Using ↑ and ↓ keys at the prompt previous commands can be retrieved. Instructions&Problems 1.3: Commands to overview what we have done during a MATLAB-session. who and whos The „who“ command gives a list of the defined variables. The „whos“ command gives beside the listing informations about the size of the variables. MATLAB can keep a diary of what is displayed on the screen. Insert at the prompt diary filename (e.g. myname_today). Then perform some calculations or operations. Using diary off stops the recording. The resulting ASCII-file with name today can later be edited (hint: graphical output is not stored!). The stored values and results cannot by read by MATLAB on a later occasion. The command save filename v1 v2... stores the variables of the list in a file filename.mat in the current directory. Clear all and retrieve by load filename. With the command delete filename.mat the file filename.mat is deleted. Instructions&Problems 1.4: In case you want to delete variables of your actual session use the command: clear varname Varname denotes the variable which should be deleted. In case you use clear without any varnames all the variables are deleted. Be carefull!

2. Numbers, Functions, Operations Instructions&Problems 2.1: Some variables as ans are predefined. Here a list of other predefined variables: eps returns the machine accuracy (distance of 1 and the closest representable floating point number realmax returns the largest floating point number the computer can handle realmin returns the smallest non-zero floating point number the computer can handle pi returns π inf is defined as 1/0 i,j are defined as √-1 (imaginary unit) pi belongs to the predefined variables; type in pi to obtain ans = 3.1416. Breitbach, Math. Tools and Simulation, SS2008

4

pi is given here with 4 decimals. What about the representation of numbers? Results are usually displayed in a short floating point format with 4 decimals. Clicking f i l e and p r e f e r e n c e s... it is possible to select other representations with more decimals or the exponential formats.

Perform the following inputs: a=pi and then pi=10. Calculate 10*pi and 10*a. Then clear a and pi and type in the whos command. Type in a and then pi. Numbers may be complex numbers. The general form is a + b*i, where i is the imaginary unit with the property i*i = -1 (try it!) Commands related to complex numbers z: abs(z), angle(z) A polar representation of z can be plotted with compass(z). Instructions&Problems 2.2: There are build in function like sin, cos tan exp etc..Typing sin(1) yields the result: ans = 0.8415 Observe that the argument of the sin-function is interpreted as given in radians! Inserting sin(pi/2) gives ans = 1 Note that the natural logarithm is abreviated by log. What is the result of the expression log(exp(1))? Use the help command to find out more about the build in functions. There are directories elfun (elementary functions) and help specfun (special functions). Type help elfun or help specfun.

The functions may be used for complex arguments. Verify the famous Euler relation for some values of x: exp(i*x) = cos(x) + i*sin(x) (e.g. exp(i*pi) = -1) Calculate sqrt(-2). Try the logarithm of a negative number e.g. –5: log(-5); do the same for other numbers. Give an interpretation of the result. Instructions&Problems 2.3: A polynomial, p(x), of degree n is stored as a row of values ((n+1) entries) included by brackets: p(x) = an*xn + an-1*xn-1 +......+a0 ⇒ [an an-1 an-2 ..... a0]

Breitbach, Math. Tools and Simulation, SS2008

5

Define the polynomials p1 and p2: p1 = [3 -3 -6] and p2 = [5 -5 5 -5] What is the value of p1 in x=5? Type in y = polyval(p1,5) to obtain the result y = 54. Multiplication of two polynomials is performed via conv(p1,p2): p3 = conv(p1,p2) The result is p3 = [15 –30 0 0 -15 30] The roots of a polynomial p1 can be determined using the functions roots(p1). Find the roots of p1 and p2! The function polyder(p1) evaluates the derivative of p1 with respect to x. Try it for p1, p2, p3. The expression [q, r] = deconv(p1,p2) divides the polynomial p1 by p2. The quotient polynomial is represented by q and the remainder is r. Calculate p1/p2 and also p2/p1. The command [u v k] = residue(p1,p2) returns a partial fraction expansion: p1( x) u (1) u (i ) = + ... + + .... + k ( x) p 2( x) x − v(1) x − v(i ) The residues are returned in u, the pole locations in v and the quotient polynomial in k.

5x 3 − 5x 2 + 5x − 5 5 25 20 = x+ + Verify that 2 3 9( x − 2) 9( x + 1) 3x − 3x − 6 What is the result of

3x 2 − 3x − 6 ? 5x 3 − 5x 2 + 5x − 5

( =

1 . 2 x + 0 .6 − 0 .6 ) + x −1 x 2 +1

Instructions&Problems 2.4: MATLAB has relational operators: == equal to ~= not equal to > greater than < less than Like the C language, MATLAB does not have a logical or boolean data typ. A false result is connected with zero (0), a true result with a nonzero value (1 as a rule). Insert x = (5>4). The result is x=1. Typing x=(4>5) gives x=0. x is now a logical variable; in case that the expression is true the value is 1 otherwise 0. It is possible to check logical expressions. Define a = 0 and b = 0, then evaluate (a ==b); type c = sin(pi), then evaluate (a == c) The logical „Or“ is represented by | and the logical „and“ is represented by & the logical „not“ is represented by ~ . The logical exclusive „Or“ is expressed by xor

Breitbach, Math. Tools and Simulation, SS2008

6

Truth Tables Inp uts L1 L2 0 0 0 1 1 0 1 1

And L1 & L2 0 0 0 1

Or L1 | L2 0 1 1 1

Xor Xor(L1,L2) 0 1 1 0

Not ~ L1 1 1 0 0

Evaluate the following expressions: a) (4<3)|(10<15) b) ((4<3)|(~(10<15)))|(2>4) c) (5<7)&(7>3) d) (3>2)&(7<5)

Instructions&Problems 2.5: Working with sets: Define two sets of numbers: A=[1 2 3 4 5 6 7 8 9 10] B=[3 6 9 12] The following important operations can be performed on sets: union(A,B): the union of sets A, B is returned intersect(A,B): returns the intersection of A, B setdiff(A,B): returns the set elements belonging to A but not to B Use the operations for the two sets A,B given above (what about setdiff(A,B) and setdiff(B,A)?)

3. Symbolic Calculations Beside numerical analysis it is possible to perform symbolic calculations (fractions, algebraic operations, integration, differentiation etc.). Instructions&Problems 3.1: Calculate 1/2 + 1/3 ; then type at the prompt: sym(1/2 + 1/3).

π

16 + cos( ) 3 in the usual way, then using sym(expression). Compute the expression 3 8+4 Instructions&Problems 3.2: Application of factor and expand: Define symbolic variables a,b,c. They are defined via the commands a=sym(‘a‘);b=sym(‘b‘);c=sym(‘c‘); or abreviated syms a b c. Then evaluate (a + b + c)2 using expand(...).

Breitbach, Math. Tools and Simulation, SS2008

7

Use factor(...) to find a simpler form of the algebraic expressions a3 + 3a2b + 3 ab2 + b3 and a6 + a4 –a2 -1 Instructions&Problems 3.3:

Let us begin now with the integration procedure. At first define a symbolic variable x: >> x = sym(‘x‘) or >>syms x The integration command is int(f,x); determine the integral of f(t)=5*t*sin(t) 5x 3 − 5x 2 + 5x − 5 ∫ 3x 2 − 3x − 6 dx (see Instructions&Problems 2.3!)

Calculate

Use funtool to get a nice environment for evaluations and graphical representations of functions.

Instructions&Problems 3.4:

A function f(x) is given in the domain [-π, π] by f(x) = x². The function is periodic: f(x) = f(x + 2π). Use the tools of the symbolic toolbox to approximate f(x) by a Fourier series with 4 terms.

At first a short introduction into „Fourier series“: a g ( x) = 0 + ∑ al ⋅ cos(l ⋅ x) + bl ⋅ sin(l ⋅ x) is a periodic function with g(x) = g(x + 2π). It is 2 possible to adjust the coefficients a, b in such a way that an arbitrary function f(x) is sufficiently represented by sin/cos – series. The coefficients are determined via:

al =

1

π

π

∫

f ( x) ⋅ cos(l ⋅ x) ⋅ dx

and bl =

−π

1

π

π

∫ f ( x) ⋅ sin(l ⋅ x) ⋅ dx

−π

b

Use the command for the evaluation of definite integrals

∫ g ( x) ⋅ dx : a

>> int(g,x,a,b)

Taylor series: use the command syms x; taylor(sqrt(1+x),4,1) to obtain the Taylor series with 4 terms of the function 1 + x at point x=1. 4. Graphical Tools, Plotting of Functions Instructions&Problems 4.1:

We find in MATLAB a wide range of graphics commands which are flexible and easy to use. Breitbach, Math. Tools and Simulation, SS2008

8

The most simple command is plot(x,y) What about the quantities x and y? To plot a function y = f(x) it is necessary to establish a r r a y s. A row array can be defined as a series of numbers enclosed by brackets []: Type x = [-3 –1 1 2 4 6]; then y = [12 0 –4 –3 –5 21]; To plot use the command: plot(x,y) (or plot(y,x)). Type at the prompt>> x = 0:0.5:3.5; then y = x*sin(x) We will obtain an error message. x represents here an array of values and not soley a single number! The operations + - * / \ have in connection with arrays and matrices a special meaning. Type y = x.*sin(x). In case of elementwise arithmetic operations the classical arithmetic operations are + - .* ./ .\ .^. That is important when working with arrays as here done. Proceed by refining the curve: x=0:0.1:3.5; then y=x.*sin(x); then plot(x,y). Define an array z by: z=cos(x.*x); You can plot both curves y(x) and z(x) using the command: plot(x,y,x,z). Plot some of the functions available in MATLAB (e.g. Bessel-functions, besselj(1,x)) The axes can be scaled via the command axis([xmin, xmax, ymin, ymax]) if necessary. xlabel('text') writes the string text below the x-axis, in a similiar way works ylabel('text'). Draw a grid with the command grid on. Space curves in parametric form can be plotted via plot3(x1,x2,x3). x1(t), x2(t), x3(t) are functions of a variable t. Define x1(t) = cos(t), x2(t) = sin(t), x3(t) = t, then use plot3. Instructions&Problems 4.2:

Logarithmic and semi-logarithmic diagrams 1 f (x ) = 1 + 0 .2 ⋅ i ⋅ x − x 2 x=logspace(-1,1,100); y=1./(1 + 0.2.*i.*x - x.^2); loglog(x,abs(y)) semilogx(x,angle(y))

Use the subplot(n,m,i) - command to generate multiple plots on one screen.

5. Vectors, Matrices and Linear Equations Instructions&Problems 5.1: Generate a vector v having the components 0, 0.5, 1.0, 1.5, 2.0.

There are so called row vectors and column vectors. So we have two possibilities: Type v = [0 0.5 1.0 1.5 2.0] (a „row“ vector) Breitbach, Math. Tools and Simulation, SS2008

9

and w = [0; 0.5; 1.0; 1.5; 2.0] (a „column“ vector, numbers separated by semicolons) Of course a vector may consist of complex numbers. To change a row vector into a colum vector or vice versa use the apostrophe (prime): z=v‘ We have access to the components of a vector via commands like v(3). Using such command you can of course change the components. Operations on vectors: dot(a,b) dot (scalar) product of two vectors a, b cross(a,b) cross product of two vectors a,b (number of elements must be three!) A matrix is an array of numbers with n rows and m columns: Typing A = [1 3 2;3 6 7] creates a (2x3) matrix (The convention is that the first index n gives the total numbers of rows and the second the total numbers of columns ) Other possibilities (row by row): A=[1 3 2 3 6 7] Try B=[1 3 2 ... 3 6 7] or specify A(i,j) element by element. The first index i is the row index, the second index j gives the column. Instructions&Problems 5.2:

Elementary operations with matrices: Addition/Subtraction: A+B, A-B, only possible if coincident numbers of rows and columns respectively. The addition/subtraction is element by element. Define matrices A, B, C (A,B are (2x3)-matrices, C=(3x2)-matrix) ⎡1 2⎤ ⎡1 3 7⎤ ⎡ 4 2 9⎤ ⎥ ⎢ A= ⎢ ⎥ B = ⎢2 6 1⎥ C = ⎢9 8⎥ 5 1 6 ⎣ ⎦ ⎣ ⎦ ⎢⎣5 3⎥⎦ Calculate A+B, B-A, try A+C and then A+C‘ The apostrophe ‘ denotes the t r a n s p o s i t i o n (change of rows and columns): ⎡1 9 5⎤ C' = ⎢ ⎥ ⎣2 8 3⎦ For a square matrix is m=n (same number of row and columns)

Hint: In case that the matrix elements consists of complex numbers the transposed matrix contains the conjugate complex numbers! Create an arbitray matrix with complex entries and use the transposition operation. Multiplication The product matrix C of two matrices A, B (C = A*B) is defined in case that the number of columns of the first matrix A is equal to the number of rows of the second matrix B. Breitbach, Math. Tools and Simulation, SS2008

10

The element cij ist given by

ci j

=

n

∑a l =1

il

⋅ bl j . n denotes here the number of

columns of A and also the number of rows of B. Calculate A*C; try A*B and A*B’ Some functions on matrices: rank(A) gives the rank of A, that is the number of linearly independent rows and columns det(A) gives the determinant of the matrix A (A must be a square matrix!) inv(A) gives the inverse of a square matrix A diag(A) gives a column vector with elements of the main diagonal of A trace(A) gives the sum of diagonal elements of a square matrix A expm(A) = I + A +A²/2! + A³/3!+.....(Matrix exponential function, I denotes the unit matrix having diagonal elements of value 1, all other elements are zero)

For the coefficients of an antisymmetric matrix holds aij =-aji. The definition implies that the diagonal elements are all zero. Show that the cross product of two vectors c = a × b can be obtained via an antisymmetric ⎡ 0 − a z a y ⎤ ⎡b x ⎤ ⎡c x ⎤ ⎢ ⎥ ⎢ ⎥ matrix A times b: c = A ⋅ b = ⎢c y ⎥ = ⎢ a z − a x ⎥ ⋅ ⎢⎢b y ⎥⎥ 0 ⎢− a y a x ⎢⎣ c z ⎥⎦ 0 ⎥⎦ ⎢⎣bz ⎥⎦ ⎣ Using antisymmetric matrices it is possible to perform rotation operations: Try to solve the follwing problem: ⎡1⎤ Given a vector ⎢⎢1⎥⎥ . The vector is rigidly rotated about the z-axis by 90° ⎢⎣1⎥⎦

⎡− 1⎤ to obtain the vector ⎢⎢ 1 ⎥⎥ . Determine a matrix D such that ⎢⎣ 1 ⎥⎦

⎡1⎤ ⎡− 1⎤ ⎢ 1 ⎥ = D ⋅ ⎢1⎥ ⎢⎥ ⎢ ⎥ ⎢⎣1⎥⎦ ⎢⎣ 1 ⎥⎦

Solution: ⎡0 ⎤ The unit vector along the rotation axis is ⎢⎢0⎥⎥ . The corresponding antisymmetric matrix is ⎢⎣1 ⎥⎦ ⎡0 − 1 0 ⎤ given by ⎢⎢1 0 0⎥⎥ . The matrix performing rotation by 90° (or π/2) about the z-axis is then ⎢⎣0 0 0⎥⎦ Breitbach, Math. Tools and Simulation, SS2008

11

⎛⎡ 0 − π / 2 0⎤ ⎞ ⎡0 − 1 0⎤ ⎜⎢ ⎟ ⎥ exp m⎜ ⎢π / 2 0 0⎥ ⎟ = ⎢⎢1 0 0⎥⎥ ⎜⎢ 0 0 0⎦⎥ ⎟⎠ ⎣⎢0 0 1⎦⎥ ⎝⎣

Instructions&Problems 5.3:

A linear equation system can be written in the following form: ⎡ a11 . . . a1n ⎤ ⎡ x1 ⎤ ⎡ b1 ⎤ ⎢ . . . . . ⎥⎥ ⎢⎢ . ⎥⎥ ⎢⎢ . ⎥⎥ ⎢ = ⋅ ⎢ . . . . . ⎥ ⎢ . ⎥ ⎢ . ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎣a m1 . . . . a mn ⎦ ⎣ x m ⎦ ⎣bm ⎦ Solve the linear equation system:

⎡1 ⎢1 ⎢ ⎢6 ⎢ ⎣1

4 0 1 2

2 3 5 1

5⎤ ⎡ x1 ⎤ ⎡ 3⎤ ⎥ ⎢ ⎢1 ⎥ ⎥ 1⎥ ⎢ x 2 ⎥ = ⎢ ⎥ ⋅ ⎢ 2⎥ 0⎥ ⎢ x3 ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ 2⎦ ⎣ x 4 ⎦ ⎣ 4⎦

Given the linear equation system shown above, create an other one having the same solution vector x.

Solve the linear equation system: ⎡14 ⎤ ⎡ 1 2 4 1⎤ ⎡ x1 ⎤ ⎢ ⎥ ⎢14 13 14 5⎥ ⎢ x ⎥ ⎥ ⋅ ⎢ 2 ⎥ = ⎢ 73⎥ ⎢ ⎢44⎥ ⎢ 9 8 8 3⎥ ⎢ x 3 ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎣30⎦ ⎣ 8 6 4 2⎦ ⎣ x 4 ⎦

What is the rank of the matrix? Is ist possible to compute the inverse matrix? Create the augmented (4x5) matrix C with the right hand side as the fifth column, then use the command ROE = rref( C ) to create the matrix ROE. If ROE contains rows with all zeros, the system contains redundant information If ROE contains a row where all elements but the last are zero, the system contains a contradiction, no solution exists. A unique solution is found in the last column of ROE. Breitbach, Math. Tools and Simulation, SS2008

12

Symbolic methods to solve linear equations. Try the following command [x1,x2] = solve( ‘ a*x1 + b*x2 = c, d*x1 + e*x2 = f ‘ , ‘ x1, x2‘)

to obtain the solution for x1, x2 in general terms. Instructions&Problems 5.4:

The command eig(A) calculates the eigenvalues and the eigenvectors of a matrix A. Typing [V, L]=eig(A) generates a full matrix V containing in the columns the eigenvectors. L is a diagonal matrix with the eigenvalues as components. The task connected with eigenvectors and eigenvalues is to find vectors and scalar values for which the following relation holds: ⎡ a11 ⎢ . ⎢ ⎢ . ⎢ ⎣a n1

. . . .

. a1n ⎤ ⎡ x1 ⎤ ⎡ x1 ⎤ ⎢.⎥ ⎥ ⎢ ⎥ . . ⎥ ⎢.⎥ = λ ⋅⎢ ⎥ ⋅ ⎢.⎥ . . ⎥ ⎢.⎥ ⎢ ⎥ ⎥ ⎢ ⎥ . a nn ⎦ ⎣ x n ⎦ ⎣ xn ⎦

6.Script Files Instructions&Problems 6.1:

What is a *.m file or script file? Type at the prompt>> welcome. Then at the prompt >>edit welcome.m It is possible to assemble commands in a file filename.m. By calling filename the commands are executed one after the other. Add the following line in the EDITOR/DEBUGGER to the welcome.m-file: disp(´I will do my best in this MATLAB course`); (Hint: An expression enclosed by apostrophes is a string.) That is the principal way to write programs in MATLAB. Organize the statements in *.m files. They are not restricted to such simple things as we did in our welcome.m-file. Insert edit. In the EDITOR/DEBUGGER environment create the following file (close every statement with a semicolon!): a=13; b=5; c=7; V=a*b*c; S=2*(a*b + b*c + a*c); Save the file as geocalc.m (or so). It is a good programming practice to reflect carefully the names of our programs! Then leave the editor and type geocalc to obtain values of V and S. Of course this program is related to a very very special situation. But before we generalize the program we will lern about the documentation of a program. The help function is a very mighty instrument also in connection with private programs. Edit geocalc.m. Than type the following statements in line one and two: % The program geocalc computes the volum and the surface of a brick element Breitbach, Math. Tools and Simulation, SS2008

13

%with sides a=13, b=5, and c=7. Save and at the MATLAB prompt>> call help geocalc Lines starting with „%“ are interpreted as non executable comment lines. The first lines of a private script file can be shown via the „help“. In case that you have written a useful program you should not foreget to add comments as clear as possible. How to generalize the programm? Erase the three lines specifying a,b,c.Save as geocalcnew and insert at the prompt: a=3; b=5; c=7; geocalcnew.

Instructions&Problems 6.2:

The program is now more flexible to use. But it is possible to introduce a user defined function: function [outarg1, outarg2,....] = functionname (inarg1, inarg2,......) % Comments ........... executable code ........... return The function is called functionname(inarg1, inarg2,....) Organize geocalc as a function geocalcfun(a,b,c). Observe that the variables used in the statements of the function are l o c a l. Add the following statement in function geocalcfun and create the enhanced function geocalcfun_1: a = 2*a, b = 2*b, c = 2*c. Save, quit the EDITOR and specify a=5, b=8, c=9, then call [Vol,Surf] = geocalcfun_1(a,b,c). Check the values of a, b, c.

Instructions&Problems 6.3:

A “while loop” is an arrangement of statements that is repeated as long as some condition is satisfied. The general form is while expression ......... executable code ......... end Type in: i=0;while (i<10); i=i+1; disp(i);end Organise the 5 statements in a script-file named exwhile.m. A “for loop” is an arrangement of statements that is repeated as long as some condition is satisfied. The “for loop” has the form For index = expression (as a rule we have as expression: first:increment:last) Breitbach, Math. Tools and Simulation, SS2008

14

........... executable code ........... end Calculate the sum 3+6+9+12+15 +...+3*j + .....303: sum = 0; for i=1:101; sum=sum+3*i;end or sum = 0; for i=0:3:303; sum=sum+i;end

Instructions&Problems 6.4:

The if construct has the following form if expression_1 ......... executable code ......... elseif expression_2 ......... executable code .........

elseif expression_i ...... executable code .........

else ......... executable code ......... end Organize the following program via a function script file: Given the price structure of some articles depending on the bought amount. Number of articles 0 - 10 11 - 30 31 - 50 > 50

Price per article 40 Euro 37 Euro 34 Euro 31 Euro

function price_fun(number) Breitbach, Math. Tools and Simulation, SS2008

15

%Calculation of total prices if number>50 pt = number*31; fprintf('total price is %4.0f elseif number>30 pt = number*34; fprintf('total price is %4.0f elseif number>10 pt = number*37; fprintf('total price is %4.0f else pt = number*40; fprintf('total price is %4.0f end

times 31 Euro = %7.0f Euro\n',number,pt);

times 34 Euro = %7.0f Euro\n',number,pt);

times 37 Euro = %7.0f Euro\n',number,pt);

times 40 Euro = %7.0f Euro\n',number,pt);

Remarks to format output: %d → integer %e → exp-format %f → floating point format \n → skip to a new line What do you think about the following variant of the program? function price_fun(number) %Calculation of total prices if number>0 pt = number*40; fprintf('total price is %4.0f elseif number>10 pt = number*37; fprintf('total price is %4.0f elseif number>30 pt = number*34; fprintf('total price is %4.0f else pt = number*31; fprintf('total price is %4.0f end

times 40 Euro = %7.0f Euro\n',number,pt);

times 37 Euro = %7.0f Euro\n',number,pt);

times 34 Euro = %7.0f Euro\n',number,pt);

times 31 Euro = %7.0f Euro\n',number,pt);

Instructions&Problems 6.4: Other branching realisations can be performed via switch/case constructions:

Switch (switchexpression) case (caseexpression_1), ......... executable code ......... case (caseexpression_2), ......... executable code ......... case (caseexpression_i), ......... executable code Breitbach, Math. Tools and Simulation, SS2008

16

......... otherwise, ......... executable code ......... end The case expressions of the case code blocks are compared successively with the switch expression. If caseexpression_i is equal to switchexpression the corresponding code block is executed. Then there is a jump to the end statement of the switch/case construct. The „otherwise“ block is optional. The expressions may be either numerical or string variables. In the following for a piecewise defined function f(x) the switch construct is applied.

DOMAIN 0 <= x <1 1 <= x <2 2 <= x <3 x >= 3

RANGE f(x) = x f(x) = x² f(x) = x + 2 f(x) = 5

function y = func(x) for k=1:length(x) t=x(k); boolexp=1; switch boolexp case ((t>=0)&(t<1)) y(k)=t; case ((t>=1)&(t<2)) y(k)=t*t; case ((t>=2)&(t<3)) y(k) = 2 + t; case(t>=3) y(k)=5.0; end end

7. Differential Equations

Instructions&Problems 7.1: 1

Solve numerically the differential equation y ′ = f ( x, y ) with f ( x, y ) = x ⋅ y 3 The simplest and most straightforward procedure is the Euler-Cauchy method, the algorithms are given in the following:

Breitbach, Math. Tools and Simulation, SS2008

17

y 2 = y ( x2 ) ≈ y1 + h ⋅ tan(ϕ ) = y1 + h ⋅ y ′( x1 ) = y1 + h ⋅ f ( x1 , y1 )

Specify h; i=1 ↓ xi = xa; yi=ya (Initial conditions) ↓ i ← i+1 ↓ xi ← xi-1 + h yi ← yi-1 + h ⋅ f (xi-1, yi-1) no ↓ xi > xb ? ↓ yes stop

%Euler-Cauchy-Method %Numerical solution of the DEQ y' = x*y^(1/3) %Solution for initial conditions y(1)=1 is: y(x) = ((x²+2)/3)^(3/2) Breitbach, Math. Tools and Simulation, SS2008

18

xa=1; xb=5; ya=1;i=1; h=0.1; x(1)=xa;y(1)=ya; while (x<xb) i=i+1; x(i) = x(i-1) + h; y(i) = y(i-1) + h*x(i-1)*y(i-1)^(1/3); end yex = ((x.^2+2)./3).^(3./2); plot(x,y,x,yex);grid on;

Instructions&Problems 7.2: MATLAB provides powerful routines to solve differential equations. The user has to prepare a special function as script file which is linked to the DEQ solver routine.

Differental equations must be organized as systems of first order differential equations. Let’s assume that we have the two differential equations: y1′′′ + a ⋅ y 2′ + b ⋅ y1 = f1 (t ) y 2′′ + c ⋅ y1′′ − d ⋅ y 2 = f 2 (t ) Define now: z1 = y1 , z 2 = y1′ , z 3 = y1′′ , z 4 = y 2 , z 5 = y 2′ And the two equations can be rewritten in the following form: z1′ = z 2 z ′2 = z 3

z 3′ = f 1 (t ) − a ⋅ z 5 − b ⋅ z1

or

⎡ z1′ ⎤ ⎡ 0 ⎢z′ ⎥ ⎢ 0 ⎢ 2⎥ ⎢ ⎢ z 3′ ⎥ = ⎢− b ⎢ ⎥ ⎢ ⎢ z ′4 ⎥ ⎢ 0 ⎢⎣ z 5′ ⎥⎦ ⎢⎣ 0

0 ⎤ ⎡ z1 ⎤ ⎡ 0 ⎤ 0 ⎥⎥ ⎢⎢ z 2 ⎥⎥ ⎢⎢ 0 ⎥⎥ 0 0 0 − a ⎥ ⋅ ⎢ z 3 ⎥ + ⎢ f1 (t )⎥ ⎥ ⎥ ⎢ ⎥ ⎢ 0 0 0 1 ⎥ ⎢z4 ⎥ ⎢ 0 ⎥ 0 − c d 0 ⎥⎦ ⎢⎣ z 5 ⎥⎦ ⎢⎣ f 2 (t )⎥⎦

1 0

0 1

0 0

z 4′ = z 5 z 5′ = f 2 (t ) − c ⋅ z 3 + d ⋅ z 4 Note that the derivatives of the zi(t) are written explicitly as functions of the zi(t). The following example is worked out to demonstrate the solution procedure. It is the differential equation of a pendulum whose motion is not restricted to small angles:

ϕ&& + ( g / l ) ⋅ sin(ϕ ) = 0

Breitbach, Math. Tools and Simulation, SS2008

19

function [phidot] = pendulum(t,phi) %Differential equation of a simple pendulum % Call [t,phi]=ode23('pendulum',[0 tend],[phi0,omega0]) l=2.5; g=9.81;

%length (m) %gravity acceleration

%Initialisation of the vector of derivatives phidot=[0;0]; %Differential equation system of first order phidot(1) = phi(2); phidot(2) = -(g/l)*sin(phi(1));

Instructions&Problems 7.3:

Evaluate the oscillation periods of a pendulum for large amplitudes. Write a single script file for the calculation and the plotting of the period as a function of the amplitude. %Study the oscillation periods of a pendulum in case of %larger amplitudes %use pendulum.m and ode-routine %pendulum parameters clear; l = 2.5; %m length g = 9.81; %m/s² acceleration of gravity Tchar=2*pi*sqrt(l/g); for ii=1:50 amp0(ii)=ii*pi/75; [t,phi]=ode23('pendulum',[0,10],[amp0(ii);0]); ki=1; while (t(ki)< 0.8*Tchar) ki=ki+1; end kj=ki; while (phi(kj)
Breitbach, Math. Tools and Simulation, SS2008

20

Instructions&Problems 7.4:

A simple introduction example for the using of the SIMULINK program. Integrate the piecewise defined linear function

1 Repeating Sequence

s Integrator

Scope output

Scope input

Solve the pendulum problem with the help of SIMULIMK:

A simulink modeling of the pendulum problem.

Breitbach, Math. Tools and Simulation, SS2008

21

Instructions&Problems 7.5:

SIMULINK has capabilities to create subsystems. It is demonstrated via a very simple model for the heat load of a house. The temperature behaviour of the air inside is influenced by a heater, by losses via transmissions through the walls and by air exchange (infiltration rate).

The differential equation for the system is as follows: M ⋅ c p ⋅ T&in = Q& Heater + m& air ⋅ c p ⋅ (Tout − Tref ) − m& air ⋅ c p ⋅ (Tin − Tref ) −

1 ⋅ (Tin − Tout ) Req

or

M ⋅ c p ⋅ T&in = Q& Heater + ( m& air ⋅ c p +

1 ) ⋅ (Tout − Tin ) Req

M : Air mass inside the house [kg] c p : Specific heat of air (constant pressure) [J/kg/°C]

Tin : Temperature inside the house [°C] Tout : Outside temperature [°C] Tref : A reference temperature (e.g. 0 °C) [°C] : Heating rate [W] Q& heater

m& air : Infiltration rate (air exchange rate) [kg/s] Req : Thermal resistance of the house boundaries [°C/W]

The construction of the SIMULINK model is straight foreward.

Breitbach, Math. Tools and Simulation, SS2008

22

Selecting the SIMULINK-objects inside the box ⇒ Clicking “edit” then “Create subsystem” then the following figure is obtained:

The subsystem box is the container for the connected objects enclosed by the rectangle. Only the input and output gates of the subsystem are visible. Clicking the subsystem the internal structure is shown. Out1 is here the indoor temperature. Instructions&Problems 7.6:

How to solve the equation x = sin( x ) + cos( x) − 0.5 with the help of SIMULINK? This question leads to so called “algebraic loops” or “algebraic constraints”. Building the model as shown below will not be successful, because solution procedures for this type of problems in SIMULINK are based on iterations.

Breitbach, Math. Tools and Simulation, SS2008

23

⇒ Will probably not work!!!!

Appropriate model using the “Algebraic Constraint”- object

Example of using “Algebraic Constraint”: Electrical network with diode

It is characterized by the following set of three equations for the quantities u a (t ) , u d (t ), id (t ) : (2 differential equations and one nonlinear equation)

du a (t ) u (t ) = − a + id (t ) dt R u g (t ) = R g ⋅ id (t ) + u d (t ) + u a (t )

C⋅

id (t ) = I 0 ⋅ (exp(u d (t ) / VT ) − 1) SIMULINK-model: Breitbach, Math. Tools and Simulation, SS2008

24

Results:

%detini.m: Initialisation Rg=10e3; R=100e3; C=4000e-12; I0=20e-9; VT=25e-03; fs=3000;

Breitbach, Math. Tools and Simulation, SS2008

25

Simulation Example A A Simple Dynamic Absorber In 1911 H. Frahm has invented a device for stabilisation of rocking oscillations of ships. This device is now called a dynamic absorber. Based on a simple idea it is a great invention in vibration control, frequently used in applications.

Let´s start with a single degree of freedom system, excited by a harmonic force f (t ) = F0 ⋅ sin (ω ⋅ t ) .

The Primary System

This is the primary system. The primary mass vibrates with two frequencies, the natural kp and the forced frequency of vibration ω . frequency of vibration ω n = mp Our objective is to eliminate the component of the forced vibration of the primary mass mp. This is achieved by attaching a secondary single degree-of-freedom system, with mass ms and spring ks, to the primary system as shown below.

The Controlled System The problem is to choose the right parameters, ms and ks, of the secondary system. The equations of motion for the controlled system are

⎡m p ⎢0 ⎣

0 ⎤ ⎡ &x&p ⎤ ⎡k + k s ⋅⎢ ⎥ + ⎢ p ⎥ m s ⎦ ⎣ &x&s ⎦ ⎣ − ks

− ks ⎤ ⎡x p ⎤ ⎡F ⎤ ⋅ ⎢ ⎥ = ⎢ 0 ⎥ ⋅ sin (ω ⋅ t ) ⎥ k s ⎦ ⎣ xs ⎦ ⎣0⎦

Breitbach, Math. Tools and Simulation, SS2008

26

Simulink Model of the controlled system SIMULINK-MODEL of a dynamic vibration absorber (forced vibrations of a two degree of freedom system)

Ramp

1

xpdot

s Integrator

sin

10

sin

Gain

Scope force

xp 1

-14

s Integrator1

(kp + ks)/M

Sum

10 Scope xp

ks/m -10 -ks/m

1 Sum1

s Integrator2

xsdot

1

xs

s Integrator3

2 ks/M

Scope xs

Breitbach, Math. Tools and Simulation, SS2008

27

Simulation Example B The Motion of a Structure Consisting of Two Hinged Bars Carrying Two Masses m, M The bar L = a+b has a pivot point which is taken as the origin of a cartesian coordinate system

Coordinates of the particles (m,M): x1 = − a ⋅ sin (ϕ )

x 2 = b ⋅ sin (ϕ ) + l ⋅ sin (θ )

y1 = − a ⋅ cos (ϕ )

y 2 = b ⋅ cos(ϕ ) + l ⋅ cos(θ )

Velocities obtained by differentiation:

x&1 = − a ⋅ cos(ϕ ) ⋅ ϕ&

x&2 = b ⋅ cos(ϕ )⋅ϕ& + l ⋅ cos(θ )⋅θ&

y&1 = a ⋅ sin (ϕ ) ⋅ ϕ&

y& 2 = − b ⋅ sin (ϕ )⋅ϕ& − l ⋅ sin (θ ) ⋅θ&

Kinetic Energy: Mass m:

Mass M: Ekin , M =

(

)

(

)

1 1 2 2 ⋅ m ⋅ x&1 + y&1 = ⋅ m ⋅ a 2 ⋅ϕ& 2 2 2

Ekin , m =

(

1 1 2 2 ⋅ M ⋅ x&2 + y& 2 = ⋅ M ⋅ b 2 ⋅ ϕ& 2 + l 2 ⋅ θ& 2 + 2 ⋅ b ⋅ l ⋅ cos(ϕ − θ ) ⋅ ϕ& ⋅ θ& 2 2

)

Total Kinetic Energy: T = Ekin ,m + E kin,M =

m ⋅ a2 + M ⋅ b2 2 M ⋅ l 2 &2 ⋅ϕ& + ⋅θ + M ⋅ b ⋅ l ⋅ cos(ϕ − θ ) ⋅ ϕ& ⋅θ& 2 2

Breitbach, Math. Tools and Simulation, SS2008

28

Potential Energy:

V = − (m ⋅ g ⋅ y1+ M ⋅ g ⋅ y 2 ) = m ⋅ g ⋅ a ⋅ cos(ϕ ) − M ⋅ g ⋅ (b ⋅ cos(ϕ ) + l ⋅ cos(θ )) Lagrangian Function ( L = T –V ):

L = T

− V

m ⋅ a2 + M ⋅ b2 2 M ⋅ l 2 &2 = ⋅ϕ& + ⋅ θ + M ⋅ b ⋅ l ⋅ cos(ϕ − θ ) ⋅ ϕ& ⋅ θ& 2 2 + (M ⋅ b − m ⋅ a ) ⋅ g ⋅ cos(ϕ ) + M ⋅ l ⋅ g ⋅ cos(θ )

Equations of Motion derived from the Lagrangian Function : ⎛ ∂L ⎞ d ⎛ ∂L ⎞ ⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ = 0 dt ⎝ ∂q& ⎠ ⎝ ∂q ⎠

for the var iables q ≡ ϕ and q ≡ θ

Two Differential Equations for ϕ(t) and θ(t): Defining terms A, B, C, RS, RT:

A = m ⋅ a2 + M ⋅ b2

B = M ⋅ b ⋅ l ⋅ cos(ϕ −θ )

C = M ⋅l2

RS = − M ⋅ b ⋅ l ⋅ sin (ϕ −θ ) ⋅ θ& 2 − (M ⋅ b − m ⋅ a ) ⋅ g ⋅ sin (ϕ ) RT = M ⋅ b ⋅ l ⋅ sin (ϕ −θ ) ⋅ ϕ& 2 − M ⋅ l ⋅ g ⋅ sin (θ ) the following equations are obtained:

ϕ&& =

C ⋅ RS − B ⋅ RT A⋅C − B2

θ&& =

A ⋅ RT − B ⋅ RS A⋅C − B2

Files related to this example

The following function program is to use with the ode-routines: function [zdot]=trebuch1(t,z) %trebuch1 m = 0.25;

%kg; mass to move

Breitbach, Math. Tools and Simulation, SS2008

29

M L a b l g

= = = = = =

6; 2.5; 1.75; L-a; 0.43; 9.81;

%kg; "big" mass %m ; length of the main bar %m ; lever length to small mass %m ; lever length to big mass %m ; length of the pendulum bar %m/s²; gravity acceleration

zdot=[0;0;0;0]; %z(1) =

phi; z(2) = dphi/dt; z(3) = theta; z(4) = dtheta/dt

zdot(1) = z(2); %Precalculate the long terms A=m*a^2+M*b^2; B=M*b*l*cos(z(1)-z(3)); C=M*l^2; RS=-(M*b*l*sin(z(1)-z(3)))*z(4)^2 - (M*b-m*a)*g*sin(z(1)); RT=M*b*l*sin(z(1)-z(3))*z(2)^2 - M*l*g*sin(z(3)); zdot(2) = (C*RS - B*RT)/(A*C-B^2); zdot(3) = z(4); zdot(4) = (A*RT - B*RS)/(A*C-B^2);

The following script file evaluates the kinetic energy of mass M. We start with the formulation in cartesian co-ordinates. The velocity components are expressed via the angles variables. Algebraic transformations yield the term Ekin,M.

Script file to calculate and simplify the kinetic energy: %Calc_EkinM connected with the trebuchet problem syms M b l xdot ydot phi theta phidot thetadot Ekin; xdot=b*cos(phi)*phidot + l*cos(theta)*thetadot; ydot=-b*sin(phi)*phidot - l*sin(theta)*thetadot; Ekin=1/2*M*(xdot^2 + ydot^2); simple(Ekin)

Breitbach, Math. Tools and Simulation, SS2008

30

Simulation Example C Energy Balance of a House The energy balance of a house is determined by transmission heat losses and losses by air exchange. In Instructions&Problem 7.5 a simple model of a house is discussed. Embed this model into an enlarged model with the following features:

The outdoor temperature varies sinuslike about a mean value. The heating system is controlled via a bang bang servo (Thermostat). Realise it in SIMULINK via Relay under “Discontinuities”

Use the following initialisation file for the parameter of the system: (Use it as a “model pre-load function”!) % HOUSE_DATA % This script runs in conjunction with the % househeating-simulink-file. % ------------------------------% Define the house geometry % ------------------------------% House length = 12 m lenHouse = 12; % House width = 10 m widHouse = 10; % House height = 2.75 m htHouse = 2.75; % Roof pitch = 45 deg pitRoof = 30/(180/pi); % Number of windows = 6 numWindows = 8; % Height of windows = 1.5 m htWindows = 1.5; % Width of windows = 1.5 m widWindows = 2.0; windowArea = numWindows*htWindows*widWindows; Breitbach, Math. Tools and Simulation, SS2008

31

wallArea = 2*lenHouse*htHouse + 2*widHouse*htHouse + ... 1/cos(pitRoof)*widHouse*lenHouse + ... tan(pitRoof)*widHouse*widHouse - windowArea; % ------------------------------% Define the k-values (W/m²K) % ------------------------------% Wall kWall = 0.4; RWall = 1/(kWall*wallArea); % Glass windows kWindow = 2.5; RWindow = 1/(kWindow*windowArea); % ------------------------------% Determine the equivalent thermal resistance for the whole building % ------------------------------Req = RWall*RWindow/(RWall + RWindow); % ------------------------------% c = cp of air (273 K) = 1005.4 J/kg-K % ------------------------------c = 1005.4; % ------------------------------% Enter the heating rate QHeat (W) % ------------------------------QHeat = 12000; % ------------------------------% Air change flow rate Mdot (kg/s) % ------------------------------Mdot = 0.1; % ------------------------------% Determine total internal air mass = M % ------------------------------% Density of air at sea level = 1.2250 kg/m^3 densAir = 1.2250; M = (lenHouse*widHouse*htHouse+tan(pitRoof)*widHouse^2*lenHouse/4)*densAir; %Initial indoor temperature TinIC=20.0;

Breitbach, Math. Tools and Simulation, SS2008

32