MATLAB - EDB/ 第1 章 行列の入力と加工 MATLAB は行列の演算を基本とするソフトです.そこでまずは行列を作ることから始めましょう. その後で行列の要素をみたり,行列を加工する練習をしましょう. 1.1 ベクトルと行列要素の作り方 ここでは ...

  • Published on
    15-May-2018

  • View
    220

  • Download
    4

Transcript

i I MATLAB iii 1 11.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2 162.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3 263.1 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 II MATLAB 37 4 334.1 MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.5 . . . . . . . . . . . . . . . . . . . . . . . . . . 424.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 ii4.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5 GUI 505.1 . . . . . . . . . . . . . . . . . . . . . . . . . . 505.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.3 GUI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 625.4 GUI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.5 MATLAB . . . . . . . . . . . . . . . . . . . 725.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 III MATLAB 74 6 756.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 766.3 Phase Portrait DemoPL . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786.4 DufFix . . . . . . . . . . . . . . . . . . . . . . . . . . . . 816.5 DemoBF . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 7 867.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 867.2 Poincare . . . . . . . . . . . . . . . . . . . . . . . . . . . 877.3 Rossler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 8 Switched Dynamical Systems with Moving Border 968.1 Formation of the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 968.2 The Poincare map and periodic solution . . . . . . . . . . . . . . . . . . . . . . . 998.3 Phase Portrait DemoPL . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1038.4 DufFix . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1038.5 DemoBF . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 A 104 B 105B.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105B.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106B.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 IMATLAB 1 1MATLAB 1.1 1.1.1 [ ; ] 1.1 1. A=0:102. t=0 : pi/2: 2*pi3. B=[1 2 3; 4 5 6]4. A=[B; 7 8 9]5. H=Hello world!x = xmin : x : xmax[xmin, xmin + x, xmin + 2x, , xmax] x 1 t=0 : pi/30 : 2*pi; plot(t, sin(t)) 1 2[] 1.1 3. 2 3 2 3 B =1 2 34 5 6 B 4. B 3 [7 8 9] 3 3 A A =1 2 34 5 67 8 9 A=1:100; 1.1 5. 1.1.2 [ ][] 1 (array constractor) [ ] (concatenation) 1.1 1. A = [5] , B =2 1 11 5 11 1 2, C =1 + 2j23 5j , D =1 5j2 + 4j3 3j4 + 2j5 j2. A =2 4 6 8 1010 8 6 4 2 , B =1 3 5 7 99 8 7 6 53. A =[Solving our problems with MATLAB.], B = [I like A] 1. A = [5] A = 5 1 1 A = [5] B = 5 2 A==B 1 31 A = B 0 3. j 1 i j i 2 sqrt(2) 4. A = [2 : 2 : 10 ; 10 : 2 : 2]A = [(2 : 2 : 10) ; (10 : 2 : 2)] [A B] [A ; B]5. A B I like Solving our problems with MATLAB.1.2 1.2.1 MATLAB (array) *1 1.1 2 MATLAB NaN Not a Number 0/0 (matrix)*2 (scalar)*3 m n m n m n *1 C Pascal Fortran MATLAB *2 matrices *3 i j 1 4 1.1 MATLAB a = 5.5, b = 2 pi+ 3j1 1 pi i j NaN 1 n a = [1 2 3 4 5] m 1 a =123 m n A =1 2 35 3 1 0 0 A = [ ]zeros(0,4) 0 4 zeros(3,0) 3 0 eye(3) ones(3,4) 1 zeros(3,4) 0 rand(3,4) magic(4) 4 4 pascal(4) Pascal 1 5A =a11 a12 a1na21 a22 a2n.........am1 am2 amn= [apq] (1.1)apq (element) (p, q) *42 Avec = vec (A) =a11...am1a12...am2...a1n...amn=12..................r1r(1.2) A 1 1 (array) A 1 apq = r (p, q) r r = (q 1)m+ p (1.3)MATLAB apq r ( )apq A (p, q)r A (r)2 1 (1.3) *4 aij i j MATLAB p q 1 61.2.2 1 M M=magic(5)M =17 24 1 8 1523 5 7 14 164 6 13 20 2210 12 19 21 311 18 25 2 9M = [mpq] , p, q = 1, 2, , 5 mpq 2 3 M( 2, 3 )7 M( 12 ) (1.3) 12 = (3 1)5 + 2 M( 12 : 15 ) M( [12 13 14 15] ) p 2 M( 2, : ) M( end, : ) 2 4 M( 2 : 4, : ) 1 7 M( 4 : 2, : )Empty matrix: 0-by-5 0 5 q 3 M( : , 3 ) M( : , end ) 3 4 M( : , 3 : 4 )2 3 4 M( 2 , 3 : 4 ) 2 4 3 4 N N=M( 2 : 4 , 3 : 4 )1.3 1 1 (1.1) (1.2) M( : ) size( ans ) ans= 25 12 25 1 2 M=magic(5) M( 2, : )=[ ] 1 8 2 3 M=magic(5) M( 2 : 3, : )=[ ] 3 M=magic(5) M( :, 3 )=[ ] 2 3 M=magic(5) M( :, 2 : 3 )=[ ] 2 2 20 M=magic(5) M( 2 : 2 : 20 )=[ ] X X=1 : 10 X( 10 : -1 : 1 ) M=magic(5) M( [1,5],: )=M( [5,1],: ) M=magic(5) M( [2, 5,1],: ) M=magic(5) N=[ 1 2 3 4 5 ] M(2, : )=N M=magic(5) N=[ 1 2 3 4 5 ] M(2, : )=N 1 9 1.2 MATLAB MATLAB diag rot90 90 tril triu fliplr flipud N M=magic(5) N=[ 1 2; 3 4] M(2:3, 3:4 )=N N M=magic(5) N=[ 1 2; 3 4] M([1, 3], [3, 4] )=N[ ] 1.2 1.4 ` A1, A2, , A`MATLAB 3 3 3 4 A=ones(3) 1 10 A( : , : , 2)=2*ones(3) A( : , : , 3)=3*ones(3) A( : , : , 4)=4*ones(3) size( A ) ans= 3 3 42 A 4 MATLAB (dimension) 3 2 A( : , : , 2) 1 1 2 1 [A( : , 1 , 1), A( : , 1 , 2)]3 1.3 whos whos(A) Name Size Bytes Class A 3 3 4 228 double array Grand total is 36 elements using 228 bytes1.5 clear wild card * clear AA 1 11 1.3 size size(A)ans= 2 3 4 2 3 4 ndims ndims(A)ans= 3 whos whosName Size Bytes Class 1.6 1.6.1 1. MATLAB 2. [ ]3.4. Hello 1. 2 A(2, 3) 2 3 1 A(12) (1.3) 2. 3. 4. size 1 12arraychar numericdoublesparseunit8cell struct 1.1 1.6.2 MATLAB Data Types There are six fundamental data types(classes) in MATLAB, each one a multidimensional array.The six classes are double, char, sparse, init8, cell, and struct. The two-dimensional versions of thesearrays are called matrices and are where MATLAB gets its name.You will probably spend most of your time working with only two of these data types: the doubleprecision matrix (double) and the character array (char) or string. This is because all computationsare done in double-precision numbers or strings.The other data types are for specialized situations like image processing (unit8), sparse matrices(sparse), and large scale programming (cell and struct).You cant create variables with the types numeric or array. These virtual types serve onlyto group together types that share some common attributes.The unit8 data type is for memory efficient storage only. You can apply basic operations suchas subscripting and reshaping to these types of arrays but you cant perform any math with them.You must convert such arrays to double via the double function before doing any math operations.You can overload methods for the build-in data types in exactly the same way you overload amethod for an object. For example to define sort for a unit8 array, create a method(sort.m) andplace it into @unit8 within a directory on your path.The table that follows describes the data types in more detail. (struct) C C 1 13 1.4 MATLAB Class Example Descriptiondouble [1 2; 3 4] Double precision numeric array (this isthe most common MATLAB variable type).char Helloe Character array (each character is 16 bitslong). Also referred to as string.sparse speye(5) Sparse double precision matrix (2-D only).The sparse matrix stores matrices withonly few non-zero elements in a fractionof the space required for an equivalent fullmatrix. Sparse matrices invoke specialmethods especially tailored to solve sparseproblems.cell {17 hello eye(3)} Cell array. Elements of cell array containother arrays. Cell arrays collect relateddata and information of a dissimilar sizetogether.struct a.date=12; Structure array. Structure arrays havea.mat=magic(5); field names. The fields contain othera.color=red; arrays. Like cell arrays, structures collectrelated data and information together.unit8 unit8(magic(3)) Unsigned 8 bit integer array. The unit8array stores integers in the range from 0to 255 in 1/8 the memory required for adouble precision array. No mathematicaloperations are defined for unit8 arrays.UserObject inline(sin(x)) User-defined type. 1 14 [ ] [ ] (nest) [ ] { } C={[1, 2], [1 2 3;4 5 6], [7 8; 9 0; 2 1]} C 3 C{2} 1 151.7 1. [[[]]] [],[],[] [[],[]] [[]],[] {{{}}} {},{},{} {{},{}} {{}},{}2. eye(4) 4 4 0 1 0 00 0 1 00 0 0 11 0 0 0,0 0 1 00 0 0 11 0 0 00 1 0 0,0 0 0 11 0 0 00 1 0 00 0 1 0,1 0 0 00 0 1 00 1 0 00 0 0 1 shiftrotate(A) 1 3. 4 4 0 1 0 00 0 1 00 0 0 10 0 0 0,0 1 0 00 0 1 00 0 0 00 0 0 0,0 1 0 00 0 0 00 0 0 10 0 0 0,0 0 0 00 0 0 00 0 0 00 0 0 04. MATLAB Help menu Examples and Demos Demo Window Matrixmanipulation Demo Windows 95 1 3 Demo 5. MATLAB Help Window matlab\elmat compan, hadamard, hilb, vander magic, pascal Hadamard 2nhadamard(2), hadamard(4), hadamard(8) 16 2MATLAB 2.1 1. 2 2. 3. 3 MATLAB MATLAB MATLAB 2 17 2.1 a = [ar] , b =[br] A = [apq] , B = [bpq] 2.2 2.1 1. t=0:10 t*5 t.2 t /2 t2. A=[1 2 3; 4 5 6; 7 8 9] A*2 A.2 A/2 A3. A=[1 2 3; 4 5 6; 7 8 9] B=2*ones(3) A.*B A*B A./B 2*A-3*B4. A=[1 1;1 -1] kron(A, A) 2.1 1.. x=-5 : 0.1 : 5f (x) = x3 x2 5x 2 y = x.3 x.2 5 x 2 2 18 2.1 MATLAB a = [ar] a or a 1a =[1ar]a/ or \aA = [apq] A or A 1A =[1apq]A/ or \A[(apq)n] A.n a+ b = [ar + br] a+ b a b = [ar br] a b A+B = [apq + bpq] A+B AB = [apq bpq] AB AB = [apqbpq] A. B AB = [apq/bpq] A./B AB =[kapkbkq]A B Kronecker AB = [apqB] kron(A, B) AX = B X = A1B A\B XA = B X = BA1 B/A An An At A. At A 2 19 plot(x, y) grid4. Kronecker 1 11 11 11 1 =11 11 1 11 11 111 11 1 11 11 1=1 1 1 11 1 1 11 1 1 11 1 1 1 2.1 1.[1 2 3 4 5]/22.1 2 36 5 47 8 91 62 53 43.1 5 69 2 47 8 3xyz=3214.[1 2 3 4]12345.1234[1 2 3 4] 2 20 2.2 MATLAB MATLAB abs angle real imag conj sin cos tan atan atan2 4 exp log log10 rem sign sqrt 2.3 2.3.1 MATLAB 2.2 2 21A =a bc d cos cos (A) =cos (a) cos (b)cos (c) cos (d) cos *1 t=0 : pi/30 : 2*pi plot( t, sin(t))2.3.2 2.3 M=magic(5) sum(M) trace(M) sum(M) sum(rot90(M))rank1 2 34 5 67 8 9 = 2*1 (array) MATLAB (list) arrayable cos 2 22 2.3 MATLAB MATLAB det sum trace norm rank null Kernel eig poly expm logm sqrtm 1 2 34 5 67 8 9xyz =abc2.4 1. 1. 2x+ 13y 3z = 7, x+ y = 1, x+ 7z = 222. x 2y 12z = 12, 2x+ 2y + 2z = 4, 2x+ 3y + 4z = 33. x+ y + z = 5, x y z = 4, 2x+ 6y + 6z = 124. x+ 3y + z = 4, x y + z = 1, 2x+ 4y = 05. x+ 2y + z 4w + 1 = 0, x+ 2y z + 2w 1 = 0,2x+ 4y + z 5w + 1 = 0, x+ 2y + 3z 10w + 2 = 0 2 232. 78 120 264 3. 1.1 2 34 5 67 89 01 22.1 20 31 41 30 11 33.1 1 11 0 21 1 10 1 11 1 22 0 24.7 11 02 3545.7 11 02 3546.1 1 11 0 21 1 10 1 11 1 22 0 27. A =1 2 11 4 12 4 0, P =1 2 11 1 11 1 0 find PAP1 2 24RR RRRRRRRRRRb' a'a b 2.1 (a)a'a(b)a'aRCG1L1C RL2L2C 2.2 4. A=[1: 5] B =1 2 3 4 51 2 3 4 51 2 3 4 51 2 3 4 51 2 3 4 5, C =1 2 3 4 52 4 6 8 103 6 9 12 154 8 12 16 205 10 15 20 255. 2.1 aa, bb, ab, ab 6. 2.2(a), (b) ab 1 7. rank null 8. 2.3 abcd 2 25a b c d1 2 3 4 2.3 4 26 3Visualization 3.1 2 2 3.1.1 plot 3.1 1. t=0 : pi/30 : 2*pi ; plot(t, sin(t))2. t=0 : pi/30 : 2*pi ; plot(cos(t), sin(t))3. t=0 : pi/30 : 2*pi ; plot(t, cos(t), t, sin(t))4. t=linspace(0 , 2*pi, 30); plot(t, cos(t),r , t, sin(t), b*) plot plot(, ), plot(, , r )r red r 3.1 linspace(xmin, xmax, ) [xmin, xmax] 3 27 3.1 y . m o c x x r + + g * b - w : k -. hold on plot(t, sin(3*t)) hold off 3.1 1. y = exp(0.3t) sin(3t), 0 < t < 102. y = 2 sin(t) + cos(4t), 0 < t < 23. y = 2 sin(t) + cos(4t) and y = 2 cos(t) + sin(3t), 0 < t < 24. x = sin(t) and y = cos(t) + sin(3t), 0 < t < 2 axis square axis normal title Exercise 1 xlabel time ylabel voltage 3 28 3.2 2 plot(x,y) ploar(t,r) semilogx(x,y) loglog(x,y) stem(x,y) grid on/off clf axis([xmin xmax ymin ymax]) 3 plot(t, y1, t, y2, t, y3) grid ongrid off clf3.1.2 polar 3 29 3.2 1. t=0 : pi/50 : 2*pi ; polar(t, exp(-0.1*t))2. t=0 : pi/50 : 2*pi ; polar(cos(t), sin(t))3. t=0 : pi/50 : 2*pi ; polar(cos(t), 2*sin(t)+cos(4*t)) polar polar()3.1.3 semilogx 3.3 2 G (s) =1s2 + 2s+ 1y = 20 log |G| = 10 log{(1 x2)2 + (2x)2} x=logspace(-1 , 1, 100); zeta=.1 y=-10*log10((1.0-x.2).2+(2.0*zeta*x).2); semilogx(x, y, r) grid on logspace( 10 10 ) [101, 101] 100 semilogx zeta 3 30RootFigureAxes Uicontrol UimenuSurface Text LightImage Line Patch 3.1 hndl=semilogx(x, y, r) set(hndl, Color, blue) set(hndl, LineWidth, 4) 3.2 3.3 = tan1 2x1 x2 tan1 atan2 atan2( y x ) y=-180*atan2(2.0*zeta*x, 1.0-x.2)/pi;3.1.4 MATLABVersion 5 (objet oriented programming: OOP) (handling) 3.1 1. 3 3 31 figure figure figure3 Figure No. 3 2. figure(1)3. axis([-1,1,-2,2])Figure No. 1 figure 4. 3.1 MATLAB MATLAB 5. Uicontrol(user interface control) Uimenu(user interface menu) Demo help 3.1.5 Bode subplot(m,n,p) m n p x=logspace(-1 , 1, 100); y=-10*log10((1.0-x.2).2+(0.2*x).2); z=-180*atan2(0.2*x, 1.0-x.2)/pi; subplot(2,1,1); semilogx(x, y, r);grid on subplot(2,1,2); semilogx(x, z, b);grid ony -200 axis([0.1,10,-180,0]) (tick) 3 32 set(gca, ytick, [-180:30:0]) gca graphics current axes gcf Root Current Figure gca Figures Current Axes gco Figures Current Object get(gca) get(gcf)set help 3.2 3 3 3D 2D 2D3.2.1 plot3 3.4 (torus) x = (r1 + r2 cos) cos , y = (r1 + r2 cos) sin , z = r2 sin t=0:pi/20:60*pi; r=2+0.5*cos(t/30); x=r.*cos(t) y=r.*sin(t) z=0.5*sin(t/30); plot3(x,y,z,r):axis([-2,2,-2,2,-1,1]):grid onplot3 3D 3 plot( x , y , z , r ) 3 33r red 3.2.2 3.5 x, y z = sinx cos y t=-pi:pi/20:pi; x=t; y=t; z=sin(x)*cos(y); mesh(x,y,z) surf(x,y,z) mesh(z) surf(z) mesh( i, j, zij ) 3 343.2.3 3.6 colormap cool; k=4;n=2k-1; theta=pi*[-n:2:n]/n; phi=(pi/2)*[-n:2:n]/n; x=cos(phi)*cos(theta); y=cos(phi)*sin(theta); z=sin(phi)*ones(size(theta)); c=hadamard(2k); surf(x,y,z,c); axis square colorbar mesh surf (xij , yij , zij) (cij) 3.3 shading interp colormap hsv light lighting gouraud *1 help shading help lighting 3D *1 gouraud 3 35 3.3 hsv hue-saturation-value(HSV)jet a variant of HSVhot black-yellow-whitecool shades of cyan and magendagrey linear gray-scalebone gray-scale with tinge of bluecopper linear copper-toneflag alternating red, white, blue, and blackpink pastel shades of pinkwhite all white color mapprism prism color mapsummer shades of green and yellowautum shades of red and yellowwinter shades of blue and greenspring shades of magenda and yellow3.3 1. polar r = |a sinn + b sin 3n + c sin 5n|a = 7, b = 1, c = 2, n = 1.5 3 36 CG, 63 p. 16 2. 3.3 Bode Nyquist 3. 3.4 3.6 %Torus by surfcolormap cool;k=5;n=2^k-1;theta=(2*pi/3)*[-n:2:n]/n;phi=(pi)*[-n:2:n]/n;r1=1;r2=0.3;a=ones(size(phi));x=(r1*a+r2*cos(phi))*cos(theta);y=(r1*a+r2*cos(phi))*sin(theta);z=r2*sin(phi)*ones(size(theta));c=hadamard(2^k);surf(x,y,z,c);axis([-1,1,-1,1,-1,1]);axis squarecolorbarfigure(1)4. 3 2 5. 6. 7. 1 IIMATLAB 33 44.1 MATLAB MATLAB 2 1 MATLAB (script) 1 (function)MATLAB sin(theta) MATLAB 1 mMATLAB M - file m myprog.m MATLAB myprogMATLAB 4 34FileSave As 4.1 1 100 1. sum([1:100])% % % This is my first program : the addition from 1 to 100sum([1:100])wa100.mMATLAB wa100 MATLAB C 2.% This is my second program : the addition from 1 to 100wa=0;for i=1:100wa=wa+i;endwafor endwafor.mwa3.1 n wa3function a=wa3(n)% This is my third program : the addition from 1 to na=0;for i=1:na=a+i;end%% Of course you can write shortly as follows:% a=sum([1:n]) 4 35%function a=wa3(n) a a = wa3(n) 4. (recursive programming) recwa.mfunction a=recwa(n)% This is my first recursive program : the addition from 1 to nif n==1a=1;return;elsea=n+recwa(n-1);end n=1 a=1 n=n a=n+recwa(n-1) n 4.2 a=b+c; 2 for end while end 4 36 1 3 (operator) >, 4 37if else end if end if 1else 2end if 1 1elseif 2 2else 3end 4.2 2 2 2 A, B wa seki function C=WaorSeki(A,B,ch)% if ch==w then C=A+B and if ch==s then C=A*B.if ch==wC=A+B;elseif ch==sC=A*B;elsedisp(Cannot calculate!!);end w was sekiCannot calculate!! switch case end help switch 4 384.4 4.3 y 1. function TI()% turtle initialization: save as "TI.m"global U X Path degU=[0;1];X=[0;0];Path=X;deg=pi/180;U y X Path deg radian (global)2. function R(a)% rotation by a degree: save as "R.m"global U X Path degtheta=a*deg;U=[cos(theta) sin(theta); -sin(theta) cos(theta)]* U;3. function F(s)% forward s length: save as "F.m"global U X Path degPath=Path+s*U;X=[X Path];4. function J(s) 4 39% jump s length: save as "J.m"global U X Path degPath=Path+s*U;X=[X NaN Path];5. function ST()% show graphics:save as "ST.m"global U X Path degplot(X(1, :), X(2, :));axis square; axis off;set(gcf, Color,[1 1 1]);figure(1);6. No. 1%TI;for i=1:20for j=1:10F(1);R(108);endR(18);endST;7. No. 2Tinit;for i=1:10for j=1:10F(1);R(108);endR(36);endST;8. branch.m rectree function branch(length, level)% save as "branch.m"if level==1return;elseF(length);R(45); 4 40 4.1 branch(length/2, level -1);R(-90);branch(length/2, level -1);R(45);F(-length);endfunction rectree(length, level)% recursive tree:save as "rectree.m"TI;branch(length, level);ST; rectree(10, 5) 1 M-file subfunction function op()% -----------------------------------------------------------% - optical art: coded by H. Kawakami Oct. 23, 1998% - main program: This main program is replced by the problem% - you want to solve. 4 41% -----------------------------------------------------------global U X Path deg parityTI(1);for i=1:6parity=-parity; A(1); B(1,15); H([0;0]);if parity==1R(60*i);elseR(60*i+60);endendST(1)% -----------------------------------------------------------% - subfunctions A(a) and B(a, b) are used% - in the main program op()% -----------------------------------------------------------function A(a)% draw triangleglobal U X Path deg parityfor i=1:3F(a); R(parity*120);endfunction B(a,b)% draw small trianglesglobal U X Path deg parityc=a;for i=1:bR(parity*125); c=0.956*c; F(c); R(parity*120); c=0.9*c;F(c); R(parity*120); F(c);end% -----------------------------------------------------------% - The following subfunctions are commonly used% - for turtle graphics% -----------------------------------------------------------function TI(q)% turtle initializationglobal U X Path deg parityU=[0;1]; X=[0;0]; Path=X; deg=pi/180; parity=1;function ST(q)% show turtleglobal U X Path deg parityplot(X(1,:),X(2,:)); axis square; axis off;set(gcf,Color,[1 1 1]); figure(1)function R(r)% rotation by r degreeglobal U X Path deg 4 42theta=r*deg; U=[cos(theta) sin(theta);-sin(theta) cos(theta)]*U;function F(s)% forward s lengthglobal U X Path degPath=Path+s*U; X=[X Path];function H(h)% return to home positionglobal U X Path degU=[0;1]; Path=h; X=[X [NaN; NaN] Path];% -----------------------------------------------------------4.5 4.5.1 2 (van derPol) d2xdt2 (1 x2) dxdt+ x = 0 (4.1) van der Pol 3 pp. 62-66 (4.1) (4.1) dxdt= ydydt= (1 x2) y x(4.2) 2 (x, y) x y (x0, y0) 4 434.5.2 MATLAB ode454 5 (Runge-Kutta method) % van der Pol equation: the first essay.% save as vdP1.mt0=[0, 20*pi];x0=[0.1; 0.0];[t, x]=ode45(myvdP, t0,x0);plot(x(:, 1), x(:, 2));axis square% end of programmyvdPfunction dx=myvdP(t,x)% van der Pol equation% save as myvdP.mdx=[x(2); 0.5*(1-x(1)*x(1))*x(2)-x(1)];% end of function vdP1% van der Pol equation: the second essay.% save as vdP2.mt0=[0, 20*pi];x0=[0.1; 0.0];[t, x]=ode45(myvdP, t0,x0);% phase portraitfigure(1);plot(x(:, 1), x(:, 2));axis([-2.5, 2.5, -2.5, 2.5]);axis square% wave forms of x(t) and x-dot(t)figure(2); 4 44subplot(2,1,1);plot(t,x(:, 1), r);axis([0, 20*pi, -2.5, 2.5]);subplot(2,1,2);plot(t,x(:, 2), r);axis([0, 20*pi, -2.5, 2.5]);% end of program No. 1 (phase portrait) No.2 x(t) dx(t)/dt *1 (4.2) myvdP 0.5 myvdP eps eps myeps 4.5.3 ode45function vandP2% -----------------------------------------% main program% phase portrait of van der Pol equation% -----------------------------------------global hndlt=0; tmax=20*pi; h=0.1; x=[0.1; 0.0]; X=[x x];GI(1);while t 4 45% van der Pol equationfunction dx=vP(t,x)dx=[x(2); 0.5*(1.0-x(1)*x(1))*x(2)-x(1)];% 4th order Runge Kutta methodfunction [t, x]=RK(t,x,h)f1=vP(t,x);f2=vP(t+h/2,x+h*f1/2);f3=vP(t+h/2,x+h*f2/2);f4=vP(t+h,x+h*f3);t=t+h;x=x+h*(f1+2*(f2+f3)+f4)/6;% Graphics initializationfunction GI(g)global hndlhf=figure(Units,Normalized,Position,[.2 .2 .6 .6]);hndl=line(color,r,linestyle,-,linewidth,1,...erase,none,xdata,[],ydata,[]);axis([-3 3 -3 3]); box on;title(Phase Portrait); xlabel(x(t)); ylabel(x-dot);grid on; axis square;figure(1);4 h line 2 X GI() line linewidth line figure axes x main program 10 GI() while main % main program styleInit_Prob % Initialization of the problemInit_Graf % Initialization of graphicsComputation % Main routineEnd_of_Prog % Termination of program 4 global 4 464.5.4 Rossler 3 Rossler dxdt= x zdydt= x+ aydzdt= bx+ xz cz(4.3) xz Rossler band *2 (4.3) a = 0.35, b = 0.4, c = 4.5x0 = z0 = 0.0, y0 = 3.06 x 10, 8 y 8, 0 z 8function Rossler% main programglobal hndl1 hndl2 hndl3t=0; tmax=200*pi; h=0.1; x=[0.0;3.0; 0.0];L=20; X=x*ones(1,L);GI(x);while t 4 47end%end of main% Rossler equationfunction dx=fn(t,x)dx=[-(x(2)+x(3)); x(1)+0.35*x(2); 0.4*x(1)+x(1)*x(3)-4.5*x(3)];% Runge Kutta methodfunction [t, x]=RK(t,x,h)f1=fn(t,x);f2=fn(t+h/2,x+h*f1/2);f3=fn(t+h/2,x+h*f2/2);f4=fn(t+h,x+h*f3);t=t+h;x=x+h*(f1+2*(f2+f3)+f4)/6;% Graphics initializationfunction GI(x)global hndl1 hndl2 hndl3hf=figure(Units,Normalized,Position,[.2 .2 .6 .6]);hndl1=line(color,r,linestyle,.,markersize,20,...erase,xor,xdata,x(1),ydata,x(2),zdata,x(3));hndl2=line(color,r,linestyle,-,linewidth,1,...erase,none,xdata,[],ydata,[],zdata,[]);hndl3=line(color,b,linestyle,-,...erase,none,xdata,[],ydata,[],zdata,[]);axis([-6 10 -8 8 0 8]);box on;title(Phase Portrait: Rossler equation);xlabel(x(t)); ylabel(y(t)); zlabel(z(t));grid on; axis square;figure(1);3 zdata L hndl1,hndl2, hndl3 GI(x) x MATLAB Lorenz subfunction4.6 MAT LAB 2 4 48 for end 2 MATLAB i=0;for t=0:0.01:10i=i+1;y(i)=sin(t);endt=0:0.01:10;y=sin(t);2 y=zeros(1,100);for i=1:100y(i)=det(X^i);end y Mathematica C MATLAB 4 494.7 1. Hello MATLAB World 2. 3 3. 4.4 optical art 6 8 12 16 4. 4.5.2 van der Pol (t, x, y) 3 5. LRC LRC dxdt= ydydt= x ky +B cos (t+ )(4.4)k B van der Pol Rossler 6. (x, y) xn+1 = yn + 0.008(1 0.05y2n)yn + F (xn)yn+1 = xn + F (xn+1)(4.5)F (x) = x+2 (1 )x21 + x2(4.6) = 0.8 Gumowski-Mira 2 xn+1 = axn + ynyn+1 = x2n + b(4.7) a = 0.1, b = 1.7 50 5 GUI MATLAB GUI(Graphical User Interface) MATLAB OOP(object oriented programming) (event driven) 5.1 MATLAB 2 axes (GUI: Graphical User Interface) Uicontrol Uimenu (event) (event driven) *1 MATLAB 5.1 5.1 (instance) *1 (simple) 5 GUI 51RootFigureUimenuSurfaceTextLightImageLinePatchPush buttonRadio buttonCheck boxList boxEditabletextStatictextPop-upmenuSliderFrameAxes UicontrolGraphics Objects Graphical User-Interface Objects 5.1 5.1 figure uicontrol GUI uimenu Figure axes image 2-D light line patch surface text 5 GUI 52*21 1 (handle) *3 root 0, figure 1 hf=figure hf= 1 ha=axes ha= 13.0006 hf, ha figure hl=line hl= 20.0004 line figure axes line (current) current figure, current axes, current object gcf % returns the value of the current figures handle gca % returns the value of the current axes handle gco % returns the value of the current objects handlegco hndl1, hndl2, ... *2 (: ) *3 5 GUI 53 5.1 figure window axes axes figure window axes *4% Two figure windows each of which contains two axeshf1=figure(Units,Normalized,Position,[.05 .3 .4 .4]);hf2=figure(Units,Normalized,Position,[.50 .3 .4 .4]);ha1=axes(Parent,hf1,Position,[.05 .3 .4 .4]);ha2=axes(Parent,hf1,Position,[.55 .3 .4 .4]);ha3=axes(Parent,hf2,Position,[.05 .3 .4 .4]);ha4=axes(Parent,hf2,Position,[.55 .3 .4 .4]);%==== example of plots ====t=0:pi/30:2*pi;axes(ha1);plot(t,sin(t),Parent,ha1);axis([0 2*pi -1 1]);axes(ha2);plot(t,sin(t).^2,Parent,ha2);axis([0 2*pi -1 1]);axes(ha3);plot(t,cos(t),Parent,ha3);axis([0 2*pi -1 1]);axes(ha4);plot(t,cos(2*t),Parent,ha4);axis([0 2*pi -1 1]);% end of program(property) *5 get set get(, )*4 gcf gca *5 (message) 5 GUI 54 set(, )get 5.2 Hello MATLAB world! for %hello MATLABcenter=0.5; dx=0.6; dy=0.2;x0=center-dx/2; y0=center-dy/2;for i=1:9,r=rand(1,3);s=ones(1,3)-r;hf=figure(Color,[r(1),r(2),r(3)]);set(hf,Position,[250+20*i,300-20*i,300,200]);ha=axes(Position,[x0,y0,dx,dy]);set(ha, Visible,off); %axis off;ht=text(center,center, Hello MATLAB world!);set(ht,Color,[s(1),s(2),s(3)]);set(ht, HorizontalAlignment,Center);end% end of program get n root figure figure 5.2 *6 Units pixels, normalized, centimeters, inches, points root figure (default) pixelsaxes normalizednormalized 1 axis Position [left, bottom,width, height] figure pixels axes text normalized set(gcf, Units, Normalized);get *6 5 GUI 55Root WindowFigure Window(0, 0)(0, 0)(1024, 768)512 pixels384 pixelsLeftBottomHeightWidth0.50.5(1.0, 1.0)(256, 192)(0.25, 0.25) 5.2 root figure 2 5.3 grid on-off Uicontrol 2 % grid on off control by push buttont=0:pi/30:10; r=rand(1,3); s=ones(1,3)-r;hf=figure(Color,[r(1),r(2),r(3)], Units,Normalized,...Position,[.2 .2 .6 .6]);ha=axes(Position,[0.2 0.3 0.6 0.6],Box,on);set(ha,XColor,[s(1),s(2),s(3)],YColor,[s(1),s(2),s(3)]);hb=plot(t,sin(t));set(hb,Color,[s(1),s(2),s(3)],LineStyle,o);set(hf,DefaultUicontrolUnits,Normalized);Uicontrol(hf,Style,Pushbutton,Position,[.25 .1 .2 .1],...Callback,grid on, String,grid on);Uicontrol(hf,Style,Pushbutton,Position,[.55 .1 .2 .1],...Callback,grid off, String,grid off);% end of program 5 GUI 56 5.2 ButtonDownFcn ChangeFcn Children Clipping CreateFcn DeleteFcn BusyAction HandleVisibility Interruptible Parent Selected Tag Type UserData Visible String Callback *7 (callbackroutine) (inheritance) 5.2 *7 M- 5 GUI 575.2 MATLAB 3 get() set() help set 5.2.1 Figure Figure 9 *81. Style and appearance: Menubar, Name, NumberTitle, Resize, Visible, WindowStyle2. General information: Children, Parent, Position, Tag, Type, Units, UserData, Color3. Colormap: Colormap, DitherMap, DitherMapMode, FixedColors, MinColorMap, ShareCol-ors4. Rendering graphics objects: BackingStore, Renderer5. Current selections: CurrentAxes, CurrentCharacter, CurrentMenu, CurrentObject, Current-Point, SelectionType6. Callback routine and execution: ButtonDownFcn, ChangeFcn, CloseRequestFcn, Delete-Fcn, KeyPressFcn, ResizeFcn, BusyAction, Interruptible, WindowButtonUpFcn, Window-ButtonDownFcn, WindowButtonMotionFcn7. Pointer defintion: Pointer, PointerData, PointerShapeHotSpot8. Figure handles: IntegerHandle, HandleVisibility, NextPlot9. Printing: InvertHardcopy, PaperOrientation, PaperPosition, PaperPositionMode, Paper-Size, PaperType, PaperUnits figure hf=figure*8 Using MATLAB Graphics, Chapter 9. 5 GUI 58hf =1 get(hf)BackingStore = onCloseRequestFcn = closereqColor = [0.8 0.8 0.8]Colormap = [ (64 by 3) double array]CurrentAxes = []CurrentCharacter =CurrentObject = []CurrentPoint = [0 0]Dithermap = [ (64 by 3) double array]DithermapMode = manualFixedColors = [ (3 by 3) double array]IntegerHandle = onInvertHardcopy = onKeyPressFcn =MenuBar = figureMinColormap = [64]Name =NextPlot = addNumberTitle = onPaperUnits = inchesPaperOrientation = portraitPaperPosition = [0.25 2.5 8 6]PaperPositionMode = manualPaperSize = [8.5 11]PaperType = usletterPointer = arrowPointerShapeCData = [ (16 by 16) double array]PointerShapeHotSpot = [1 1]Position = [120 120 560 420]Renderer = paintersRendererMode = autoResize = onResizeFcn =SelectionType = normalShareColors = onUnits = pixelsWindowButtonDownFcn =WindowButtonMotionFcn =WindowButtonUpFcn =WindowStyle = normalButtonDownFcn =Children = []Clipping = onCreateFcn =DeleteFcn =BusyAction = queueHandleVisibility = onInterruptible = on 5 GUI 59Parent = [0]Selected = offSelectionHighlight = onTag =Type = figureUserData = []Visible = on set { } set(hf)BackingStore: [ {on} | off ]CloseRequestFcnColorColormapCurrentAxesCurrentObjectCurrentPointDithermapDithermapMode: [ auto | {manual} ]IntegerHandle: [ {on} | off ]InvertHardcopy: [ {on} | off ]KeyPressFcnMenuBar: [ none | {figure} ]MinColormapNameNextPlot: [ {add} | replace | replacechildren ]NumberTitle: [ {on} | off ]PaperUnits: [ {inches} | centimeters | normalized | points ]PaperOrientation: [ {portrait} | landscape ]PaperPositionPaperPositionMode: [ auto | {manual} ]PaperType: [ {usletter} | uslegal | a3 | a4letter | a5 | b4 | tabloid ]Pointer: [ crosshair | fullcrosshair | {arrow} | ibeam | watch |topl | topr | botl | botr | left | top | right | bottom | circle |cross | fleur | custom ]PointerShapeCDataPointerShapeHotSpotPositionRenderer: [ {painters} | zbuffer ]RendererMode: [ {auto} | manual ]Resize: [ {on} | off ]ResizeFcnShareColors: [ {on} | off ]Units: [ inches | centimeters | normalized | points | {pixels} ]WindowButtonDownFcnWindowButtonMotionFcnWindowButtonUpFcnWindowStyle: [ {normal} | modal ]ButtonDownFcnChildren 5 GUI 60Clipping: [ {on} | off ]CreateFcnDeleteFcnBusyAction: [ {queue} | cancel ]HandleVisibility: [ {on} | callback | off ]Interruptible: [ {on} | off ]ParentSelected: [ on | off ]SelectionHighlight: [ {on} | off ]TagUserDataVisible: [ {on} | off ] 5.4 MATLAB (pointer) 16 nargin ptrdemo 0 set(gcf,WindowButtonMotionFcn,ptrdemo arg); arg else 16 16 1 2 NaN function ptrdemo(arg)% Pointers demoptext = {crosshair,fullcross,arrow,watch;...topl,topr,botl,botr;...left,top,right,bottom;...circle,cross,fleur,custom};if nargin == 0clf reset;g = [.5 .5 .5];ha = axes(Units,Normal,Position,[0 0 1 1],...GridLineStyle,-,Xcolor,g,Ycolor,g);set(ha,Color,[.6 1 .6],Xtick,[1 2 3],Ytick,[1 2 3],...XLim,[0 4],YLim,[0 4]);grid on;% define a custom pointerP=ones(16)+1; P(1,:)=1;P(16,:)=1;P(:,1)=1;P(:,16)=1;P(1:4,8:9)=1;P(13:16,8:9)=1;P(8:9,1:4)=1;P(8:9,13:16)=1;P(5:12,5:12)=NaN;set(gcf,Pointer,custom,PointerShapeCData,P,... 5 GUI 61PointerShapeHotSpot,[9 9]);for y = 1:4for x = 1:4ht=text(x-0.95,y-1,ptext{y,x});set(ht,FontSize,16,VerticalAlignment,bottom);endendset(gcf,WindowButtonMotionFcn,ptrdemo arg);set(gcf,Units,Normal);elsept = get(gcf,CurrentPoint);ind = floor(pt * 4) + 1;x = ind(1);y = ind(2);if (x > 0) & (x < 5) & (y > 0) & (y < 5)set(gcf,Pointer, ptext{y,x});endend% end of program5.2.2 Axes Axes 10 *91. Style and appearance: Box, Clipping, GridLineStyle, Layer, LineStyleOrder, LineWidth,MinorGridLineStyle, SelectionHighlight, TickDir, TickDirMode, TickLength, Visible2. General information: Children, ClippingMode, CurrentPoint, Parent, Position, Selected,Tag, Type, Units, UserData3. Annotation: FontAngle, FontName, Fontsize, FontUnits, FontWeight, Title, XLabel, YLa-bel, ZLabel, XTickLabel, YtickLabel, ZTickLabel, XTickLabelMode, YTickLabelMode,ZTickLabelMode,4. Axis control: XDir, YDir, ZDir, XGrid, YGrid, ZGrid, XLim, YLim, ZLim, XLimMode,YLimMode, ZLimMode, XScale, YScale, ZScale, XTick, YTick, ZTick, XTickMode, YTick-Mode, ZTickMode, XAxisLocation, YAxisLocation5. Viewpoint: CameraPosition, CameraPositionMode, CameraTarget, CameraTargetMode,CameraUpVector, CameraUpVectorMode, CameraViewAngle, CameraViewAngleMode,View6. Scaling and aspect ratio: DataAspectoRatio, DataAspectoRatioMode, PlotBoxAspectRatio,PlotBoxAspectRatioMode, ProjectionType,7. Callback execution: BusyAction, ButtonDownFcn, CreateFcn, DeleteFcn, Interruptible8. Rendering method: DrawMode*9 Using MATLAB Graphics, Chapter 10. 5 GUI 629. Targeting axes: HandleVisibility, NextPlot10. Color: AmbientLightColor, CLim, CLimMode, Color, ColorOrder, XColor, YColor, ZColor Axes ha=axes get(ha) set(ha) 5.1 Axes 3.3 3.2 2 1 y function bode% Bode diagramx=logspace(-1,1,100);zeta=.1;y=-10*log10((1.0-x.^2).^2+(2.0*zeta*x).^2);z=-180*atan2(2.0*zeta*x,1.0-x.^2)/pi;hs1=semilogx(x,y,b);ha1=gca;set(ha1,Position,[.1 .2 .8 .7]);title(Bode Diagram);xlabel(Angular Frequency);ylabel(Gain[dB]);grid on;ha2=axes(Position,get(ha1,Position),...YAxisLocation,right,Color,none,...XColor,k,YColor,k,XScale,log,YLim,[-180,0]);hs2=line(x,z,Color,r,Parent,ha2);ylimits=get(ha2,YLim);yinc=(ylimits(2)-ylimits(1))/6;set(ha2,YTick,[ylimits(1):yinc:ylimits(2)]);ylabel(Phase[degree]);figure(1);% end of program5.3 GUI (GUI) uicontrol Style Frame, Text, Edit, Slider, Popupmenu, Listbox,Radiobutton, Checkbox, Pushbutton 5 GUI 63x-axis:y-axis:OKClearCloseAxes TextEditFigurePush button 5.3 5.4 uimenu errordlg, helpdlg, questdlg, warndlg, listdlg, msgbox, choices uigetfile, uiputfile, printdlg, uisetfont, uisetcolor miscellaneous: gtest, ginput, rbbox, waitforbuttonpress 5.5 5.3 Edit Edit OK 5 GUI 64 Clear Close 2 function ShowPts% First trial for Example 5.4% This program has a callback routine:% GetAndDrawE which is stored in% separate file. Nov. 1, 98, HKglobal xm ym hndl1xm=[-3,3]; ym=[-3,3]; %Range of CoordinateInitGUI; % Graphics initializationInitG; % Graphics initializationset(gca,ButtonDownFcn,GetAndDrawE);%end of main% Graphics initializationfunction InitGUIglobal xm ym hndl1hf=figure(Units,Normalized,Position,[.2 .2 .7 .7]);set(hf,DefaultUicontrolUnits,Normalized);set(hf, Pointer,crosshair);hupb1=uicontrol(hf,Style,Pushbutton,String,Close,...Value,0,Position,[0.775,0.15,0.1,0.1],...Callback,closereq);hupb2=uicontrol(hf,Style,Pushbutton,String,Clear,...Value,0,Position,[0.775,0.3,0.1,0.1],...Callback,delete(line));% Graphics initializationfunction InitGglobal xm ym hndl1hndl1=line(color,r,marker,.,markersize,20,...erase,none,xdata,[],ydata,[]);set(gca,Units,Normalized,Position,[.05 .2 .7 .7]);axis([xm(1) xm(2) ym(1) ym(2)]);box on;title(Rectangular Coordinate);xlabel(x-axis);ylabel(y-axis);grid on;axis square;figure(1);% end of program Callback function GetAndDrawE% This is a callback routine for Example 5.4 5 GUI 65global xm ym hndl1a=get(gca,CurrentPoint);x=a(1,1);y=a(1,2);set(hndl1,xdata,x,ydata,y);drawnow; GUI 3 1. callback routine main callback routine main subfunction callbackroutine file 2. main callback routine main default 1 callback routine file subfunction subfunction demo 3. while subfunction 2 file 1 (1) (2) (3) function ShowPts(arg)% Second trial for Example 5.4% This program has a self-callback routine.% The second style is applied. Nov. 9, 98 HiK.global xm ym hndl1if nargin==0init_all;elseshow_point;end%end of mainfunction init_allglobal xm ym hndl1xm=[-3,3]; ym=[-3,3]; %Range of CoordinateInitGUI; % Graphics initializationInitG; % Graphics initializationset(gca,ButtonDownFcn,ShowPts arg);function show_point 5 GUI 66global xm ym hndl1a=get(gca,CurrentPoint);x=a(1,1);y=a(1,2);set(hndl1,xdata,x,ydata,y);drawnow;% GUI initialization% Graphics initializationfunction InitGUIglobal xm ym hndl1hf=figure(Units,Normalized,Position,[.2 .2 .7 .7]);set(hf,DefaultUicontrolUnits,Normalized);set(hf, Pointer,crosshair);hupb1=uicontrol(hf,Style,Pushbutton,String,Close,...Value,0,Position,[0.775,0.15,0.1,0.1],...Callback,closereq);hupb2=uicontrol(hf,Style,Pushbutton,String,Clear,...Value,0,Position,[0.775,0.3,0.1,0.1],...Callback,delete(line));function InitGglobal xm ym hndl1hndl1=line(color,r,marker,.,markersize,20,...erase,none,xdata,[],ydata,[]);set(gca,Units,Normalized,Position,[.05 .2 .7 .7]);axis([xm(1) xm(2) ym(1) ym(2)]);box on;title(Rectangular Coordinate);xlabel(x-axis);ylabel(y-axis);grid on;axis square;figure(1);3 function ShowPts% Third trial for Example 5.4% This program has an endless while statement.% The third style is applied. Nov. 9, 98 HiK.global xm ym hndl1 bu x yxm=[-3,3]; ym=[-3,3]; bu=-1;x=0;y=0;InitGUI;InitG;set(gca,ButtonDownFcn,bdf)set(gca,Userdata,1);while get(gca,Userdata)==1,if bu==1, 5 GUI 67gcpnt;end;drawnow;endclear global xm ym hndl1 x y;closereq;%end of main% GUI initializationfunction InitGUIhf=figure(Units,Normalized,Position,[.2 .2 .7 .7]);set(hf,DefaultUicontrolUnits,Normalized);set(hf, Pointer,crosshair);hupb1=uicontrol(hf,Style,Pushbutton,String,Close,...Value,0,Position,[0.775,0.15,0.1,0.1],...Callback,set(gca,Userdata,-1));hupb2=uicontrol(hf,Style,Pushbutton,String,Clear,...Value,0,Position,[0.775,0.3,0.1,0.1],...Callback,delete(line));% Graphics initializationfunction InitGglobal xm ym hndl1set(gca,Units,Normalized,Position,[.05 .2 .7 .7]);axis([xm(1) xm(2) ym(1) ym(2)]);hndl1=line(color,r,marker,.,markersize,20,...erase,none,xdata,[],ydata,[]);box on;title(Rectangular Coordinate);xlabel(x-axis);ylabel(y-axis);grid on;axis square;figure(1);function gcpntglobal hndl1 bu x ya=get(gca,CurrentPoint);x=a(1,1);y=a(1,2);set(hndl1,xdata,x,ydata,y);bu=-1;function bdfglobal bubu=1; Editfunction PlotPoint(arg)% Main program for Example 5.4% This program has two self-callback routines:% Point_Draw and Put_Point, Nov. 4, 1998. 5 GUI 68% File Name: PlotPoint.mglobal xm ym huex huey hndl1 hndl2if nargin==0xm=[-3,3]; ym=[-3,3];InitGUI; % GUI initializationInitG; % Graphics initializationset(gca,ButtonDownFcn,PlotPoint PDraw);elseif arg==PDrawPoint_Draw;elseif aeg==PPointPut_Pointend%end of mainfunction Point_Draw% This is a callback routine for ButtonDownFcnglobal xm ym huex huey hndl1 hndl2a=get(gca,CurrentPoint);x=a(1,1); y=a(1,2);cx=num2str(x); cy=num2str(y);set(hndl1,xdata,x,ydata,y);set(huex,String,cx); set(huey,String,cy);drawnow;function Put_Point% This is a callback routine for PushButton:OKglobal xm ym huex huey hndl1 hndl2fx=get(huex,String);x=str2num(fx);if (xxm(2))x=xm(2);end;fy=get(huey,String);y=str2num(fy);if (yym(2))y=ym(2);end;set(hndl2,xdata,x,ydata,y);drawnow;% GUI Initializationfunction InitGUIglobal xm ym huex huey hndl1 hndl2hf=figure(Units,Normalized,Position,[.2 .2 .7 .7]);set(hf,DefaultUicontrolUnits,Normalized);set(hf, Pointer,crosshair);hupb1=uicontrol(hf,Style,Pushbutton,String,Close,...Value,0,Position,[0.775,0.15,0.1,0.1],...Callback,closereq); 5 GUI 69hupb2=uicontrol(hf,Style,Pushbutton,String,Clear,...Value,0,Position,[0.775,0.3,0.1,0.1],...Callback,delete(line));hupb3=uicontrol(hf,Style,Pushbutton,String,OK,...Value,0,Position,[0.775,0.55,0.1,0.1],...Callback,PutDraw);huex=uicontrol(hf,Style,Edit,String,,...Position,[.75 .8 .15 .05],HorizontalAlignment,Left);huey=uicontrol(hf,Style,Edit,String,,...Position,[.75 .7 .15 .05],HorizontalAlignment,Left);hutx=uicontrol(hf,Style,Text,String,x-axis:,...Position,[.68 .8 .06 .05],HorizontalAlignment,Center);huty=uicontrol(hf,Style,Text,String,y-axis:,...Position,[.68 .7 .06 .05],HorizontalAlignment,Center);% Graphics Initializationfunction InitGglobal xm ym huex huey hndl1 hndl2hndl1=line(color,r,marker,.,markersize,20,...erase,none,xdata,[],ydata,[]);hndl2=line(color,g,marker,.,markersize,20,...erase,none,xdata,[],ydata,[]);set(gca,Units,Normalized,Position,[.05 .2 .7 .7]);axis([xm(1) xm(2) ym(1) ym(2)]);box on;title(Rectangular Coordinate);xlabel(x-axis); ylabel(y-axis);grid on; axis square;figure(1);5.4 GUI function Event_driven_like_program% Example of programming style% main partglobal %definition of global variablesCreateWindows;ComputeAndDraw;CloseWindows;%-----------------------------------function CreateWindows;% Creation and initialization of graphics objectsfunction CloseWindows;% Termination of graphics objects%-----------------------------------function ComputeAndDraw% main working part of this programglobal %definition of global variableswhile 1, 5 GUI 70EventLoop;% here put the main computationend;function EventLoop;function HitButton(event);function HitKey(event);%-----------------------------------%define other functions needed for this program% end of program4.5.3 5.6 van der Pol 4.5.3 2 1 1 function PhPt% main programglobal hndl1 hndl2 hndl3 x X bu Lt=0; tmax=20*pi; h=0.1; bu=-1; L=20;x=[0.1; 0.0]; X=x*ones(1, L);InitGUI;InitG(x);set(gca,ButtonDownFcn,bdf)set(gca,Userdata,1);while get(gca,Userdata)==1,if bu==1,gcpntend;[t, x]=RK(t,x,h);X=[x X(:,1: L-1)];set(hndl1,xdata,X(1,1),ydata,X(2,1));set(hndl2,xdata,X(1,1:L-1),ydata,X(2,1:L-1));set(hndl3,xdata,X(1,L-1:L),ydata,X(2,L-1:L));drawnow;end%close(gcf);closereq;%end of main% van der Pol equationfunction dx=vP(t,x)dx=[x(2); 0.5*(1.0-x(1)*x(1))*x(2)-x(1)];% Runge Kutta methodfunction [t, x]=RK(t,x,h) 5 GUI 71f1=vP(t,x);f2=vP(t+h/2,x+h*f1/2);f3=vP(t+h/2,x+h*f2/2);f4=vP(t+h,x+h*f3);t=t+h;x=x+h*(f1+2*(f2+f3)+f4)/6;% Graphics initializationfunction InitGUIset(gcf,DefaultUicontrolUnits,Normalized);set(gcf, Pointer,crosshair);hd1=uicontrol(gcf,Style,Pushbutton,String,Close,...Value,0,Position,[0.85,0.15,0.1,0.1],...Callback,set(gca,Userdata,-1));% set(hd1, BackgroundColor,[.8 .8 .8]);hd2=uicontrol(gcf,Style,Pushbutton,String,Cls,...Value,0,Position,[0.85,0.3,0.1,0.1],...Callback,delete(line));% set(hd2, BackgroundColor,[.1 .1 .1]);% Graphics initializationfunction InitG(x)global hndl1 hndl2 hndl3hndl1=line(color,r,marker,.,markersize,25,...erase,xor,xdata,x(1),ydata,x(2));hndl2=line(color,r,linestyle,-,linewidth,2,...erase,none,xdata,[],ydata,[]);hndl3=line(color,b,linestyle,-,linewidth,2,...erase,none,xdata,[],ydata,[]);axis([-3 3 -3 3]);box on;title(Phase Portrait);xlabel(x(t));ylabel(x-dot);grid on;axis square;figure(1);function gcpntglobal hndl1 hndl2 hndl3 x X bu La=get(gca,CurrentPoint);x=[a(1,1); a(1,2)];X=x*ones(1,L);bu=-1;function bdfglobal hndl x X bubu=1;% end of program 5 GUI 725.5 MATLAB MATLAB guide GUI MATLAB Property Editor Callback Editor Visual Basic GUI human interface Simplicity, Consistency andFamiliarity demo 5.6 1. 5.2 sin, cos, tan 2. 3. 4. 5. 5 GUI /2 1 6. 3D function az(button)global AZ_PNTif nargin==0button=on;endif strcmp(button,on)set(gcf,WindowButtonDownFcn,az down);set(gcf,WindowButtonMotionFcn,az motion);set(gcf,WindowButtonUpFcn,az up);AZ_PNT=[];elseif strcmp(button,off)set(gcf,WindowButtonDownFcn,);set(gcf,WindowButtonMotionFcn,);set(gcf,WindowButtonUpFcn,);clear global AZ_PNTelseif strcmp(button,down) 5 GUI 73AZ_PNT=get(gcf,CurrentPoint);elseif strcmp(button,motion)if isempty(AZ_PNT)==0new_point=get(gcf,CurrentPoint);x=[AZ_PNT(1),new_point(1)];y=[AZ_PNT(2),new_point(2)];[azimuth,zenith]=view;position=get(gcf,Position);AZ_WINDOW=position(3:4);azimuth=azimuth-360.*(x(2)-x(1))/AZ_WINDOW(1);zenith=zenith-360.*(y(2)-y(1))/AZ_WINDOW(2);view([rem(azimuth,360),rem(zenith,360)])AZ_PNT=get(gcf,CurrentPoint);endelseif strcmp(button,up)AZ_PNT=[];end IIIMATLAB 75 61987 62 3 Fortran MatLab 6.1 Duffing d2xdt2+ kdxdt+ c1x+ c3x3 = B0 +B cos(t) (6.1)dxdt= ydydt= ky c1x c3x3 +B0 +B cos(t)(6.2) = (k, c1, c3, B0, B) (6.1) (6.2) Newton (first guess) singular non-singular non-singular 1. phase portrait Poincare 2. Newton 6 763. 6.2 6.2.1 dxdt= f(t, x, ) (6.3) x Rn R *1. f : RRn R Rn t 2 x(t0) = x0 t0 = 0 (6.3) x(t) = (t, x0, ) (6.4)Poincare map T T : Rn Rn; x 7 (2, x0, ) (6.5) (6.3) 2mT mF (x) = (2m, x, ) x = 0 (6.6) (6.6) Newton 1. x = x(0) 2. x = x(n+1) x(n) x(n+1) = x(n) + hDF (x(n))h + F (x(n)) = 0(6.7)DF (6.3) 6.2.2 (6.3) T x T T (x) x = 0 (6.8)*1 6 77 6.1 Duffing x = 1(t, x0, y0, ) x1 x1 = x2 x0y = 2(t, x0, y0, ) x2 x2 = kx2 c1x1 c3x3 +B0 +B cos(t) y01/x0 x3 x3 = x4 12/x0 x4 x4 = Px3 kx4 01/y0 x5 x5 = x6 02/y0 x6 x6 = Px5 kx6 11/B x7 x7 = x8 02/B x8 x8 = Px7 kx8 + cos(t) 021/x20 x9 x9 = x10 022/x20 x10 x10 = Px9 kx10 +Qx23 021/x0y0 x11 x11 = x12 022/x0y0 x12 x12 = Px11 kx12 +Qx3x5 021/y20 x13 x13 = x14 022/y20 x14 x14 = Px13 kx14 +Qx25 021/x0B x15 x15 = x16 022/x0B x16 x16 = Px15 kx16 +Qx3x7 021/y0B x17 x17 = x18 022/y0B x18 x18 = Px17 kx18 +Qx3x7 0where P = c1 c3x31 and Q = 6c3x1. The parameter is selected as B. 6 78g() = det(DT (x) In) = 0 (6.9)1. (tangent bifurcation) (6.9) = 1 (6.8) g(1) = 0 (x, ) Newton 2. (period doubling bifurcation) (6.9) = 1 (6.8) g(1) = 0 (x, ) Newton (6.3) (6.2) 6.1 6.3 Phase Portrait DemoPLPhase Portrait 1987 GUI function DemoPL% main programglobal hndl1 hndl2 huex huey x X bu hf B DB B0 DB0 k dk iflglobal huek hueB hueB0B=0.3; DB=0.01; B0=0.0; DB0=0.01; k=0.1; dk=0.01;ifl=1;bu=-1; m=64; h=2.0*pi/m;x=[0.0; 0.0]; X=[x x];InitGUI;InitG(x);set(gca,ButtonDownFcn,bdf)set(gca,Userdata,1);set(hf,KeyPressFcn,KeyInt);while get(gca,Userdata)==1,set(hndl1,xdata,X(1,1),ydata,X(2,1));cx=num2str(X(1,1)); cy=num2str(X(2,1));set(huex,String,cx); set(huey,String,cy);ck=num2str(k); cB=num2str(B); cB0=num2str(B0);set(huek,String,ck); set(hueB,String,cB); set(hueB0, String,cB0);drawnow;for j=0:m-1if bu==1,gcpnt;break; 6 79end;t=j*h;[t, x]=RK(t,x,h);X=[x X(:,1)];if ifl==1,set(hndl2,xdata,X(1,1:2),ydata,X(2,1:2));end;endendclosereq;%end of main% Duffings equationfunction dy=Duf(t,y)global B B0 kdy=[y(2); -k*y(2)-y(1)*y(1)*y(1)+B0+B*cos(t)];% Runge Kutta methodfunction [t, z]=RK(t,z,h)f1=Duf(t,z);f2=Duf(t+h/2,z+h*f1/2);f3=Duf(t+h/2,z+h*f2/2);f4=Duf(t+h,z+h*f3);t=t+h;z=z+h*(f1+2*(f2+f3)+f4)/6;% Graphics initializationfunction InitGUIglobal hndl1 hndl2 huex huey hfglobal huek hueB hueB0hf=figure(Units,Normalized,Position,[.1 .2 .7 .7]);set(hf,DefaultUicontrolUnits,Normalized);set(hf, Pointer,crosshair);hd1=uicontrol(hf,Style,Pushbutton,String,Close,Value,0,...Position,[0.8,0.1,0.1,0.1],Callback,set(gca,Userdata,-1));hd2=uicontrol(hf,Style,Pushbutton,String,Cls,Value,0,...Position,[0.8,0.22,0.1,0.1],Callback,delete(line));hutx=uicontrol(hf,Style,Text,String,x-axis:,...Position,[.8 .87 .08 .04],HorizontalAlignment,Center);huex=uicontrol(hf,Style,Edit,String,,...Position,[.8 .84 .1 .05],HorizontalAlignment,Left);huty=uicontrol(hf,Style,Text,String,y-axis:,...Position,[.8 .77 .08 .04],HorizontalAlignment,Center);huey=uicontrol(hf,Style,Edit,String,,...Position,[.8 .74 .1 .05],HorizontalAlignment,Left);hutk=uicontrol(hf,Style,Text,String,k:,...Position,[.8 .63 .08 .04],HorizontalAlignment,Center);huek=uicontrol(hf,Style,Edit,String,,...Position,[.8 .6 .1 .05],HorizontalAlignment,Left);hutB=uicontrol(hf,Style,Text,String,B:,... 6 80Position,[.8 .53 .08 .04],HorizontalAlignment,Center);hueB=uicontrol(hf,Style,Edit,String,,...Position,[.8 .5 .1 .05],HorizontalAlignment,Left);hutB0=uicontrol(hf,Style,Text,String,B0:,...Position,[.8 .43 .08 .04],HorizontalAlignment,Center);hueB0=uicontrol(hf,Style,Edit,String,,...Position,[.8 .4 .1 .05],HorizontalAlignment,Left);% Graphics initializationfunction InitG(x)global hndl1 hndl2 huex hueyhndl1=line(color,r,marker,.,markersize,20,erase,none,...xdata,[],ydata,[]);hndl2=line(color,b,linestyle,-,linewidth,1,erase,none,...xdata,[],ydata,[]);set(gca,Position,[0.02 0.1 .8 .8]);axis([-2 2 -2 2]);box on;title(Phase Portrait);xlabel(x(t));ylabel(x-dot);grid on;axis square;figure(1);function gcpntglobal x X bua=get(gca,CurrentPoint);x=[a(1,1); a(1,2)];X=[x x];bu=-1;function bdfglobal bubu=1;% Key interrupt functionfunction KeyInt;global hf B DB B0 DB0 k dk ifla=get(hf,CurrentCharacter);switch acase Difl=1-ifl;case bB=B-DB;case BB=B+DB;case cB0=B0-DB0;case CB0=B0+DB0;case kk=k-dk; 6 81case Kk=k+dk;otherwise%disp(other key is pressed);end6.4 DufFixfunction DufFix% main programglobal B B0 kglobal IterMax ErrMax epsxB=0.3; DB=0.01; B0=0.0; DB0=0.01; k=0.1; dk=0.01;L=1; m=128;x=[-0.32; 0.04; 1.0; 0.0; 0.0; 1.0];IterMax=10; ErrMax=100.0; epsx=0.00001;while 1,[s_ind, iterN, x, chara]=Fixed(L, m, x);switch s_indcase 0% convergencedisp([s_ind iterN B B0 x(1) x(2) chara]);case 1% divergencedisp([s_ind States are divergent]);break;case 2% so many iterationsdisp([s_ind So many iterations in the Newton method]);break;otherwise% impossibledisp([This case never occurs]);endB0=B0+DB0;end% end of main program% Find fixed point or L-periodic pointfunction [s_ind, iterN, x, chara]=Fixed(L,m,x)global IterMax ErrMax epsxs_ind=0; dist=1.0; iterN=1; p=zeros(2,1); chara=[];while dist>epsx,iterN=iterN+1;p(1)=x(1); p(2)=x(2);x=NewtonR(L,m,x);ddx=x(1)-p(1); ddy=x(2)-p(2);dist=(abs(x(1)-p(1))+abs(x(2)-p(2)))/2; 6 82if dist>ErrMaxs_ind=1;% divergencebreak;endif iterN>IterMaxs_ind=2; % sorry iterating too muchbreak;endend% characteristic multipliers[V, W]=eig([x(3) x(4); x(5) x(6)]);chara=[W(1,1) W(2,2)];% Newtons iterationfunction x=NewtonR(L,m,x)p=x(1:2); A=zeros(2,2);x=sysvar(L,m,x);A(1,1)=x(3)-1.0;A(1,2)=x(5);A(2,1)=x(4);A(2,2)=x(6)-1.0;x(1:2)=p+A\[p(1)-x(1); p(2)-x(2)];% Poincare map of statesfunction x=sysvar(L,m,x)x(3)=1.0; x(4)=0.0; x(5)=0.0; x(6)=1.0;h=2.0*pi/m;for kkk=1:Lfor j=0:m-1t=j*h;[t, x]=RK(t,x,h);endend% Runge Kutta methodfunction [t, z]=RK(t,z,h)f1=Duf(t,z);f2=Duf(t+h/2,z+h*f1/2);f3=Duf(t+h/2,z+h*f2/2);f4=Duf(t+h,z+h*f3);t=t+h;z=z+h*(f1+2*(f2+f3)+f4)/6;% Duffings equationfunction dy=Duf(t,y)global B B0 kdy=zeros(6,1); c1=0.0; c3=1.0;P= -c1-3.0*c3*y(1)*y(1) ;dy(1)=y(2);dy(2)=-k*y(2)-y(1)*y(1)*y(1)+B0+B*cos(t);dy(3)=y(4);dy(4)=P*y(3)-k*y(4);dy(5)=y(6);dy(6)=P*y(5)-k*y(6); 6 836.5 DemoBF DufFix function DufBf% main programglobal B B0 kglobal IterMax ErrMax epsxL=1; m=128; IterMax=20; ErrMax=100.0; epsx=0.00001;% A0 initial guess%B=0.46; DB=0.01; B0=0.0; DB0=0.01; k=0.2; dk=0.01;%x=[-0.654; 0.224; 1.0; 0.0; 0.0; 1.0;0 ; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];%mu=1.0;B=3.0; DB=0.01; B0=0.4988; DB0=0.01; k=0.2; dk=0.01;x=[1.9632; -2.42; 1.0; 0.0; 0.0; 1.0;0 ; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];mu=-1.0;while 1,[s_ind, iterN, x, chara]=Fixed(L, m, x, mu);switch s_indcase 0% convergencedisp([iterN B B0 x(1) x(2) chara]);case 1% divergencedisp([s_ind States are divergent]);break;case 2% so many iterationsdisp([iterN So many iterations in the Newton method]);break;otherwise% impossibledisp([This case never occurs]);endB0=B0+DB0;end% end of main program% Find fixed point or L-periodic pointfunction [s_ind, iterN, x, chara]=Fixed(L,m,x, mu)global B B0 kglobal IterMax ErrMax epsxs_ind=0; dist=1.0; iterN=1; p=zeros(2,1); chara=[];while dist>epsx,iterN=iterN+1;p(1)=x(1); p(2)=x(2); p(3)=B;x=NewtonR(L,m,x,mu); 6 84ddx=x(1)-p(1); ddy=x(2)-p(2);ddz=(B-p(3))/B;dist=(abs(x(1)-p(1))+abs(x(2)-p(2))+abs(ddz))/3;if dist>ErrMaxs_ind=1;% divergencebreak;endif iterN>IterMaxs_ind=2; % sorry iterating too muchbreak;endend% characteristic multipliers[V, W]=eig([x(3) x(4); x(5) x(6)]);chara=[W(1,1) W(2,2)];% Newtons iterationfunction x=NewtonR(L,m,x,mu)global B B0 kp=[x(1);x(2);B]; A=zeros(3,3);x=sysvar(L,m,x);A(1,1)=x(3)-1.0;A(1,2)=x(5);A(1,3)=x(7);A(2,1)=x(4);A(2,2)=x(6)-1.0;A(2,3)=x(8);A(3,1)=mu*(-x(9)-x(12));A(3,2)=mu*(-x(11)-x(14));A(3,3)=x(15)*(x(6)-mu)-x(16)*x(5)+(x(3)-mu)*x(18)-x(4)*x(17);ppp=p+A\[p(1)-x(1); p(2)-x(2);-mu*mu+mu*(x(3)+x(6))-x(3)*x(6)+x(4)*x(5)];x(1)=ppp(1); x(2)=ppp(2); B=ppp(3);% Poincare map of statesfunction x=sysvar(L,m,x)x(3)=1.0; x(4)=0.0; x(5)=0.0; x(6)=1.0;x(7)=0.0; x(8)=0.0; x(9)=0.0; x(10)=0.0;x(11)=0.0; x(12)=0.0; x(13)=0.0; x(14)=0.0;x(15)=0.0; x(16)=0.0; x(17)=0.0; x(18)=0.0;h=2.0*pi/m;for kkk=1:Lfor j=0:m-1t=j*h;[t, x]=RK(t,x,h);endend% Runge Kutta methodfunction [t, z]=RK(t,z,h)f1=Duf(t,z);f2=Duf(t+h/2,z+h*f1/2);f3=Duf(t+h/2,z+h*f2/2);f4=Duf(t+h,z+h*f3);t=t+h; 6 85z=z+h*(f1+2*(f2+f3)+f4)/6;% Duffings equationfunction dy=Duf(t,y)global B B0 kdy=zeros(6,1); c1=0.0; c3=1.0;P= -c1-3.0*c3*y(1)*y(1) ;Q=-6.0*c3*y(1);dy(1)=y(2);dy(2)=-k*y(2)-y(1)*y(1)*y(1)+B0+B*cos(t);dy(3)=y(4);dy(4)=P*y(3)-k*y(4);dy(5)=y(6);dy(6)=P*y(5)-k*y(6);dy(7)=y(8);dy(8)=P*y(7)-k*y(8)+cos(t);dy(9)=y(10);dy(11)=y(12);dy(13)=y(14);dy(15)=y(16);dy(17)=y(18);dy(10)=P*y(9)-k*y(10)+Q*y(3)*y(3);dy(12)=P*y(11)-k*y(12)+Q*y(3)*y(5);dy(14)=P*y(13)-k*y(14)+Q*y(5)*y(5);dy(16)=P*y(15)-k*y(16)+Q*y(3)*y(7);dy(18)=P*y(17)-k*y(18)+Q*y(5)*y(7);86 77.1 dxdt= f(x, ) (7.1) x x =x1...xn Rn, =1...m Rm (7.2)t = 0 x0 (7.1) x (t) = (t, x0), x0 = (0, x0) (7.3) (7.1) L *1 (local cross section) Poincare T T : ; x 7 (L (x) , x) (7.4)T *2 DT *1 1 (codimension 1) *2 T DT Tx 7 877.2 Poincare (cross section) g : Rn R g (x) = 0 (7.5) = {x Rn | g (x) = 0} (7.6) gxf 6= 0 (7.7) x0 t = L (x0) x1 x1 = x (L(x0)) = (L(x0), x0) (7.8) (7.7) x0 DT (x0) =x1x0=x0+tLx0=x0+ fLx0(7.9)g (x1) = g ( (L (x0) , x0)) = 0 (7.10) x0 gx(x0+ fLx0)= 0 (7.11)Lx0= 1gxfgxx0(7.12) (7.9) DT (x0) =x0 1gxffgxx0=[In 1gxffgx]x0(7.13) Rn1 = {(u1, , un1)} x1 = u1, x2 = u2, , xn1 = un1, xn = h (u1, u2, , un1) (7.14) Rn : Rn1 Rn; u 7 x = (u) = (u1, , un1, h(u1, , un1)) (7.15) 7 88 Rn1 :p : Rn Rn1; x 7 p (x) = (x1, , xn1) (7.16) Rn1 Poincare : Rn1 Rn1; u 7 (u) = pT (u) = u1 (7.17)D (u0) = Dp (x1)DT (x0) (u0) = [In1...0]DT (x0)In1Dh= [In1...0]DT (x0)In1 gugxn(7.18)g (u1, u2, , un1, h (u1, u2, , un1)) = 0h(u) Dh(u) = 1gxngu(7.19) (7.18) D(u) = [In1...0][In 1gx ff gx]x0In1 gugxn(7.20)7.2.1 12 dxdt= f1 (x, y)dydt= f2 (x, y)(7.21)Poincare g (x, y) = y = 0 (7.22)gx= [0, 1] ,dgdx=gx+gyDh = 0 Dh = gx= 0 (7.23) 7 89D(u) =[1 0]1 f1f20 0[x0] 10 =[1 f1f2]x010 (7.24)x0=11 1221 22 (7.25) (7.25) D (u) = 11 f1f221 (7.26)7.2.2 23 dxdt= f1 (x, y, z)dydt= f2 (x, y, z)dzdt= f3 (x, y, z)(7.27)Poincare g (x, y, z) = y = 0 (7.28) (x, z) (x, z) f2 | 6= 0 (7.29)g (x, y, z) = y = 0y = h (x, z) = 0f =f1f2f3(7.30) 7 90 f2 6= 0 Dh = [gxgz]=[0 0],gx=[0 1 0] gx f | = f2f | gx=f1f2f3[0 1 0]=0 f1 00 f2 00 f3 0I3 1f2fgx=1 0 00 1 00 0 10f1f200 1 00f3f20=1 f1f200 0 00 f3f21D(u, v) =1 0 00 1 01 f1f200 0 00 f3f21[x0]1 00 00 1 =1 f1f200 f3f21 x01 00 00 1(7.31)x0=11 12 1321 22 2331 32 33 (7.32) (7.31) D(u, v) =11 f1f221 13 f1f22331 f3f221 33 f3f223 (7.33) 1x = 0 D(u, v) =22 f2f112 23 f2f11332 f3f112 33 f3f113 (7.34) 2z = 0 D(u, v) =11 f1f331 12 f1f33221 f2f331 22 f2f332 (7.35) 7 917.2.3 34 x1 = f1 (x1, x2, x3, x4)x2 = f2 (x1, x2, x3, x4)x3 = f3 (x1, x2, x3, x4)x4 = f4 (x1, x2, x3, x4)(7.36)Poincare g (x) = x4 = 0 (7.37)f4 (x1, x2, x3, x4) 6= 0 (7.38)g (u, h (u)) = 0 (7.39)Dg (u, h (u)) =gu+gx4hu= 0 (7.40)Dh (u) = gu=[0 0 0],gx=[0 0 0 1]I4 1f4fgx=1 0 0 00 1 0 00 0 1 00 0 0 1 1f40 0 0 f10 0 0 f20 0 0 f30 0 0 f4=1 0 0 f1/f40 1 0 f2/f40 0 1 f3/f40 0 0 0(7.41)D(u1, u2, u3) =11 f1f441 12 f1f442 13 f1f44321 f2f441 22 f2f442 23 f2f44331 f3f441 32 f3f442 33 f3f443(7.42) 7 927.3 Rossler 7.3.1 Rossler Poincare Rossler Rossler dxdt= y z = f1dydt= x+ ay = f2dzdt= bx cz + xz = f3(7.43)g (x, y, z) = y = 0 (7.44) (x, z) (x, z) f2 | 6= 0 (7.45) Poincare T [u1, 0, u3] L (7.43) F1(L, u1, u3) = 1(L, u1, 0, u3) u1 = 0F1(L, u1, u3) = 2(L, u1, 0, u3) = 0F1(L, u1, u3) = 3(L, u1, 0, u3) u3 = 0(7.46) (7.46) Newton (7.46) TaylorF + Jh = 0 (7.47)F =F1(L, u1, u3)F2(L, u1, u3)F3(L, u1, u3) , h =Lu1u3 , J =f11x1 1 1x3f22x12x3f33x13x3 1(7.48) (7.47) L u1, u3 Newton 7 93L(n+1)u(n+1)1u(n+1)3 =L(n)u(n)1u(n)3 J1L(n)u(n)1u(n)3 (7.49) J ddt1x12x13x1=0 1 11 a 0b+ 3 0 1 c1x12x13x1,1x12x13x1t=0=100 (7.50)ddt1x32x33x3=0 1 11 a 0b+ 3 0 1 c1x32x33x3,1x32x33x3t=0=001 (7.51) (7.33) 7.3.2 : RosFixfunction RosFix% main program for finding periodic solutions of Rosslers equationglobal a b cglobal IterMax ErrMax epsxa=0.6; da=0.01; b=0.5; db=0.01; c=1.6; dc=0.02;%L=10.9;L=5.65; % the periodm=512;%x=[1.13;0.0;2.64;1.0;0;0;0;0;1];x=[1.44;0.0;1.76;1.0;0;0;0;0;1];IterMax=10; ErrMax=100.0; epsx=0.00001;while 1,[s_ind, iterN, L, x, chara]=Fixed(L, m, x);switch s_indcase 0% convergencedisp([iterN a L x(1) x(3) chara]);case 1% divergencedisp([s_ind States are divergent]);break;case 2% so many iterationsdisp([iterN So many iterations in the Newton method]);break;otherwise 7 94% impossibledisp([This case never occurs]);enda=a+da;end% end of main program% Find fixed point: Poincare section is fixed as x(2)=0function [s_ind, iterN, L, x, chara]=Fixed(L,m,x)global a b cglobal IterMax ErrMax epsxs_ind=0; dist=1.0; iterN=1; p=zeros(3,1); chara=[];while dist>epsx,iterN=iterN+1;p(1)=L; p(2)=x(1); p(3)=x(3);% disp([L x(1) x(2) x(3) dist]);[L, x]=NewtonR(L,m,x);ddx=x(1)-p(2); ddy=x(3)-p(3); ddz=L-p(1);dist=(abs(ddx)+abs(ddy)+abs(ddz))/3;if dist>ErrMaxs_ind=1;% divergencebreak;endif iterN>IterMaxs_ind=2; % sorry iterating too muchbreak;endend% characteristic multipliers[V, W]=eig([x(4)+x(3)*x(5)/x(1), x(7)+x(3)*x(8)/x(1);x(6)-(b*x(1)-c*x(3)+x(1)*x(3))*x(5)/x(1),...x(9)-(b*x(1)-c*x(3)+x(1)*x(3))*x(8)/x(1)]);chara=[W(1,1) W(2,2)];% Newtons iterationfunction [L, x]=NewtonR(L,m,x)global a b cp=[L; x(1); x(3)]; A=zeros(3,3);x=sysvar(L,m,x);A(1,1)=-x(3);A(2,1)=x(1);A(3,1)=b*x(1)-c*x(3)+x(1)*x(3);A(1,2)=x(4)-1.0;A(2,2)=x(5);A(3,2)=x(6);A(1,3)=x(7);A(2,3)=x(8);A(3,3)=x(9)-1.0;ppp=p+A\[p(2)-x(1); -x(2); p(3)-x(3)];L=ppp(1); x(1)=ppp(2); x(2)=0.0; x(3)=ppp(3);% Poincare map of statesfunction x=sysvar(L,m,x)x(4)=1.0; x(5)=0.0; x(6)=0.0; x(7)=0.0; x(8)=0.0; x(9)=1.0; 7 95h=L/m;for j=0:m-1t=j*h;[t, x]=RK(t,x,h);end% Runge Kutta methodfunction [t, z]=RK(t,z,h)f1=fn(t,z);f2=fn(t+h/2,z+h*f1/2);f3=fn(t+h/2,z+h*f2/2);f4=fn(t+h,z+h*f3);t=t+h;z=z+h*(f1+2*(f2+f3)+f4)/6;% Rosslers equationfunction dy=fn(t,y)global a b cdy=zeros(9,1);dy(1)=-y(2)-y(3);dy(2)=y(1)+a*y(2);dy(3)=b*y(1)-c*y(3)+y(1)*y(3);dy(4)=-y(5)-y(6);dy(5)=y(4)+a*y(5);dy(6)=(b+y(3))*y(4)+(y(1)-c)*y(6);dy(7)=-y(8)-y(9);dy(8)=y(7)+a*y(8);dy(9)=(b+y(3))*y(7)+(y(1)-c)*y(9);96 8Switched Dynamical Systems with MovingBorder8.1 Formation of the problem8.1.1 Definition of two dynamical systems and a border manifold BWe define two dynamical systems on Rn:dxdt= f (x, 1) , x Rn (8.1)anddydt= g (y, 2) , y Rn (8.2)where 1, 2 are the system parameters with appropriate dimensions.We denote the solution x(t) of Eq. (8.1) asx (t) = (t, t0, x0) (8.3)with the initial condition:x (t0) = (t0, t0, x0) = x0 (8.4)Similarly the solution y(t) of Eq. (8.2) is defined as:y (t) = (t, t0, y0) (8.5)with the initial condition:y (t0) = (t0, t0, y0) = y0 (8.6)By switching two flows of x(t) and y(t) appropriately, we would like to control the resultant flowon Rn. To this end, let us define a codimension-one manifold B by the pre-image of the function:B = {x Rn : (x, t) = 0} (8.7) 8 Switched Dynamical Systems with Moving Border 97where : Rn R R; (x, t) 7 (x, t) (8.8)The submanifold B will be temporally called the border manifold. We also assume that B is timevaring periodically with period L: (x, t+ L) = (x, t) , x B, t R (8.9)8.1.2 Switching action on the border BLet us construct a flow z(t) starting from the initial point P Rn at t = 0 and terminatingat the final point Q Rn at t = L. Switching action between two flows will be designed tooccur whenever our resultant flow meets the border manifold B. In [0, L] the number of switchinginstances is assumed as nL. As an example of z(t), let the initial flow is defined by Eq. (8.1) , andthe switching instances are0 < 1 < 2 < < nL < L (8.10)Now we can define the total trajectory z(t) as follows:1. At the first inteval 0 t 1 we have:z(t) = x (t) = (t, 0, P ) 0 t 1 (8.11)2. At the time t = 1 the point x(1) gives the initial point y (1) = P1 = x(1) for the flowy(t). The the solution enters into the seecond inteval 1 t 2.z(t) = y (t) = (t, 1, P1) 1 t 2 (8.12)3. Now we find the initial point P2 = y(2, 1, P1) for the flow x(t). The solution is obtainedbyz(t) = x (t) = (t, 2, P2) 2 t 3 (8.13)4. Similarly, we can define the trajectory z(t) consecutively at every interval k t k+1 k =3, , nL 1, and we obtain the solution in the final interval with the initial point PnL1.Note that in the final interval the flow is defined by Eq. (8.1) (or Eq. (8.2) ) , if nL iseven(or odd). Hence the final point Q is determined by the final flow at t = L.Since the border B moves periodically, we can define our resultant flow z(t) in every interval[0, L]. Hence z(t) can be defined in t [0, ).Under these circumstances we would like to investigate the dynamic behaviors of the total system. 8 Switched Dynamical Systems with Moving Border 98iRCLE vSW1SW2iRCLE vSW1SW2(a) (b) 8.1 Linear resonant circuits. (a) with a voltage source E. (b) without E.8.1.3 A simple exampleAs a simple example, we formulate the dynamics in the case of linear systems with a linear timevaring border.Let us consider linear circuits shown in Fig. 8.1 .The circuit equations are given by:Cdvdt= i vRLdidt= v + E(8.14)for the circuit (a) andCdvdt= i vRLdidt= v(8.15)for the circuit (b).For the border movement, we define (v, t) = av b (t) c = 0 (8.16)As concrete examples of b(t) we illustrate two cases:b (t) = cos2Lt (8.17)andb (t) =1Lt, t [0, L] (8.18)Note that the latter ramp function is discontinous at t = 0 or t = L. At this discontinous instancethe border v =1ab (t) +cajumps from v =1 + cato v =ca. 8 Switched Dynamical Systems with Moving Border 99This is the limiting case ( 0) of the functionb (t) =1L t, t [0, L ]1(t L) , t [L , L](8.19)This discontinuity causes the change of the number of switching nL, that is nL becomes nL 1in the discontinous case.8.2 The Poincare map and periodic solution8.2.1 Definition of the Poincare mapAs the border movement is periodical, we can use this periodicity for defining the Poincaremap under the suitable condition. To define the Poincare map, the vector field at t = 0 andt = L must be the same. Therefore, the number of switching nL has to be even when (x, t) iscontinous in time. Under this condition our total system is invariant at every instance of time:t = kL, k = 0, 1, . Hence we can define the Poincare mapping T as follows:T : Rn Rn; x (0) 7 x (L) (8.20)8.2.2 Periodic solutionthe case of continuous borderFor a moment we assume that the Poincare map is well defined and nL = 2. Assuming also theexistence of a periodic solution, we would like to calculate this periodic solution. See Fig. 8.2 .Following the solution illustrated in Fig. ??, we have the relations which give the periodicsolution:x1 = (1, 0, x0)x2 = (2, 1, x1)x0 = (L, 2, x2)(8.21)with the conditions: (x (1) , 1) = 0 (x (2) , 2) = 0(8.22)These relations contains 3n+ 2 unknown variables (x0, x1, x2, 1, 2). The number of relationsis also 3n+ 2. Hence we can solve these relations.the case of ramp borderIn particular, in the case of the ramp border Eq. (8.18) we can choose t = 0 at the discontinuousinstance. See Fig. 8.3 . In this case we have the periodicity conditions: 8 Switched Dynamical Systems with Moving Border 100 x, t)=0t 1t 1t=0 t=Lt 2t 2section(a) (b)x0x1x2 8.2 A simple periodic solution. (a) Intersection between the border and the solution .(b) Poincare section. x, t)=0t 1t 1t=0 t=Lsection(a) (b)x0 x1 8.3 A periodic solution with discontionous border. (a) Intersection between the borderand the solution . (b) Poincare section.x1 = (1, 0, x0)x0 = (L, 1, x1) ( (1, 0, x0) , 1) = 0(8.23)Or equivalently, we havex1 = (1, 0, x0)x0 = (L 1, 0, x1) (x1, 1) = 0(8.24) 8 Switched Dynamical Systems with Moving Border 101Assume for a moment the border function as: (x, t) = a x b (t) c = 0, t [0, L] (8.25)differential of periodic solutionAs an example of the differential of the Poincare map, we calculate the derivative of the followingrelations:x1 = (1 (x0) , 0, x0)x2 = (L 1 (x0) , 0, (1 (x0) , 0, x0)) ( (1 (x0) , 0, x0) , 1 (x0)) = 0(8.26)That is, we would like to obtain the derivative:x2x0.Firstly, from the first relation of Eq. (8.26)x1x0=t1x0+x0= f ()1x0+x0(8.27)Secondary, from the third relation of Eq. (8.26)x0=xx1x0+t1x0=x(f ()1x0+x0)+t1x0= 0(8.28)Hence we have1x0= 1xf () +txx0(8.29)Lastly, we obtainx2x0=t( 1x0)+x1x1x0= g () 1x0+x1(f ()1x0+x0)={g () + x1f ()}1x0+x1x0=x1x0+1xf () +t{g () x1f ()}xx0(8.30) 8 Switched Dynamical Systems with Moving Border 102In the 2-dimensional case with the border Eq. (8.25), our formula becomes as follows:x0=xx1x0+t1x0= ax1x0 b (1) 1x0= a(f1x0+x0) b (1) 1x0= {a1f1 + a2f2 b (1)} 1x0+ ax0= 0(8.31)Hence1x0=1b (1) a1f1 a2f2[a1 a2]1x101x202x102x20=1b (1) a1f1 a2f2[a11x10+ a22x10a11x20+ a22x20](8.32)x1x0={g () + x1f ()}1x0+x1x0=x1x0+g1 () + 1x11f1 () +1x21f2 ()g2 () + 2x11f1 () +2x21f2 () 1x0=x1x0+1b (1) a1f1 a2f2g1 () + 1x11f1 () +1x21f2 ()g2 () + 2x11f1 () +2x21f2 ()[a11x10+ a22x10a11x20+ a22x20](8.33)where we putx0 =x10x20 , x1 =x11x21 , f =f1f2 , g =g1g2x0=1x101x202x102x20 , x1=1x111x212x112x21(8.34) 8 Switched Dynamical Systems with Moving Border 1038.2.3 Numerical methodNow let us consider our numerical method for solving Eq. (8.24) under the border function Eq.(8.25). We have to find the roots of the following equations:F1 (x0, x1, 1) = x1 (1, 0, x0) = 0F2 (x0, x1, 1) = x0 (L 1, 0, x1) = 0F3 (x0, x1, 1) = a x1 b (1) c = 0(8.35)These equations can be solved by Newtons method. To this end we must calculate the Jacobianmatrix of Eq. (8.26) , which is obtained as:F1x0F1x1F11F2x0F2x1F21F3x0F3x1F31= x0(1, 0, x0) In f ( (1, 0, x0) , )In x1(L 1, 0, x1) g ( (L 1, 0, x1) , )0 a b (1)(8.36)8.3 Phase Portrait DemoPL8.4 DufFix8.5 DemoBF104 A1. Windows F1 dos F2 dos F3 BSD F1 Ctrl + Alt + Del F1 Windows2. MATLAB Windows MATLAB 3. MATLAB MATLAB exit4. Windows Windows 5. 105 B1. 2. B.1 3 1. (Turtle graphics: TG)2. G0(s) a a b c bG0 (s) =s+ a(s+ b) (s+ c)c G0(s) d G0(s) e G0(s) 3. B 106B.2 A4 : 3.5cm, : 2.5cm, : 3cm, : 2cm 1B.3 4 1231998 8 23B.3.1 1: 4 1. TI2. R(a)3. F (a)4. ST 3 1 120 3 B 107 % draw a triangle by step by step method TI; F(1); R(120); F(1); R(120); F(1); ST % draw a triangle by using for end statement TI; for i=1:3 F(1); R(120); end ST 3 3 2 function A(a) % draw a triangle for i=1:3 F(a); R(120); end % function B(a,b) % draw b number of triangles A(a); R(-120); c=a; for i=1:b R(125); c=0.956*c; F(c); R(120); c=0.864*c; F(c); R(120); F(c); B 108-1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.81 B.1 end 6 % draw triangles TI; for i=1:6 B(1,15); H([0; 0]); R(60*i); end H(a) function H(a) % return to the home position global X U Path U=[0;1]; Path=a; X=[X Path]; B.1 MATLAB B 109B.3.2 2 51 8 23 G0 (s) =s+ 51(s+ 8) (s+ 23)D (s) = (s+ 8) (s+ 23) = s + 31s+ 184 % step response of the open transfer function num=[0 1 51]; den=[1 31 184]; step(num, den) B.20 0.1 0.2 0.3 0.4 0.5 0.6 0.700.050.10.150.20.250.3Time (secs)Amplitude B.2 % bode diagram num=[0 1 51]; den=[1 31 184]; bode(num, den) B.3 B 11010-1 100 101 102 103-60-40-200Frequency (rad/sec)Gain dB10-1 100 101 102 103-30-60-90-120-1500Frequency (rad/sec)Phase deg B.3 % bode diagram num=[0 1 51]; den=[1 31 184]; nyquist(num, den) B.4B.3.3 3 B 111-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3-0.2-0.15-0.1-0.0500.050.10.150.2Real AxisImag Axis B.4

Recommended

View more >