참고도서 Matlab Tutorial 참고도서 : Mastering Matlab 7, 4, Hanselman and Littlefield, Prentice Hall. Mastering AutoCAD 5 and AutoCAD LT 5, 4, Omura, Sybex. Autocad,, 김영숙, 세진사. Prof. Youngseuk Keehm ( 김영석 ) Tel: 4-85-857 Office: 자연과학관 3 Email: keehm@kongju.ac.kr 조교 : 이민희 (45 호, -636-4383) nd Semester, 5 What you can do in Matlab? MATLAB = MATrix LABoratory. Elementary function, data analysis, D and 3D graphics, interpolation, grid. Signals, Fourier transforms, special functions 3. Least-squares, matrix function, integration, differential functions 4. Statistics, stochastic simulations 5. Image processing, Movie 6. CAD 개요, 좌표계, 기본도형및응용 7. 개체스냅및편집, 문자와치수기입 8. 레이어의이해와이용 9. 도면배치및출력. 면구조 (Surface) 와렌더링 (Rendering). 3차원면의재질및광원효과 Matlab - Introduction. Starting Matlab who(s), ls, pwd, cd, diary. Datatypes 3. Creating variables 4. Vector/Matrix Variables and Accessing 5. Display format 6. plotting 7. Vector/Matrix operators 8. Max, Min, Mean, Sum 9. Transpose, flip. ones, zeros, eye, cascading, extracting. load data. cell structure 3. Basic matrix calculations Basics of Matlab. Review of Matlab basics 컴퓨터구조, 파일시스템, OS who(s), ls, pwd, cd, diary. 변수의선언 스칼라, 벡터, 행렬 help datatypes 3. 변수들에서필요한요소의추출 4. Vector/Matrix 연산자들 + - *.* /./ \.\ ^.^ 5. Max, Min, Mean, Sum, Std Basics of Matlab cont. 6. Transpose, flip 7. ones, zeros, eye, cascading, extracting 8. 매트랩데이터의로드와저장 9. 매트랩기본함수들 (elfun). Loop. If, else, end 문의구조
파일시스템과메모리 변수의선언. 컴퓨터구조, 파일시스템 ( 하드디스크 ), Ram - 데이터를저장하거나읽어올때는하드디스크를 access 하는데매트랩메인윈도우의 current directory 에서확인한다. 매트랩을실행하면자동으로 C:\Matlab7.\Work라는디렉토리에서시작한다. 이것을변경할경우 cd 라는명령어를사용하고, 현재자기가어디있는지확인하고싶은경우는 pwd 를쓴다. 현재디렉토리에어떤파일이있는가를보고싶은경우는 ls 나 dir 을사용한다. - 현재쓰고있는변수들은메모리상에위치하는데이를매트랩에선 workspace 라고한다. 지금어떠한변수를가지고있는지확인할때는 whos 라는명령어를쓴다.. 변수의선언 스칼라, 벡터, 행렬, integer, single-precision, double-precision ( help datatypes 로가능한변수의타입을확인할수있다.) 어떠한데이터타입인지지정하지않고변수를만드는경우는자동으로 double precision으로생성된다. (matlab에선 double type). 변수의선언 ( 계속 ) - 변수를선언하는방법은여러가지가있는데 () 파일에서읽어오는방법, () 계산식으로만드는방법 (3) 초기값을주고일단정의하는방법등이있다. () matlab data file (*.ma 이있는경우 : load mydata.mat 일반 text (ascii) 파일인경우 : load, textread, dlmread, fscanf () 계산식으로만드는방법 ( 예 ) >> x=linspace(, *pi, ); >> y=sin(x); (3) 초기값으로정의하는방법 (ones, zeros) ( 예 ) >> x=zeros(,) % x 행렬을만들고 으로초기화 >> y=*ones(4,) % 4x 행렬을만들고 으로초기화 >> z=eye(5) % 5x5 의단위행렬 (I) 생성 변수에서필요한부분추출 변수에서필요한부분추출 3. 변수의추출 () 필요한범위를 : 를이용해추출하는방법 >> a = rand(,); % x random number 행렬 >> b = a(:5, 4:6); % 행에서 4행까지, 4열에서 6열까지추출 >> c = a([ 4 5], [3 8 9]) %,4,5행과 3,8,9열추출 >> d = a(:, [9 8 6]) % 9,8,6번째열의모든원소를이순서로배열 () 논리식으로추출하는방법 >> a = [ 3-8 9 - -5]; >> b = a( a<) = [ -8 - -5] % b는 3개의요소를가지는 vector >> b = a( find(a<) ) % 위와같은결과 >> ind = a< = [ ]; % 원래 vector와같은사이즈, boolean >> ind = find(a<) = [3 6 7]; % 실제인덱스를리턴 3. 변수의추출 () 논리식으로추출하는방법 ( 계속 ) >> a = [ 3-8 9 - -5]; >> b = a( a> & a<5 ) = [ 3 ] % 논리연산 AND &: AND, : OR, ~=: NOT EQUAL; ==: EQUAL; >, <, >=, <= >> c = [ -5 4 9-3 - 4 3 ; 3-4 -4 4-8 9 -]; >> d = c( [c(,:)> & c(,:)<], :) = [ 9 ; -4 -]; c의 행이 보다크고 행이 보다작은모든열의집합 * 행렬에서논리식으로추출할경우특히주의를하여야한다. 변수에서필요한부분추출 3. 변수의추출논리식으로추출하는방법 ( 예제문제 ) >> a = [ 3 4 5 6 7 8 9 ;..4 -...4.5.6.3. -.]; () 행은시간이고 행이해수의상대적높이라고할때, 해수의상대적높이가음수인경우를제외한데이터 b를구하라. >> b = a( [a(,:)>=], :); * 행과열의논리식을따로써주는것에주의 () 시간이 3과 8사이에있으면서, 해수의높이가. 보다높은경우인데이터 c를구하라 >> b = a( [a(,:)>=3 & a(,:)<=8 & a(,:)>.], :); 행렬의연산 4. Vector/Matrix 연산자들 + - *.* /./ \.\ ^.^ 연산자앞에마침표가있는경우에는각요소들간의연산으로정의가된다. ( 예 ) 3 5 A =, B = 4 3 3 4 5 A* B =, A.* B = 4 ( 참고 ).. A/ B = A* B.4..3.9 A \ B = A * B..7 ( 질문 ) 왜나눗셈 operator 가두개일까? (/ 와 \)
간단한통계함수 5. 통계 operator : min(), max(), sum(), mean(), median(), var(), std() : bar(), bar3(), hist(), hist3(), boxplot() : rand(), randn(), unidrnd(), unifrnd() ( 문제 ) A 가행렬인경우 max(a) 는? ( 문제 ) A 가행렬인경우한번에평균값을구하는방법? ( 문제 ) mean 과 median 의차이는? ( 문제 ) bar 와 hist 의차이는? 데이터의로드와저장 6. 데이터 ( 파일 ) 의입출력 (a) 현재메모리 (workspace) 에있는변수, 데이터의저장 : save() 를이용한다. ( 예 ) >> save mydata.mat data data data3 : 이경우에는 mydata.mat 는 binary file 이된다. >> save mydata.txt data4 ascii : 이경우는 acsii 파일즉아무뷰어로나볼수있는텍스트파일이된다. (b) 데이터파일을읽어오는경우 : load() 를이용한다. >> load mydata.mat % matlab format file >> load mydata.txt % ascii file ( 참고 ) ascii file 을읽을경우데이터포맷의열갯수가모든행에서일치해야만한다. ( 다른예 ) textread, dlmread, fscanf (help 를사용하여그용도를익힐것 ) 예제문제 현재까지배운내용을다음예제를통해서익힌다. 주어진 data (hwdata.ma 를로드해서다음을실습해본다. (a) 데이터를읽어서포함되어있는변수와크기를알아보라. (b) 변수중하나를 ascii format 으로데이타파일로저장하고, notepad 로그내용을확인해보라. (c) 변수중하나를 plot, scatter, hist 를이용해서그래프를그려보라. (d) 변수중에서행렬로되어있는것을골라일정깊이에있는것만을선택적으로추출해보라. (e) 추출된데이터의기본적인통계학적값들을구해보라. D & 3D graphics 매트랩이공학이나수치해석에서강력한툴이된이유중의하나가내장되어있는다양하고매우유용한그래픽함수때문이다. 실제로계산이나데이터처리를한이후에는그것을어떻게 visualization 할것인가는이제수치해석자체만큼이나중요한문제가되었다. 이번장에서는기본적으로매트랩이제공하는여러가지 graphic function 들에대해알아보고실제예제와과제을통해자기것으로만드는것을목적으로한다. 매트랩이제공하는방대한양의그래픽함수들에대해세세히다루기는불가능하므로각자개인이 graphics 에어떠한함수와유틸리티가있는지알아보기바람. >> help graphics ( 또는 graphd, graph3d, specgraph) Cross plots Graphics Handler Cross plot 이란 차원에 x 와 y 값을가지는변수들을선, 점, 도형등으로표시하는가장기본적인것을말한다. 이범주에드는가장대표적인함수들이 plot(), scatter(), plotyy(), ezplot(), semilogx(), semilogy(), loglo) 등이있다. ( 예 ) >> plot(x, y, r ) % line plot with red line >> scatter(x, y, s, c) % 점으로표현된그림, s는각점의크기, c는칼라 >> plotyy(x,y, x,y) % 두개의그래프를한그림에표시 >> ezplot( sin(x), [ 6*pi]) % 문자로표현된함수의그림그리기 >> semilogy(x, y, ro ) % 데이터를빨간원으로표시, y축 log scale * 각함수의 help 를살펴보면예제와옵션이자세히나와있다. 이러한그래픽들은아주다양한옵션들 ( 예를들어선의굵기, 좌표축의최소, 최대값, 각축에있는tic mark를어디다찍을지, 각축의라벨의폰트와사이즈등등 ) 을가지고있는데이를함수자체에서인자로모두넣기란불가능하다. 그래서 Figure나 Graph를그리는모든함수는그것들의 handler ( 제어자 ) 를넘겨주는데다음의예를살펴보자. >> h = figure; >> h = plot(x, y, g-- ) 여기서 h과 h는각각그림자체와그속에있는graph를제어할수있는제어자가된다. 이제어자가어떠한속성들을갖고있는지보려면, >> get(h) 이제어자를이용해그림의속성중하나를바꾸려면 ( 예를들어선두께 ), >> set(h, linewih, 3) 와같이쓰면된다. >> axis % 현재축의최대, 최소값표시또는셋팅 3
Graphics Handler Graphics Handler >> get(h) % Figure Alphamap = [ ( by 64) double array] BackingStore = on CloseRequestFcn = closereq Color = [.8.8.8] Colormap = [ (64 by 3) double array] CurrentAxes = [34.] CurrentCharacter = h CurrentObject = [] CurrentPoint = [ ] DockControls = on DoubleBuffer = on Position = [4 58 56 4] Renderer = painters >> get(h) % Plot Color: [ ] EraseMode: 'normal' LineStyle: '-' LineWih:.5 Marker: 'none' MarkerSize: 6 MarkerEdgeColor: 'auto' MarkerFaceColor: 'none' XData: [ 3 4 5 6 7 8 9 ] YData: [x double] ZData: [x double] BeingDeleted: 'off' ButtonDownFcn: [] Children: [x double] Graphics Handler Subplot >> get(gca) % Axes Color = [ ] FontName = Helvetica FontSize = [] FontUnits = points FontWeight = normal LineStyleOrder = - LineWih = [.5] Position = [.3..775.85] TickLength = [..5] TickDir = in XTick = [ ( by ) double array] XTickLabel = XTickLabelMode = auto YTick = [ ( by 9) double array] YTickLabel =. Children = [35.] 하나의그림판에여러개의그래프를표시해야하는경우가종종생긴다. # of data Value 8 6 4 -.5 - -.5 - -.5.5.5 Value - - -3..4.6.8 Time 이러한경우 subplot() 을사용한다. >> subplot(,, ) >> hist(value, ) >> subplot(,, ) >> plot(time, value,. ) Surface Plot surf() Surface Plot surf() 면구조를그리는경우의가장쉬운예는 Topographic map이다. 이경우에차원 (x,y) 의지점마다높이가주어진다. 이러한경우는 surf() 가가장대표적으로쓰이는함수이다. z 라는변수가 x 의크기를가지고각각의위치는높이를나타낸다고하자. 이경우각각의좌표는등간격을가진다. (dx=dy=cons >> z = peaks(); >> surf(z) >> shading interp >> camlight >> lighting phong >> print djpeg99 a.jpg 앞의예제와같은면구조를나타내는또하나의범용적인방법은등고선으로나타내는것이다. 다음의예제는 contour() 을이용해등고선도를그린것이다. 등고선도의높이를나타내기위해서 colorbar() 를이용하여그림의우측에등고선색의높이를표시하고있다. >> z = peaks(); >> contour(z, [-:.5:] ) >> colorbar 9 8 7 6 5 4 3 6 4 - -4 3 4 5 6 7 8 9-6 4
Surface Plot 비등간격데이타 비등간격데이타 - 예제 앞의예제에서는주어진자료가등간격으로 (Grid data) 주어져있기때문에그값을그냥이용해서면구조나등고선도를그릴수있었다. 하지만많은경우데이타가등간격이아니므로, 이경우는등간격의데이타를생성하는작없이필요하다. Griddata() 또는 interp() 이용. < 주어진데이타 > 이경우에는주어진데이타가 x, y, z 세개일것이다. 그림에서보듯이 x, y가등간격이아니므로먼저등간격의 x,y을만든다. 9 >> [x,y]=meshgrid(::, ::); 8 그다음엔주어진자료들로부터 (x,y) 에서 7 의예측값들을구한다. 6 >> z = griddata(x, y, z, x, y); 5 또는 4 >> z = interp(x, y, z, x, y); 3 Original Data 비등간격, incomplete data griddata 이용후의결과 3 4 5 6 7 8 9 그림화일과같은 차원구조 3 차원데이타 Graphics 앞에서우리는지각과같은면의높이가주어진상태에서 relief map 이나등고선도를그리는방법에대해서알아보았다. 우리가잘아는그림파일들도 (bmp, jpg, tiff) 이와똑같이 차원구조를가진다. 그림파일을표시하는함수는 imagesc(), pcolor() 와같은것이있고 차원 data 에는언제나사용가능하다. 그림파일을특수하게다루는함수에는 imread(), imshow() 등이있다. 4 6 8 4 6 8 앞에까지다룬예제는실제는 차원 data 이다. 여기서는실제 3 차원자료, 즉 3 차원 grid 의각점마다데이타가주어진경우를다루어보자. Data 는 3 차원에서구슬을채워넣은공극구조 (random pack of spheres) 를염수로포화시키고전류를흘려서전류장을구한결과이다. >> load curr.mat >> whos curr Name Size curr xx >> subplot(,,) >> imagesc(z) >> subplot(,,) >> pcolor(z) 8 6 4 >> subplot() >> slice(curr, [],[],[]) >> shading flat >> axis tight >> axis image >> colorbar 4 6 8 3 차원데이타 Graphics NaN 앞의그림은가장자리면의전류를표시하기때문에실제어떻게전류가흘렀는지보기가힘들다. 다음과같이하면보다직관적으로이해하기쉬운그림을그릴수있다. >> slice(curr, [ ],[ ],[4 4 7]) >> shading flat >> axis tight >> colorbar NaN은특별한데이터타입으로Not-A-Number의약자이다. 실제로정규적인데이타에서빠진데이터나에러가난데이타를-999나다른의미없는숫자로대체하는데이와비슷한역할을한다. 이러한 NaN을포함한데이타를위해특별한함수들이있다. nanmean(), nanmedian(), nanstd(), nanvar(), nanmin(), nanmax(), nansum(), isnan() 5
Linkaxes() Graphic Handler revisited >> ax() = subplot(,,); >> plot(vshale,depth), axis ij >> ax() = subplot(,,); >> help linkaxes를이용하여 >> plot(porosity,depth), axis ij y축만을링크하는방법을알아보라. >> linkaxes(ax) >> h=plot(vshale); >> axis ij >> set(gca, fontname, tahoma,fontsize,) >> xlabel('vsh') >> ylabel('depth (f') x 4..4.6.8 x 4..4 x 4..4 3 3 Depth (f.6.6 4 4..8.8 5 5.4 6 6.. 7 7.6.4.4.6.6 8 8.8..4.6.8 Vsh.8.8.5.5..4.6.8..4.6.8 Graphic Handler revisited Graphic Handler revisited >> h=plot(vshale,depth); >> axis ij x 4. >> h=plot(vshale); >> axis ij x 4. >> set(gca, fontname, tahoma, fontsize,).4 >> set(gca, fontname, tahoma,fontsize,).4 >> xlabel('vsh').6 >> xlabel('vsh').6 >> ylabel('depth (f') >> set(h, color, r ) Depth (f.8 >> ylabel('depth (f') >> set(h, color, r ) >> set(gca, color, b ) Depth (f.8...4.4.6.6.8..4.6.8.8..4.6.8 Vsh Vsh Graphic Handler revisited 많이쓰이는 handler s properties >> h=plot(vshale); >> axis ij >> set(gca, fontname, tahoma,fontsize,) >> xlabel('vsh') >> ylabel('depth (f') >> set(h, color, r ) >> set(gca, color, b ) >> set(gcf, color, g ) Depth (f x 4..4.6.8. >> set(gca, position, ) % 그림판안에서 plot의상대적위치 >> set(gca, fontname, times, fontsize, ) % x,y축의레이블폰트와사이즈 >> set(gca, xtick, [ ]) % x축의 tic 마크의빈도와위치 >> set(gca, xticklabel, { One, Two, Three }) % tic label이특수한문자일때 >> set(gcf, renderer, zbuffer ) % 면구조렌더링방법 (painter, zbuffer, opengl).4.6 나머지는각자알아볼것..8..4.6.8 Vsh 6
그림창의메뉴들 Plot 과같이많이쓰이는함수들 File new, open, close, save, save as, generate m-file, preferences, export setup, page setup, print setup, print preview, print Edit copy figure, copy options, figure properties, axes properties View figure toolbars, camera toolbars, plot edit toolbars, property editor Insert x,y,z label, title, legend, line, arrow, etc. Tools edit plot, zoom, rotate, axes properties, basic fitting, data statistics title() xlabel() ylabel() text() gtext() grid on/off hold on/off axis() subplot() axes() ginput() view() Matlab Function: m-file Matlab Function: m-file Fortran 의 subroutine 이나 function, 또는 C/C++ 의 method function 과같이동일한작업수행이요구되는일에는 function 을만드는것이능률적이다. 매트랩에서우리가흔히쓰는많은명령어들은이런한 function 들이며대부분이.m 이라는확장자를가진텍스트 (ASCII) 파일들이다. 이러한파일을생성, 수정하는데에는 edit 명령어가쓰인다. 예를들어서평균을구하는 mean 함수의경우를살펴보자. >> edit mean 을치면 edit window 가뜨고 mean() 함수의내용이창에표시된다. function y = mean(x,dim) %MEAN Average or mean value. % For vectors, MEAN(X) is the mean value of the elements in X. For % matrices, MEAN(X) is a row vector containing the mean value of % each column. For N-D arrays, MEAN(X) is the mean value of the % elements along the first non-singleton dimension of X. % % MEAN(X,DIM) takes the mean along the dimension DIM of X. % % Example: If X = [ % 3 4 5] % % then mean(x,) is [.5.5 3.5] and mean(x,) is [ % 4] % % Class support for input X: % float: double, single Matlab Function: m-file 함수의선언 % % See also MEDIAN, STD, MIN, MAX, COV. % Copyright 984-4 The MathWorks, Inc. % $Revision: 5.7.4. $ $Date: 4/3/9 6:6:6 $ if nargin==, % Determine which dimension SUM will use dim = min(find(size(x)~=)); if isempty(dim), dim = ; end y = sum(x)/size(x,dim); else y = sum(x,dim)/size(x,dim); end 함수의맨처음은항상다음과같은형태를가진다. function y = mean(x,dim) 이파일이함수임을나타냄 두개의 input 변수를가짐 함수의이름이 mean 계산결과가 y 라는변수로출력 예를들어두개의데이타 x,y를받아서특별한그림을그리는함수를우리가만든다고하자. 그러면그함수는다음과같은형태를지닐수있다. function myplot(x,y) 이경우는output이화면에그림으로출력되기때문에따로변수 ( 계산결과 ) 를돌려줄필요가없어서 y가생략되어있다. 7
함수의설명문 함수의설명문 함수의선언부바로다음에빈줄이없이 % 기호로시작되는설명문이등장한다. function y = mean(x,dim) %MEAN Average or mean value. % For vectors, MEAN(X) is the mean value of the elements in X. For % matrices, MEAN(X) is a row vector containing the mean value of % each column. For N-D arrays, MEAN(X) is the mean value of the % elements along the first non-singleton dimension of X. % 이부분에는대개이함수의설명과그사용법이들어간다. 앞에서배운 help를치면바로이부분이여러분의command window에표시된다. 다시말하면함수의 help 에해당되는부분이다. % 로시작되는것은설명문이라고해석된다. (Fortran에서 C로시작하는 comment 문과같은역할 ). 하지만함수선언부분과이설명문사이에빈줄이없다는사실에주의하라. 만약빈줄이있으면이부분은help를이용할때표시되는않는다. 설명문을따라가다보면빈줄이있고설명문이계속된다. 따라서 help 를이용했을때, 아랫부분은표시되지않는다. % % Class support for input X: % float: double, single % % See also MEDIAN, STD, MIN, MAX, COV. % Copyright 984-4 The MathWorks, Inc. % $Revision: 5.7.4. $ $Date: 4/3/9 6:6:6 $ 여기에서 Copyright 와버젼등의관계를숨어있는설명문에포함하고있다. 함수의내용 (implementation) 함수의내용 (implementation) dim = min(find(size(x)~=)); 처음이문을접하면상당히복잡해보인다. 하나하나씩살펴보면, size(x) : size() 는input 변수의사이즈를알려주는함수이다. 만약 A가 3x4x의 3차원행렬이면size(A) 는 [3 4 ] 를 return한다. find( size(x) ~= ) : find는 size(x) 가 이아닌경우의그첨자 (index) 를 return 한다. 현재 A의경우는모두가크기가 이아니므로 [ 3] 을리턴한다. min(find(size(x)~=)): 찾아진첨자중제일작은것을리턴한다. 여기서는 이제일작다. 그러므로 방향이최종계산값이다. A=[ 3 5] ( 행벡터 ) 의경우, 위의표현은어느방향을리턴하는가? size(a) = [ 4] find( size(a) ~= ) = [] % 방향은사이즈가 이므로 min(find(size(a)~=)) = % 값이하나밖에없으므로 를리턴하고최종결과는 방향이된다. dim = min(find(size(x)~=)); 이제왜이렇게복잡하게쓰는지를어는정도이해할수있다. 만약 3차원에있는어떤벡터가 (xx) 의사이즈를가지고있고, 평균값을구하는함수가항상방향을명시해야한다면이의경우는 >> mean(a,3) 이라고항상써야하고이는불편하기그지없다. 따라서처음으로사이즈가 이아닌방향을찾아서평균을구하는것이논리적으로맞다는것이다. 만약 A 가 (xx5) 의크기를갖는행렬이라면 >> mean(a) 의결과값의사이즈는? ( x x ) 함수의내용 (implementation) 함수의내용 (implementation) function y = mean(x,dim) if nargin==, 우리가함수을선언할때, 이함수는두개의 input변수를필요로한다. 하지만함수안에서한개만넣어도되는경우를처리하고있으면한개만넣어도된다. 그것이첫번째에나오는 nargin== 의부분이다. 이말은input 변수가하나뿐일때, 즉 x 뿐일때어떻게계산을수행하는가를묻는말이다. 이경우를잠시생각해보자. Case : x=[,, 3, 5] ( 행벡터 ). 이경우user는이네개값의평균을원할것이다. Case : x=[, 5, 7, 9] ( 열벡터 ). 이경우도이네개값의평균을원할것이다. Case 3: x=[ ; 4 5] ( 행렬 ). 이경우평균을구하는방법은두가지가있다는방법을우리는배웠다. ( 방향, 방향 ). 그래서방향을정하지않은경우는 방향이라고가정한다. 물론, if 문을써서각각의경우에방향을정해주면되지만이것을간단히하면, if isempty(dim), dim = ; end 위에서우리는평균을구하는방향을정했는데, 또방향을구하는하나의명령문이들어가있다. 이경우는어떤경우인가? 만약 user 가 input 변수로사이즈가 인변수를 input 으로넣었다고하자. 사이즈가 인변수가과연어떻게존재하는가? >> A=ones(,) 과같은경우가해당되는데이와같은문을쓰는 user 는없다. 하지만우리가어떤계산을수행하였는데결과가없을때매트랩은공집합, 즉사이즈가 인변수를리턴한다. 만약우리가그결과를평균함수에대입하는경우가가장빈번히발생하는경우이다. 따라서이경우도처리를해주어야하는데여기서는 isempty() 라는함수를쓴다. 만약크기가 인변수의경우이함수값은 이된다. 그경우위의 dim = min(find(size(x)~=)); 는 의방향을리턴하므로우리는이경우를 의방향으로명시할필요가있다. 8
함수의내용 (implementation) If-else-end 문과 Loop 문 if nargin == y = sum(x)/size(x,dim); else y = sum(x,dim)/size(x,dim); end 그다음으로넘어가보자. 이제우리는 user 가방향을지정하지않은경우방향을모두구했다. 따라서평균값은 y = sum(x)/size(x,dim); 와같이모든값을더해서그갯수로나누어주면된다. User 가방향을정한경우는그방향으로합을구하고, 그방향으로의갯수 ( 사이즈 ) 로나누어주면이함수의기능은끝이난다. 이제만약크기가 인값을 input 변수로넣었을경우를다시생각해보면, 그방향은앞에서 로잡았다. 그런데이값들의합도 이고갯수도 이므로결국 / 의꼴이되어서 error 가생긴다. 이러한 error 도이함수가처리해야하는가는함수를만든사람의생각에따라다르다. 현재매트랩함수 mean() 은이경우에 divided-by-zero error 를낸다. 실제로하나의함수을만들어보기전에알고리듬의기본이되는두가지구조 조건문 (if) 과루프반복에대해서간단히배워보자. 조건문 (if-else-end) 계산식이조건등에따라틀려질경우 if 문을쓰게되는데, 다른언어에서와같이그구조는비슷한다. 하지만 matlab 의경우는반드시 end 로그조건식을끝내야한다. if(a>) if(a>=) a=a*; else a=a*5; end elseif(a==) a=; else a=-5*a; end If-else-end 문과 Loop 문 함수만들기 예제 ( 벡터의내적 ) 다음은반복계산을하는경우를예로들어보자. Fortran 에서는 do-loop 를많이쓰는데매트랩에서는 for-loop 가표준이다. 다음의예는 개의 x 값에대해서그값은로그값의사인값을구하는 loop 의예제이다. 여기에서는실제로함수를만들어보기로하자. 여러분들이잘아는벡터의내적을구하는프로그램을짜보기로하자. 실제함수를만들기전에우리가알아야할것들을나열해보면, y=zeros(,); for i=: y(i) = sin( log (x(i)) ); end 많은경우에매트랩의벡터나행렬 operator 를이용해서 loop 를없앨수가있는데이방법이계산시간이적게걸리므로장려되는방법이다. y=sin( log (x) ); 는앞의루프를이용한식과같은결과를제공한다. 그리고 while() 문이라는것도있는데 while( 조건 )-end 사이에있는것을조건이참인동안반복실행한다. () Input 변수는몇개이어야하는가? () 벡터의내적계산이정상적으로수행되기위해선어떤조건이필요한가? (3) 리턴해야할변수가있는가? 있다면그사이즈는얼마인가? (4) 특수하게다루어야할경우가어떤것들이있을까? 벡터의내적을구하는함수 벡터의내적을구하는함수 여기에서는실제로함수를만들어보기로하자. 여러분들이잘아는벡터의내적을구하는프로그램을짜보기로하자. 실제함수를만들기전에우리가알아야할것들을나열해보면, () Input 변수는몇개이어야하는가? () 벡터의내적계산이정상적으로수행되기위해선어떤조건이필요한가? (3) 리턴해야할변수가있는가? 있다면그사이즈는얼마인가? (4) 특수하게다루어야할경우가어떤것들이있을까? function c=dot(a, b) % Inner product of two vectors % if( length(find(size(a)~=)) > length(find(size(b)~=)) > ) disp( Input parameters should be vectors ); end if( length(a) ~= length(b) ) disp( Two vectors should have the same size ); end c=. for i=:length(a) c=c+a(i)*b(i); 또는 c= a(:) * b(:); end 9
y Review on Functions/Parameters covered m-file revisited edit(), nargin, find(), isempty(), if-elseif-else-end, for-end, while length(), disp(), sqrt(), input() 위의함수중알지못하는것은 help를이용하여익히기바람. 그리고조건문중에서 switch 문의구조도익히기바람. 앞에서배웠던 m-file의경우는함수를정의하고그알고리듬을프로그래밍하는방법이다. m-file 은꼭함수 (function) 일필요는없다. 단순히명령어들의나열로도실행이가능하다 (script file). < 예 > Sine Function.5 % sine 함수를 plot 하는 script t=linspace(,8*pi,); plot(t, sin(, r ); axis([ 8*pi -.5.5]); grid on xlabel( Time ), ylabel( sin(time) ) title( Sine Function ) sin(time).5 -.5 - -.5 5 5 5 Time m-file revisited 수치알고리듬의 implementation 하지만앞의경우에있어선입력변수가없으므로여러가지다른경우에반복계산을하기가매우불편하다. 만약우리가시간의시작값과끝값을바꾸고싶다면다음과같이함수를만들면된다. function sineplot(stime, etime) t=linspace(stime, etime, ); plot(t, sin(, 'r'); axis([stime, etime, -.5,.5]); grid on xlabel('time'), ylabel('sin(time)'), title('sine Function') -.5 >> sineplot(, 4*pi) 4 6 8 sin(time).5.5 -.5 - Sine Function Time 수식으로주어지는알고리듬을계산하는예제를몇가지살펴보자. ( 예제 ) 기울기가 a 이고 y절편이 b인직선을 [x x] 구간에서도시하는매트랩프로그램을작성하라. ezplot이나 symbolic math를이용해서analytic하게그림을그리는경우도있겠지만, 많은경우수치해석에서는많은점들을이어그려서연속된함수의그래프를완성한다. 이점들은 (x, y) 와같이주어진다는사실은우리는알고있다. 독립변수인 x를어떤구간에서정의하고 y값은직선의방정식을이용해서구해보자. 위의직선의방정식은y = ax + b 이다. x = linspace(x, x, ); y = a*x + b; 수치알고리듬의 implementation 수치알고리듬의 implementation 그리고이좌표들을 plot으로도시하면된다. plot(x, y) 이를정리해서함수로만들면, 리턴하는변수는필요없고입력변수는 4개 ( 기울기, 절편, x 시작값, x 끝값 ) 가된다. function drawline(a, b, x, x) x=linspace(x, x, ); y=a*x + b; plot(x, y); grid on xlabel( x ), ylabel( y ) 5 5-5 - -5 ( 예제 ) 반지름이 r 이고그중심이 (x, y) 인원을도시하는매트랩프로그램을작성하라. 일단우리는원의방정식을알고있다. ( x x) + ( y y) = r 하지만이경우x와 y의범위를구하기가어려우므로, 극좌표계를쓰는것이훨씬효과적이다. 중심이 (, ) 이고반지름이 r인원의방정식을각 (θ) 으로표현하면, x = r cos( θ ) y = r sin( θ ) θ r ( r cosθ, r sinθ ) - -5 >> drawline(, -, -, ) - -8-6 -4-4 6 8 x
y y 수치알고리듬의 implementation 수치알고리듬의 implementation 중심을 (x, y) 으로이동하는 앞의내용을함수로정리하면, 것은단지평행이동이기때문에, x = r cos( θ ) + x y = r sin( θ ) + y 그러면, 우리는원의 x와 y의좌표를각으로표현할수있다. ( r cosθ + x, r sinθ + y) r θ (x, y) function plotcircle(r, x, y) theta = linspace(, *pi, ); x = r * cos( theta) + x; y = r * sin( theta) + y; plot(x, y), grid on 9 8 7 6 5 theta = linspace(, *pi, ); x = r * cos(theta) + x; y = r * sin(theta) + y; axis image xlabel('x'), ylabel('y') 4 3 >> plotcircle(6,, 5) - 4 6 8 x 수치알고리듬의 implementation 복잡한수식을가진함수 - Viscoelasticity 앞의함수의경우원의경계가Axis의경계와맞물려보기가좋지않으므로 axis의최대, 최소값을조정할필요가있다. function plotcircle(r, x, y) theta = linspace(, *pi, ); x = r * cos( theta) + x; y = r * sin( theta) + y; plot(x, y), grid on axis image xlabel('x'), ylabel('y') axis( [min(x)-, max(x)+, min(y)-, max(y)+] ) 8 6 4-4 - 4 6 8 x 지각의구성물 ( 암석, 토양 ) 이완전선형탄성체이면매질에주어진응력과그응력으로인한변형률은선형으로연관된다. 그경우우리가잘아는탄성파의속도는 M μ Vp =, Vs =, ρ ρ 가된다 (M 과 μ 는암석의탄성계수이다 ). 하지만실제암석은응력을작용했을때일부는열에너지로소멸된다. 이는마치암석이유체처럼점도를가진것처럼행동하는것으로해석될수있다. Elasticity ( 탄성이론 ) 에서이러한효과를첨가한것이 Viscoelasticity( 점성탄성이론 ) 이며, 여러가지모델이있지만표준선형고체 (Standard linear solid) 모델을많이쓴다. E E >> plotcircle(6,, 5) η 복잡한수식을가진함수 - Viscoelasticity 복잡한수식을가진함수 - Viscoelasticity 앞의그림과같이스프링과데시팟 (dashpot 출입문등에문이갑자기닫히지않게하기위해붙여놓는장치 ) 의조합으로이루어진모델은응력이어떻게작용하는가에따라그반응이달라진다. 가령우리가이모델을빠른속도로당긴다고하면데시팟은잘움직이지않아서탄성이전혀없는벽처럼반응할것이다. 반면에우리가매우느린속도로이모델을당길경우데시팟은저항이거의없이움직이기때문에이계는두개의스프링만있는모델처럼행동할것이다. 이말을다시풀이하면우리가매질을어떻게진동시키는가에따라서그반응이달라진다, 즉매질의탄성파의전파속도가달라진다는말이된다. 일반적으로매질을빠르게진동시키는경우, 즉탄성파의주파수가높은경우탄성파의속도는빨라진다. 이렇게주파수가달라짐에따라탄성파의속도가변하는것을 dispersion( 파의분산 ) 이라고한다. 그리고또하나의현상으로이러한매질을진행할때파는그에너지를잃어버린다. 이를파의감쇄 (attenuation) 이라고한다. 앞의결과를자세히이야기하기는너무복잡하므로, 다음요약된그림으로보면 x 축은주파수이고두개의그래프는탄성계수와감쇄계수 (/Q) 이다. 탄성파의속도는탄성계수에비례하므로주파수가커지면속도가커진다. 한편감쇄는공명주파수주위에서제일크다. 공명주파수는암석의속도변화가가장큰곳이다.
복잡한수식을가진함수 - Viscoelasticity 복잡한수식을가진함수 - Viscoelasticity 앞의내용을적당한근사를이용해서수식으로풀어보면다음과같은식을갖는다. 참고로감쇄가일어나는경우탄성계수는복소수가되는데이것의실수성분을 Re[], 허수성분을 Im[] 로표시한다. 앞의수식을 M 가, M 가 4, 밀도가 이라고하면다음과같은그림이나온다. 이제우리는이것들을입력변수로해서이그림을그리는함수를만들어보자. M(ω) = E ( E + iωη) M M + i ω M M E + E + iωη = ω r M + i ω M M ω r V = Re[ M ( ω)] ρ Re[ M ( ω)] = Q Im[ M ( ω)] ω dv = V V V dω g ω/ω r 복잡한수식을가진함수 - Viscoelasticity 복잡한수식을가진함수 - Viscoelasticity 그림에서보듯이우리의독립변수 (x) 는 ω/ω r 이다. 따라서우리는이값을 omega 라고하자. 그림에서보듯이이주파수는로그스케일로그려져있다. 그래서 linspace 대신에 logspace 를이용해보자. omega = logspace(-,, ); %. 에서 까지 다음은이각각의값에대해서일단 M 을계산해야한다. 입력변수 M, M 를 M, Minf 라고하자. 앞의그룹속도 (Vg) 를정리해보면, V g V = dv V ω d ω 미분값을수치적으로구하기위해 diff() 라는함수를이용하는데, 다음그림과같이유한차분법 (finite-difference) 의일종이다. M = Minf * ( M + i*omega*sqrt(m*minf) )./ ( Minf + i*omega*sqrt(m*minf) ); f f f 3 f 4 f 5 f 6 x x x 3 x 4 x 5 x 6 이렇게하고나면 omega 와그길이가같은 M 벡터가가생겼다. 이벡터는복소수임에주의하라. 우리가구해야하는속도, 감쇄계수는 M 의함수로 V = sqrt( real(m)./ density ); Qinv = imam)./ real(m); 이제탄성파의그룹속도만이남았다. 그런데이수식에는아직우리가배우지않은미분값이들어있다. df dx ( f f) ( x ) at= x x diff(f) = [f-f, f3-f, f4-f3, f5-f4, f6-f3] 주의할점은벡터의길이가하나준다는사실이다. 이는미분이기울기이고기울기는두점이필요하므로당연한결과이다. 복잡한수식을가진함수 - Viscoelasticity 복잡한수식을가진함수 - Viscoelasticity 그래서 diff() 를써준후에마지막값을하나더해야만우리가원하는길이의벡터를구할수있다. 한가지주의할점은여기서 Δx 는 omega 의차이인데, 우리가 logspace 를썼기때문에이것이일정하지않다는사실이다. 따라서 x 축값의차이를각각계산해야한다. dv = diff(v); dv = [dv dv(end)]; % 성분하나를더한다. domega = diff(omega); domega = [domega domega(end)]; dvdw = dv./ domega; Vg = (V.* V)./ ( V omega.* dvdw); 주의할점은연산자들이모두각각의성분값의연산으로주어진다는것을명심해야한다. 정리하여함수를만들어보자. function plotviscoelas(m, Minf, density) omega = logspace(-,, ); M = Minf * ( M + i*omega*sqrt(m*minf) )./ ( Minf + i*omega*sqrt(m*minf) ); V = sqrt( real(m)./ density ); Qinv = imam)./ real(m); dv = diff(v); dv = [dv dv(end)]; domega = diff(omega); domega = [domega domega(end)]; dvdw = dv./ domega; Vg = (V.* V)./ ( V omega.* dvdw); semilogx(omega, V, 'b'), hold on semilogx(omega, Qinv, 'r'), semilogx(omega, Vg, 'k:') legend('v','/q','vg') xlabel('{\omega} / {\omega}_r') ylabel('value')
복잡한수식을가진함수 - Viscoelasticity >> plotviscoelas(,4,) Value.5 V /Q Vg.5.5 - - ω / ω r 구간에따라함수가다른경우 π 다음과같은함수가있다고하자. y = sin( x) : x < π π y = cos( x ) : x < π 5π y = + tan( x π ) : π x < 4 주어진구간 [, 5π/4] 에서이함수를그리는 m-file을작성하라. x =linspace(, 5*pi/4, ); y = zeros( size(x) ); for i=:length(x) if(x(i)<pi/) y(i)=sin( x(i) ); elseif(x(i)<pi) y(i)=-cos( x(i) - pi/ ); else y(i)=+tan(x(i)-pi ); end end plot(x,y) 구간에따라함수가다른경우 m-file 추가예제 하지만각각의성분마다 if를써서계산하는것은매우비효율적이다. x =linspace(, 5*pi/4, ); y = zeros( size(x) ); y(x<pi/) = sin( x(x<pi/) ); y(x<pi & x>=pi/) = -cos( x(x<pi & x>=pi/) - pi/ ); y(x>=pi) = +tan( x(x>=pi) - pi); plot(x,y) 이보다더좋은방법은논리곱을이용하는방법이다. x =linspace(, 5*pi/4, ); y = sin(x).* (x<pi/) + ( -cos(x-pi/) ).*(x>=pi/ & x<pi) + (+tan(x-pi)).*(x>=pi); plot(x,y) 이문제에서는우리가과제에서다루었던문제를이용해서 m-file 을만들어보자. ( 과제 #의 6번문제 ) 5개의컬럼으로구성된 ASCII 데이타가주어지고그파일이름을입력변수로한다. 각각의컬럼은 x좌표, y좌표, 비소, 카드늄, 납의농도이다. 함수가행해야할기능은 () 빠진데이타를찾아서그데이타를주어진데이타들에서모두제거하기, () 각농도의평균및분산값계산하기, (3) 각농도들간의상관관계계산하기, (4) 각농도들을격자화자료로만들기, (5) 각농도들의면구조를 surf를이용해출력하기 (subplot을이용해세가지그래프를한그림에서출력 ), (6) 또다른그림을그리는데카드늄의농도로 surf 그래프를그리되, 그높이를나타내는칼라는납의농도를이용하기. 함수의정의 먼저우리는입력변수및출력변수가무엇인지파악해야한다. 입력변수 : 데이타파일의파일이름 (string type) 출력변수 : 각농도의평균값, 분산, 그리고상관관계 데이타파일읽기와빠진데이타제거 첫번째수행해야할일은데이타를읽는것이다. 데이타는 ascii 파일로포맷이우리가원하는데로주어져있다고가정하자. data = load( filename ); 출력변수를각각의값으로출력하려면, 9 개의변수가필요해너무복잡해지므로평균, 분산, 상관관계를벡터및행렬로하여 3 개로만드는것이좋다. 그러면함수의선언부분이정해진다. 편의상함수의이름은 hwprob6 으로하자. function [ m, sd, ccoef ] = hwprob6( filename ) 그다음은이데이타들중에빠진데이타를제거해보자. 빠진데이타들은 3 이라는값으로항상주어진다고가정하자. 과제에주어진데이타는비소의농도에만빠진데이타가있었지만일반적으로모든데이타에서일어날수있으므로모든데이타를다비교해야한다. ind = find( data(:,3) < ^ & data(:,4) < ^ & data(:,5) < ^ ); data = data( ind, : ); 3
평균, 분산, 그리고상관계수구하기 빠진데이타를제거한이후에평균, 분산, 그리고상관계수를구해보자. 출력변수들을 m, sd, ccoef 로이미정해놓았으므로이변수들을써야한다. m = mean( data( :, 3:5) ); sd = std( data( :, 3:5) ); ccoef = corrcoef( data( :, 3:5) ); 물론위와같이하지않고각각의농도값들을다른변수로지정해서각각평균, 표준편차, 상관계수들을구해도된다. 하지만이렇게구하는것이함수가보다간단해진다. 이는프로그래머의선택에달려있다. 평균과분산은그방향을지정하지않은경우컬럼방향으로평균을구한다는사실을이미배웠다. 이는우리가원하는방향이므로따로방향을지정해주지않아도된다. 상관계수의경우도행렬을입력변수로넣는경우각컬럼들의상관계수를구하므로이또한우리가원하는바이다. 격자화자료구하기 이제그래프를그리기위해서격자화자료를만들어보자. 앞에서배웠듯이격자화자료를만들기위해서먼저격자화된좌표 (x, y) 을먼저만들어야한다. 그리고격자화된좌표의갯수가원래좌표들의갯수와비슷한것이좋다. 격자좌표의갯수를먼저구해보자. 원래의좌표들의갯수는 ndata = length( data(:,) ); 그리고 x와 y의좌표들의최대최소값은 minx = min( data(:,) ); maxx = max( data(:,) ); miny = min( data(:,) ); maxy = max( data(:,) ); 그래서데이타가들어가는 box의크기는 lx = (maxx - minx); ly = (maxy - miny); 만약 lx에 nx의데이타가들어가면 ly에는 nx*(ly/lx) 만큼의데이타가들어간다. 따라서총갯수는 nx*nx*ly/lx 가되고이갯수는원래데이타의갯수와비슷해야한다. 격자화자료구하기 격자화자료구하기 앞의결과를종합하면 nx = sqrt( ndata * lx / ly); 가되는데, 이것은갯수이므로반드시양의정수가되어야한다. 따라서 nx = round ( sqrt( ndata * lx / ly) ); ny = round ( nx * ly / lx ); 문제가조금복잡해지긴하지만, 우리는한가지더생각해야할일이있다. 아까구한 box의경계에서는데이타없어서 griddata를이용할경우그값을정하지못하는경우가있다는것을배웠다. 그래서우리는아래그림과같이격자데이타를구하려고한다. y 의경우도마찬가지지만우리는등간격을이용하기로했으므로 (dx=dy) y축의격자데이타는다음과같이주어진다. y = [miny+dx/ : dx : miny+dx*ny + *dx/3]; [x, y] = meshgrid(x, y); 이제격자데이타의 x, y 좌표가주어졌으므로각각의격자자료를구해보자. as = griddata( data(:,), data(:,), data(:,3), x, y); cd = griddata( data(:,), data(:,), data(:,4), x, y); pb = griddata( data(:,), data(:,), data(:,5), x, y); minx nx maxx dx = lx / (nx+); x = [minx+dx/ : dx : minx+dx*nx + *dx/3];?? surf() 를이용한 plot surf() 를이용한 plot 이제모든자료가구해졌고, 이를이용해단지 plot을하면된다. figure; subplot(3,,); surf(x, y, as); xlabel('x'); ylabel('y'); title('arsenic'); shading interp; camlight; lighting phong; axis tight subplot(3,,); surf(x, y, cd); xlabel('x'); ylabel('y'); title('cadnium'); shading interp; camlight; lighting phong; axis tight 마지막은또다른그래프를그려서카드늄의농도로 surface 를그높이 color 는납의농도로표시하는것이다. figure; surf(x, y, cd, pb); xlabel('x'); ylabel('y'); title('cadnium'); shading interp; camlight; lighting phong; axis tight colorbar subplot(3,,3); surf(x, y, pb); xlabel('x'); ylabel('y'); title('lead'); shading interp; camlight; lighting phong; axis tight 4
결과 m 파일 결과 m 파일 - 계속 function [ m, sd, ccoef ] = hwprob6( filename ) data = load( filename ); ind = find( data(:,3) < ^ & data(:,4) < ^ & data(:,5) < ^ ); data = data( ind, : ); m = mean( data( :, 3:5) ); sd = std( data( :, 3:5) ); ccoef = corrcoef( data( :, 3:5) ); ndata = length( data(:,) ); minx = min( data(:,) ); maxx = max( data(:,) ); miny = min( data(:,) ); maxy = max( data(:,) ); lx = (maxx - minx); ly = (maxy - miny); nx = round ( sqrt( ndata * lx / ly) ); ny = round ( nx * ly / lx ); dx = lx / (nx+); x = [minx+dx/ : dx : minx+dx*nx + *dx/3]; y = [miny+dx/ : dx : miny+dx*ny + *dx/3]; [x, y] = meshgrid(x, y); figure; subplot(3,,); surf(x, y, as); xlabel('x'); ylabel('y'); title('arsenic'); shading interp; camlight; lighting phong; axis tight subplot(3,,); surf(x, y, cd); xlabel('x'); ylabel('y'); title('cadnium'); shading interp; camlight; lighting phong; axis tight subplot(3,,3); surf(x, y, pb); xlabel('x'); ylabel('y'); title('lead'); shading interp; camlight; lighting phong; axis tight figure; surf(x, y, cd, pb); xlabel('x'); ylabel('y'); title('cadnium'); shading interp; camlight; lighting phong; axis tight colorbar as = griddata( data(:,), data(:,), data(:,3), x, y); cd = griddata( data(:,), data(:,), data(:,4), x, y); pb = griddata( data(:,), data(:,), data(:,5), x, y); Fourier Transform ( 푸리에변환 ) Fourier Transform ( 푸리에변환 ) 푸리에변환이란무엇인가? 쉽게말하면시계열 (time series) 의주파수영역으로의 mapping이다. 예를들어서시간에변화에따른해수높이의변화를측정했다고하자. 우리가알고있는조석의변화를생각해보면이시계열은아마도 시간정도의주기를가지고있다는사실을유추해볼수있다. Height Height 5 5 5 5 5 3 35 4 Day 5 5 3 4 5 6 7 8 Day 앞의시계열데이타를푸리에변환한결과가오른편에나타나있다. 여기로부터무엇을알수있는가?. 주기가 인데이타 : 백그라운드데이타. 주파수가.7 (/day) 인데이타 : 이는주기가약 365일인데이타 3. 주파수가.9 (/day) 인데이타 : 주기가약.53, 즉 시간 5분인데이타 - 실제조석 Real amplitude Real amplitude Real amplitude x 5 5 5-5 x 4 FT of time series -5-5 - -5 5 5 /day x 5 8 6 4 3 4 5 6 7 /day x -3.88.9.9.94.96.98 /day 기본지식 기본지식 시간 ( 주기 ) 주파수 sec, day hertz, /day 공간 ( 파장 ) 파수 m, micron /m, /micron.5 cos(π 그렇다면시계열함수가두가지의주기를가지는함수로구성된다면어떻게될까? 오른쪽에서보듯이이제시계열함수가주기가 인코사인과주기가 인코사인의합으로주어진다..5.5 cos(π + cos(π 예를들어 cos(π 함수의경우를살펴보자. 이함수는주기가 인함수이므로푸리에변환을하여주파수도메인에서살펴보면그에너지가주파수가 인경우 (/) 에만존재하고나머지는전부 이다. -.5-4 6 8 5 4 3 FT cos(π 이를푸리에변환을하면주파수가 인부분과주파수가.5 (=/) 인부분에서그에너지를가지고나머지에서는 이다. 다시말하면푸리에변환을하면그시계열이가지는에너지의주된주파수성분을알수있다. -.5 - -.5 4 6 8 5 4 3 FT[ cos(π + cos(π] * 주기 : 반복이일어나는시간차이. - -4-3 - - 3 4 cos(π : t가 일때의값이 일때반복된다. cos(π : t가 일때의값이 일때반복된다. - -.5.5.5 5
푸리에급수 복소수와극좌표계 이러한푸리에변환과푸리에급수는어떻게다른가? 푸리에급수 (Fourier series) 는우리가알고있는함수로표현하기어려운함수를우리가잘알고있는함수들의합으로표현하는방법의일종으로 sin() 함수와 cos() 함수를이용한다. 다음의결과가알려져있다. 다시말해서어떠한복잡한함수라고하더라도무한개의사인과코사인함수만있으면그각각의진폭을조정하면정확하게그함수를표현할수있다는것이다. 이것이왜중요한가? 우리가어떠한함수를얻었을때그함수를우리가잘아는사인이나코사인같은함수로근사할수있으면그함수의도함수및근을구하는데매우유용하게쓰일수있다. 푸리에변환은이렇게주어진시계열을푸리에급수로표현했을때그각각의주파수가가지는에너지의크기라고생각할수있다. y( = n= { a cos( n + b cos( n } n n 실제로푸리에변환을알아보기전에복소수와극좌표에대해서되집어보기로한다. 복소수 (Complex number) : z=x+iy x: 실수성분 (real part = Re[z]) y: 허수성분 (imaginary part = Im[z]) 이것을극좌표계 (polar coordinate) 에서표시하면, z = z (cosθ, sinθ ) i z = z e θ z = x + y y θ = phase = tan x Im z θ ( z cosθ, z sinθ ) Re 복소수와극좌표계 푸리에변환의정의 복소켤레 (complex conjugate) iθ * = x iy = z e z = zz* = x + y z 이를이용하면각각의성분은 z + z * z z * x =, y = i Euler formula iθ e = cosθ + i sinθ e iθ = cos( θ ) + i sin( θ ) = cosθ i sinθ iθ e + e cosθ = iθ e e sinθ = iθ iθ *note i e cos πft = positive frequency i + e π ( f ) t π ( f ) t negative frequency Fourier transform of t, f : 실수인독립변수 G(f) : 주파수도메인에서의푸리에변환된함수 ( 복소수 ) : 시간도메인에서의시계열함수 ( 복소수 ) Euler formula 를써서분해하여보면 g 가실수라도 G 는복소수가될수있다. G( f ) = FT{ } = G( f ) = e iπft e iπft = cos(πf i sin(πf 역푸리에변환의정의 푸리에변환의역푸리에변환 Inverse Fourier transform of G(f) t, f : 실수인독립변수 G(f) : 주파수도메인에서의푸리에변환된함수 ( 복소수 ) : 시간도메인에서의시계열함수 ( 복소수 ) Euler formula 를써서분해하여보면 G 가실수라도 g 는복소수가될수있다. = IFT{ G( f )} = f ) = i G( f ) e πft i G( f ) e πft = G( f )cos(πf + i G( f )sin(πf df G(f) 가 의푸리에변환이면 if t is continuous IFT{ G( f )} = IFT{ FT{ } } = + / { t ) + t )} if t is discontinuous 라는결과가다음이충분조건일때성립한다.. 의불연속점이유한개이고. 는유한구간에서유한개의최대, 최소값을가지며, 3. 가절대적으로적분가능해야한다. 참고로위의조건은충분조건이다. 따라서위의조건이충족되면항상만족하고위의조건이충족되지않는경우에위의식이성립하는경우도있다. 6
매트랩에서푸리에변환 이제실제로매트랩에서푸리에변환을시도해보자. 푸리에변환을하는함수는 fft() 와 ifft() 이다. 이함수는discrete Fourier Transform을수행하는함수로주어진것이시계열이던공간함수이건관계없이푸리에변환를한다. 다시말해서주어진시계열이나공간함수의간격 ( 또는 dx) 을 로가정하고계산한다. 따라서주어진시계열을푸리에변환을하여주파수계열을얻었다면그각각의값들의위치를따로계산해주어야한다. 매트랩에서푸리에변환 여기서우리는재미있는사실을발견할수있다. t=[,.,.4,,.] f=[, 5,, 5,,5] 따라서우리가촘촘히데이타를얻으면주파수사이의거리는넓어진다. 그리고 t=[,.,.4,,] f=[,.5,.,.5,,5]. (N-) (N-) t f 결과적으로다음의결론을얻는다.. 주파수사이의간격은시계열의전체길이가길수록짧아진다.. 데이타를얻을수있는주파수의대역은시계열간격이길수록짧아진다. 시간도메인과주파수도메인은서로반대되는성질이있다. 만약 t=[,.,.4,,] 이라면, f=[,.5,,.5,,5] 이된다. 나이퀴스트주파수 (Nyquist f) 다음의사인함수을제대로기술하려면한주기당몇개의데이타가필요할까? 아래그림에서보듯이주기당 개 ( 빨간원 ) 의데이타값만있을경우는제대로함수를표시하기힘들다. 4개 ( 파란 x) 인경우는상대적으로잘표현하고있다. 따라서하나의주기에 개의데이타가있는경우보다더작은경우는그주파수는표현이불가능하다. 이때의주파수를나이퀴스트주파수라한다..5 -.5 sin(x) 나이퀴스트주파수 (Nyquist f) f Nyq = n t = [,,,..., n ] f = [,,,..., n n n 그런데실제주파수의영역은나이퀴스트주파수의두배이다. 하지만우리는나이퀴스트주파수다음에는실제데이타가없다는사실을알고있다. 그래서나이퀴스트와그이후의값들은실제로음의주파수가된다. 음의주파수가왜존재하는가? 이는푸리에변환의수학적정의때문에도출되는값으로물리적으로의미는없다. n n n + t = [,,,..., n ] f = [,,...,,,,,..., ] n n n n n n ] - 이러한작업을해주는매트랩함수가있는데이것이 fftshift() 이다. -6-4 - 4 6 x 푸리에변환매트랩예제 짝함수 (Even) 와홀함수 (Odd) 코사인함수를푸리에변환을해보자. t=linspace(, *pi, ); g=cos(; =t()-t(); df=//(-); f=[-5:499]*df; G=fft(g); G=fftshift(G); subplot(,,) plot(t, g); subplot(,,) plot(f, G);.5 -.5 - - 5 5 5 3 35 4 5 4 3 -.6 -.4 -...4.6.8 어떤함수 가 y축을중심으로대칭이면그함수는짝함수 (even function) 이다. 다시말하면 -= 이다. 어떤함수 가원점을중심으로대칭이면그함수는홀함수 (odd function) 이다. 다시말하면 -= - 이다. 홀함수는차수가홀수 (, 3, 5) 인경우라서홀함수라고하고짝함수는차수가짝수 (, 4, 6) 인함수이라서짝함수라고부른다..8.6.4. -. -.4 -.6 -.8 - - -.5.5.9.8.7.6.5.4.3.. - -.5.5 7
짝함수와홀함수의성질 a O( df = ; O( df = ; E ( E ( = even O ( O ( = even a E( df = E( df E( df = E( df a a a E( O( = odd Theorem 어떠한함수이든그함수는짝함수와홀함수의합으로표현될수있다. E( = / { + } E( = E( O( = / { } O( = O( = E( + O( 짝함수와홀함수의이용 다시푸리에변환의정의로돌아가보자. G( f ) = FT{ } = = E( G( f ) = = O( G( f ) = G( f ) = E( cos( πf = G( f ) : 실수짝함수 G( f ) = i O( sin( πf = G( f ) : 허수홀함수 = ie( G( f ) = i E( cos(πf + G( f ) = i E( cos( πf = G( f ) : 허수짝함수 = io( G( f ) = i O( cos(πf + G( f ) = E( cos(πf i E( sin(πf = O( cos(πf i O( sin(πf = i O( sin(πf E( sin(πf = i E( cos(πf O( sin(πf = O( sin( πf = G( f ) : 실수홀함수 E( cos(πf O( sin(πf e iπft Matlab 을이용한선형대수계산 앞에서우리는매트랩이행렬연산을중심으로그개발이시작되었다는것을배웠다. 그러면여기서일반적인선형행렬연산을매트랩을이용해서어떻게수행하는지몇가지의예를살펴보고자한다. 벡터와행렬의 Norm The p-norm of a vector x x = p x p i i / p is computed by norm(x,p). This is defined by any value of p >, but the most common values of p are,, and. The default value is p =, which corresponds to Euclidean length. The p-norm of a matrix A, Ax p A = max p x x p / p can be computed for p =,, and by norm(a,p). Again, the default value is p =. Ax=b 의풀이 가장일반적인선형방정식 Ax=b 의풀이를생각해보자. 만약 A 가 (m,n) 의크기를가지는행렬인경우, () m > n : 주어진식이해보다많다 (overdetermined). 이경우에모든식을만족시키는해가없는경우가이므로, 실제구해진해는차이값 ( 에러 ) 들의합을최소화하는해를구한다. () m = m : 식과해의갯수가일치하는경우로이경우에는대부분답을구할수있다. 다만주어진식들중에서같은것이있는경우답이정할수없다. (3) m < n : 주어진식이해보다적을경우로 (underdetermined), 하나의해를구할수없고관계식이주어진다. 매트랩에서는 / 이나 \ 연산자를이용해서구한다. Ax=b 의풀이 만약행렬 A 가정사각행렬 (square matrix) 인경우에는정확한 Ax=b 의풀이는 A 의역행렬을구하여, 양변에곱하면된다. x = inv(a) * b; 하지만행렬 A 가역행렬을가지지않는경우, 다시말해서 A 의행렬식이 인경우, 다시말해서 A 의행렬의열들이모두서로독립이아닌경우는역행렬을가지지않고따라서해는존재하지않는다. 그런데행렬식이 은아니지만아주 에가까운경우를생각해보자. A =, 4 B = 4. 위의 B 의경우는수학적으로역행렬을가지지만, 실제컴퓨터에서는역행렬을계산할수도있고, 못할수도있다. 이렇게행렬이얼마나 singluar 에가까운가를알아보는데는 condition number 를쓴다. >> cond(b) 8
Ax=b 의풀이 만약행렬 A 에어떤행렬 B 를곱했더니그결과가단위행렬이되었다고하자. Ax = b, BA = x 3 = x, 3 Bb = 당연히행렬 B 는 A 의역행렬이다. 그런데이렇게역행렬을곱하고나면미지수 x 앞에있는행렬이단위행렬이되고실제그해들은식의오른쪽에있는 Bb 가바로답이된다. 따라서, 선형방정식의여러가지풀이방법은행렬 A 에어떤값을곱해서그결과가근사적으로단위행렬이되게하는것이그기본목적이다. Eigenvalue and Eigenvector 사각행렬을다룰때, 우리는흔히고유값과고유벡터란용어를흔히쓴다. 이것이무엇인지알아보자. 고유값은어떤행렬 A가 Ax = λx 를만족시킬때그때의스칼라값 λ를고유값그리고벡터 x를고유벡터라고한다. 이렇게하기위해서우리는 A 를고유값변환이라는것을이용해행렬을분해하게된다. 그결과만제시하면, AV = VΛ, A = VΛV 여기서중간에주어진행렬은대각행렬이된다. 그리고여기서우리는한가지의사실을알아야하는데, 벡터의좌표변환을기억해보면, 다음과같이된다. x = Tx 비슷한방법으로행렬의좌표변환을알아보면다음과같이된다. A = TAT Eigenvalue and Eigenvector Eigenvalue and Eigenvector 따라서앞의결과로부터, 우리는 A 라는행렬은대각행렬 Λ 을 3 차원좌표변환을한값이라는사실을알수가있다. 다음예제를살펴보자. 여기서우리는 (σ xx =4 σ yy =8 σ yx =3.46) 이라는것에서시작해서주응력을구해보자. 그러면이대각행렬의성분들은무엇이되는가? 그것은바로아까말했던고유값들이다. 그리고좌표변환에이용된행렬 V 는그행렬의고유벡터들의모임이다. 다시이를정리하면무엇이되는가? 어떠한행렬을고유분해를하면 ( 회전변환 ) 우리는대각행렬을얻게된다. 이내용은우리에게매우익숙하다. 응력의경우를예를들어보자. 우리는응력이텐서, 다시말해서행렬이라는사실을알고있다. 이것을회전변환을해서대각성분만남는다는말은무엇인가. 응력장의좌표를회전시켜서주응력 ( 대각성분 ) 만남긴다는말이다. 다시말해서우리가원하는 eigenvalue 는주응력성분이되고, eigenvector 는그주응력방향이된다. 3 MPa MPa MPa y x -6 τ (σ, τ) 6 σ 매트랩에서 eig 라는명령어가행렬의고유치와고유벡터를구하는함수이다. MPa σ = 6-4 cos(6 ) = 4 MPa τ = 4 sin(6 ) = 3.46 MPa Eigenvalue and Eigenvector Symbolic Math 이번에는고유치와고유벡터를구해서, (σ xx =4 σ yy =8 σ yx =3.46) 의결과로부터주응력을계산해보자. t 8 3.46 σ = 3.46 4 t 8 λ 3.46 σx λx = x = 3.46 4 λ (8 λ)(4 λ) 3.46 = λ = or 이계산은매트랩에서 ei) 란함수로수행할수있다. >> [v d]=ei[8 3.46; 3.46 4]) v =.5 -.87 -.87 -.5 d =.. 또다른매트랩에서유용한 tool중의하나는 symbolic math의기능이다. 이는실제숫자를이용하지않고문자로계산하는방법이다. 예를들어 >> solve( 'a*x^ + b*x + c = ', 'x') //a*(-b+(b^-4*a*c)^(/)) //a*(-b-(b^-4*a*c)^(/)) 이는여러분이잘아는이차방정식의근의공식이다. 그러면 3차방정식의근은? >> solve( 'a*x^3 + b*x^ + c*x + d = ', 'x') /6/a*(36*b*c*a-8*d*a^-8*b^3+*3^(/)*(4*a*c^3-c^*b^- 8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)-/3*(3*a*c-b^)/a/(36*b*c*a-8*d*a^- 8*b^3+*3^(/)*(4*a*c^3-c^*b^-8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)-/3*b/a -//a*(36*b*c*a-8*d*a^-8*b^3+*3^(/)*(4*a*c^3-c^*b^- 8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)+/3*(3*a*c-b^)/a/(36*b*c*a-8*d*a^- 8*b^3+*3^(/)*(4*a*c^3-c^*b^-8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)- /3*b/a+/*i*3^(/)*(/6/a*(36*b*c*a-8*d*a^-8*b^3+*3^(/)*(4*a*c^3-c^*b^- 8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)+/3*(3*a*c-b^)/a/(36*b*c*a-8*d*a^- 8*b^3+*3^(/)*(4*a*c^3-c^*b^-8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)) -//a*(36*b*c*a-8*d*a^-8*b^3+*3^(/)*(4*a*c^3-c^*b^- 8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)+/3*(3*a*c-b^)/a/(36*b*c*a-8*d*a^- 8*b^3+*3^(/)*(4*a*c^3-c^*b^-8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)-/3*b/a- /*i*3^(/)*(/6/a*(36*b*c*a-8*d*a^-8*b^3+*3^(/)*(4*a*c^3-c^*b^- 8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)+/3*(3*a*c-b^)/a/(36*b*c*a-8*d*a^- 8*b^3+*3^(/)*(4*a*c^3-c^*b^-8*b*c*a*d+7*d^*a^+4*d*b^3)^(/)*a)^(/3)) 9
Symbolic Math 앞에서보듯이그답이무지하게복잡하다. 이러한경우에는그식을다른변수로치환할수있다. >> sol = solve( 'a*x^3 + b*x^ + c*x + d = ', 'x'); 그런후에 a,b,c,d 의값을정하고답을구하고싶으면, >> a=; b=; c=; d=-; >> eval(sol). -.5 +.866i -.5 -.866i 이경우 sol 은식을가진 symbolic object 가되고, 여기에값을대입해해를구하고싶을땐, 각각의변수에값을대입하고 eval() 을이용해서그함수를구하면된다. Symbolic Math 단순히식을푸는것 (solve) 뿐만이아니라다른수식계산과매트랩에서가능한데그중몇가지만살펴보면, 함수의미분과적분 >> diff( sin(x), x ) cos(x) >> int('sin(a*x)','x') -/a*cos(a*x) >> int('sin(x)','x',, pi) >> int('sin(a*x)','x',, pi) -(cos(pi*a)-)/a