수치해석 161009 Ch3. Programming with MATLAB
3.1 M- 파일 (1/5) 스크립트파일 - 일련의 MATLAB 명령어를구성되어저장된 M- 파일이다. %scriptdemo.m g=9.81; m=68.1; cd=0.25; t=12; v=sqrt(g*m/cd)*tanh(sqrt(g*cd/m)*t) >> scriptdemo v = 50.6175
3.1 M- 파일 (2/5) 함수파일 - function 이라는단어로시작하는 M- 파일이다. function outvar = funcname(arglist) %helpcomments statements outvar = value; outvar = 출력변수의이름 funcname = 함수의이름 arglist = 함수의인수목록 helpcomments = 사용자가제공하는함수에관한정보 statements = value를계산하여그값을 outvar에배정하는문장
예제 3.2 (1/3) Q. 번지점프하는사람의자유낙하속도를함수파일을사용하여구하라. dv cd 2 gm = g v gcd v( t) = tanh t dt m cd m function v = freefallvel(t, m, cd) % freefallvel: bungee velocity with second-order drag % v=freefallvel(t,m,cd) computes the free-fall velocity % of an object with second-order drag %input: % t=time(s) % m=mass(kg) % cd=second order drag coefficient (kg/m) %output: % v=downward velocity (m/s) g=9.81; %accelearation of gravity v=sqrt(g*m/cd)*tanh(sqrt(g*cd/m)*t);
예제 3.2 (2/3) Q1. 68.1 kg 인사람의 12 초후의속도를구하려면 >> freefallvel(12, 68.1, 0.25) ans = 50.6175 Q2. 100 kg 인사람의 8 초후의속도를구하려면 >> freefallvel(8, 100, 0.25) ans = 53.1878
예제 3.2 (3/3) 도움설명을불러내려면다음과같이입력한다. >> help freefallvel freefallvel: bungee velocity with second-order drag v=freefallvel(t, m, cd) computes the free-fall velocity of an object with second-order drag. input: v = downward velocity (m/s)
3.1 M- 파일 (3/5) - 함수 M- 파일은 2개이상의결과를반환할수있다. 예 ) 벡터의평균과표준편차의계산 function [mean, stdev] = stats(x) n=length(x); mean=sum(x)/n; stdev=sqrt(sum((x-mean).^2/(n-1))); >> y=[8 5 10 12 6 7.5 4]; >> [m,s] =stats (y) m = 7.5000 S = 2.8137
3.1 M- 파일 (4/5) 부함수 (subfunctions) - 함수가다른함수를부를수있다. 이러한함수는 M- 파일을구분하여작성할수도있고, 한개의 M 파일에포함시킬수도있다. function v= freefallsubfunc(t, m, cd) v=vel(t, m, cd); 주함수 function v=vel(t, m, cd) g=9.81; v=sqrt(g*m/cd)*tanh(sqrt(g*cd/m)*t); 부함수
3.1 M- 파일 (5/5) >> freefallsubfunc (12, 68.1, 0.25) ans = 50.6175 >> vel (12, 68.1, 0.25)??? Undefined command/function vel.
3.2 입력 - 출력 (1/4) input 함수 - 사용자로하여금명령창에서직접입력하도록한다. m = input ('Mass (kg): ') name = input ('Enter your name: ', 's') % 문자열을받는경우
3.2 입력 - 출력 (2/4) disp 함수 - 어떤값을손쉽게나타낸다. disp(' ') disp('velocity (m/s): ') % 문자열을나타내는경우
3.2 입력 - 출력 (3/4) fprintf 함수 - 정보를표현할때추가적인제어를제공한다. fprintf ('format', x, ) % 포맷코드와제어코드를넣어서나타내는경우 >> fprintf( The velocity is %8.4f m/s\n, velocity) The velocity is 50.6175 m/s
3.2 입력 - 출력 (4/4) 표 3.1 fprintf 함수에서사용하는포맷코드와제어코드 포맷코드설명 %d %e %E %f %g 정수포맷 e를사용하는과학포맷 E를사용하는과학포맷소수포맷 %e 나 %f 중간단한포맷 제어코드설명 \n \t 새로운줄로시작탭 3장 MATLAB 프로그래밍
예제 3.3 (1/2) Q. 예제 3.2 에서와같이번지점프하는사람의자유낙하속도를계산하라. 입출력으로 input 과 disp 함수를사용하라. function freefalli % freefalli: interactive bunge velocity % freefalli interactive computation of the free-fall velocity of an object % with second-order drag. g=9.81; % acceleration of gravity m=input('mass(kg):'); cd=input('drag Coefficient(kg/m):'); t=input('time(s):'); disp(' ') disp('velocity (m/s):') disp(sqrt(g*m/cd)*tanh(sqrt(g*cd/m)*t)) minusvelocity= -sqrt(g*m/cd)* tanh(sqrt(g*cd/m)*t); fprintf('the velocity is %8.4f m/s\n', minusvelocity) name = input ('Enter your name: ', 's'); disp('your name is'); disp(name)
예제 3.3 (2/2) 명령창에서다음과같이입력한다. >> freefalli Mass (kg): 68.1 Drag coefficient (kg/m): 0.25 Time (s): 12 Velocity (m/s): 50.6175 The velocity is -50.6175 m/s Enter your name: Kim Your name is Kim
3.3 구조프로그래밍 (1/11) 명령을연속적으로수행하지않는것을허용하는구문 - 판정 ( 또는선택 ): - 루프 ( 또는반복 ): 판정 [if 구조 ] ): 판정에기초를둔흐름의분기점이다. ): 반복을허용하는흐름의루프이다. if condition statements
3.3 구조프로그래밍 (2/11) function grader(grade) if grade >= 60 disp('passing grade:') disp(grade) >> grader(95.6) passing grade: 95.6000 < 파일편집기 > < 명령창 >
3.3 구조프로그래밍 (3/11) Cf) 에러함수 error (msg) function f = errortest(x) if x==0, error(' zero value encountered'), f=1/x; >> errortest(10) ans = 0.1000 >> errortest(0)??? Error using ==> errortest zero value encountered
3.3 구조프로그래밍 (4/11) [ 논리조건 ] value1 relation value2 ~ (Not) - 논리적부정을나타낼때사용한다. & (And) - 두식에서논리적곱을나타낼때사용한다. (Or) - 두식에서논리적합을나타낼때사용한다.
3.3 구조프로그래밍 (5/11) < 복잡한논리식의단계별계산 >
3.3 구조프로그래밍 (6/11) [if else 구조 ] [if elseif 구조 ] if condition statements1 else statements2 if condition1 statements1 elseif condition2 statements2 else statementselse
예제 3.4 (1/3) Q. 내장함수인 sign 함수와같은기능을갖도록 if else 구조를사용하여 mysign 함수를작성하라. 풀이 ) 내장함수인 sign 함수의기능을알아보자. >> sign(25.6) ans = 1 >> sign(-0.776) ans = -1 >> sign(0) ans = 0
예제 3.4 (2/3) function sgn = mysign(x) % mysign(x) returns 1, -1, and 0 for positive, negative, and zero values, respectively. if x > 0 sgn = 1; elseif x < 0 sgn = -1; else sgn = 0;
예제 3.4 (3/3) 명령창에서다음과같이확인할수있다. >> mysign(25.6) ans = 1 >> mysign(-0.776) ans = -1 >> mysign(0) ans = 0
3.3 구조프로그래밍 (7/11) 루프 [for 구조 ] for index = start:step:finish statements
예제 3.5 (1/2) Q. 순차곱셈을위한 for 루프 풀이 ) factorial 함수와같은기능을갖도록 for 루프를사용하여프로그램을작성한다. >> factorial(5) ans = 120
예제 3.5 (2/2) function fout = factor(n) % computes the product of all integers from 1 to n. % x = 1; for i=1:n x = x * i; fout = x; 명령창에서다음과같이확인할수있다. >> factor(5) ans = 120
3.3 구조프로그래밍 (8/12) [ 벡터화 ] i = 0; for t = 0:0.02:5 i = i + 1; y(i) = cos(t); 위와같은 for 루프를다음과같이벡터화할수있다. t = 0:0.02:5; y = cos(t);
3.3 구조프로그래밍 (9/11) [ 메모리사전할당 ] t = 0:.01:5; for i = 1:length(t) if t(0)>1 y(i)=1/t(i); else y(i) =1; 벡터화를이용하여메모리영역을미리할당 (ones 나 zeros 주로이용 ) t = 0:.01:5; y = ones(size(t)); for i = 1:length(t) if t(0)>1 y(i)=1/t(i);
3.3 구조프로그래밍 (10/11) [while 구조 ] while condition statements 5 >> x = 8; >> while x > 0 x = x - 3; disp(x) 2-1 while 과 사이의조건이참인동안에만반복한다. (while 루프의종료는시작위치에서조건이맞지않을경우발생한다.)
3.3 구조프로그래밍 (11/11) [while break 구조 ] while (1) statements if condition, break, statements >> x = 8; >> while(1) x = x - 3; disp(x); if x<0, break, x = 5 x = 2 x = -1 while 과 사이의조건이참인동안반복하나, while 루프내의어떤위치에서도종료할수있다.
3.4 내포화와들여쓰기 내포화 (nesting) - 다른구조안에구조를배치하는것이다.
예제 3.6 (1/3) Q. 2 차방정식을풀기위해내포화를이용함 f ( x) = ax + bx+ 2 c x = b ± b 2 4ac 2a
예제 3.6 (2/3) function quadroots(a, b, c) % quadroots : roots of a quadratic equation % quadroots(a, b, c) : real and complex roots % of quadratic equation % input: % a = second-order coefficient % b = first-order coefficient % c = zero-order coefficient % output: % r1: real part of first root % i1: imaginary part of first root % r2: real part of second root % i2: imaginary part of second root
예제 3.6 (2/3) if a == 0 %special cases if b~= 0 %single root r1=-c/b else %trivial root error(' Trivial solution. Try again! ') else d=b^2 4*a*c; %determinant if d>=0 % real roots r1 = (-b+sqrt(d))/(2*a) r2 = (-b-sqrt(d))/(2*a) else %complex roots r1 = -b/(2*a) i1 = sqrt(abs(d))/(2*a) r2=r1 i2=-i1
예제 3.6 (3/3) 명령창에서다음과같이확인할수있다. >> quadroots(1,1,1) r1 = -0.5000 i1 = 0.8660 r2 = -0.5000 i2 = -0.8660
예제 3.6 (3/3) 명령창에서다음과같이확인할수있다. >> quadroots(1,5,1) r1 = -0.2087 r2 = -4.7913 >> quadroots(0,0,0)??? Error using ==> quadroots Trivial solution. Try again!
3.5 M- 파일로의함수전달 (1/8) 무명함수 - M- 파일을만들지않고간단한함수를생성할수있게한다. 명령창에서다음과같은구문을사용한다. fhandle = @(arglist) expression >> fl = @(x,y) x^2 + y^2; >> fl (3,4) ans = 25
3.5 M- 파일로의함수전달 (2/8) inline 함수 - Matlab 7 이전에서무명함수와같은역할수행. funcname = inline('expression',`'var1', 'var2', ) >> f1 =inline('x^2 + y^2', 'x', 'y') f1 = inline function: f1(x,y) = x^2 + y^2 >> f1(3,4) ans = 25
3.5 M- 파일로의함수전달 (3/8) function 함수 - 다른함수에작동하는함수. 예 ) fplot (fun, lims) >> vel= @(t) sqrt(9.81*68.1/0.25)* tanh(sqrt(9.81*0.25/68.1)*t); >> fplot(vel, [0 12])
예제 3.7 Q. 어떤범위에서함수의평균값을구하기위한 M- 파일만들기 gm gm vt () = tanh t C d C d 풀이 ) t=0 에서 t=12 까지의범위에서함수값을그래프로그릴수있다. >> t=linspace(0, 12); >> v= sqrt(9.81*68.1/0.25)*tanh(sqrt(9.81*0.25/68.1)*t); >> mean(v) ans = 36.0870
3.5 M- 파일로의함수전달 (4/8) function favg = funcavg(a, b, n) % input: % a= lower bound of range % b= upper bound of range % n= number of intervals % output: % favg = average value of function x = linspace(a,b,n); y=func(x); favg=mean(y); function f=func(t) f=sqrt(9.81*68.1/0.25)*tanh(sqrt(9.81*0.25/68.1)*t);
3.5 M- 파일로의함수전달 (5/8) function favg = funcavg(f, a, b, n) % input: % a= lower bound of range % b= upper bound of range % n= number of intervals % output: % favg = average value of function x = linspace(a,b,n); y=f(x); favg=mean(y); function f=func(t) f=sqrt(9.81*68.1/0.25)*tanh(sqrt(9.81*0.25/68.1)*t);
3.5 M- 파일로의함수전달 (6/8) 명령창에서다음과같이확인할수있다. >> vel= @(t) sqrt(9.81*68.1/0.25)*tanh(sqrt(9.81*0.25/68.1)*t); >> funcavg(vel, 0, 12, 60) ans = 36.0127 >> funcavg(@sin, 0, 2*pi, 180) ans = -6.3001e-017
3.5 M- 파일로의함수전달 (7/8) 매개변수의전달 - 매개변수에새로운값을취할때편리함. - function 함수의마지막입력인수에 varargin 추가함. [ funcavg 의수정 ] function favg = funcavg (f, a, b, n, varargin) x = linspace(a,b,n); y = f(x, varargin{:}); favg = mean(y)
3.5 M- 파일로의함수전달 (8/8) 명령창에서다음과같이확인할수있다. >> vel= @(t, m, cd) sqrt(9.81*m/cd)*tanh(sqrt(9.81*cd/m)*t); >> funcavg(vel, 0, 12, 60, 68.1, 0.25) %m=68.1, cd=0.25 ans = 36.0127 >> funcavg(vel, 0, 12, 60, 100, 0.28) %m=100, cd=0.28 ans = 38.9345
3.6 사례연구 : 번지점프하는 번지점프하는사람의속도 속도 (1/3) 수학적모델 : dv cd 2 = g v dt m Euler 법 : v i+1 = v i + dv dt i t
3.6 사례연구 : 번지점프하는 번지점프하는사람의속도 속도 (2/3) function v = velocity1(dt, ti, tf, vi) % velocity1: Euler solution for bungee velocity % v = velocity1(dt, ti, tf, vi) % Euler method solution of bunge jumper velocity % inputs: % dt = time step (s), ti = initial time (s), tf = final time (s) % vi = initial value of depent variable (m/s) % output: % v = velocity at tf (m/s) t = ti; v = vi; n = (tf - ti)/dt; for i=1:n dvdt= deriv(v); v = v + dvdt *dt; t = t + dt; v = v; function dv = deriv( v) dv = 9.81 (0.25 / 68.1)*v^2;
3.6 사례연구 : 번지점프하는 번지점프하는사람의속도 속도 (3/3) 명령창에서다음과같이확인할수있다. >> velocity1(0.5, 0, 12, 0) %0.5초간격, 24번을계산 ans = 50.9259 >> velocity1(0.01, 0, 12, 0) %0.01초간격, 1,200번을계산 ans = 50.6239 >> velocity1(0.001, 0, 12, 0) %0.001초간격, 12,000번을계산 ans = 50.6181