Chap. 상미분방정식의해법 CAE 기본개념소개 Euler법 Heun 법 중점법 Runge-Kutta법 1 Chap. 미분방정식 상미분방정식 상미분방정식 (Ordnar Dfferental Equaton; ODE) One-step method Euler 법 (Euler s method) Heun 법 (Heun s method) 중점법 (Mdpont method) Runge-Kutta 법 (RK method) Multstep method Adaptve Runge-Kutta method Chap. 미분방정식 1
Euler s Method 3 Chap. 미분방정식 Euler 법에서의오차해석 Talor 급수전개 절단오차 Ex) P. 59 예제.1 ' 4e.8t.5 = at t=, step sze =1 ( t) 4 1.3 ( e.8t e.5t ) e.5t 4 Chap. 미분방정식
Euler 법 M Fle functon [t,] = Eulode(ddt,tspan,,h) % [t,] = Eulode(ddt,tspan,,h): % uses Euler's method to ntegrate an ODE % nput: % ddt = name of the M-fle that evaluates ODE % tspan = [t, tf] where t and tf = ntal and % fnal values of ndependent varable % = ntal value of dependent varable % h = step sze % output: % t = vector of ndependent varable % = vector of soluton for dependent varable t = tspan(1); tf = tspan(); t = (t:h:tf)'; n = length(t); t = tspan(1); tf = tspan(); t = (t:h:tf)'; n = length(t); % f necessar, add an addtonal value of t % so that range goes from t = t to tf f t(n) < tf t(n+1) = tf; n = n+1; end = *ones(n,1); %preallocate to mprove effcenc for = 1:n-1 %mplement Euler's method (+1) = () + ddt(t(),())*(t(+1)-t()); end 5 Chap. 미분방정식 Heun s method Euler 법 : 구간의시작점에서의도함수를구간전체에적용함으로인해오차발생 시작점에서의도함수및끝점에서의도함수이용 예측자 (predctor) 와수정자 (corrector) 사용 predctor corrector 6 Chap. 미분방정식 3
Heun s method predctor : corrector : 개선된결과를얻기위한 Heun 법수정자의그래픽표현 수정자의수렴에대한종료판정 : Ex) P. 6 예제. 7 Chap. 미분방정식 중점법 (Mdpont method) predctor : corrector : 8 Chap. 미분방정식 4
5 9 Chap. 미분방정식 Runge-Kutta 법 n = 1, a 1 = 1: Euler 법 : n 차 Runge-Kutta 법 ( 보통 차 /4 차많이사용 ) 고차도함수를구하지않고도 Talor 급수방법이가지는정확도를가짐 h x f h ), ( 1 1 Chap. 미분방정식 Runge-Kutta 법 ( 차 ) 차 Runge-Kutta 표현상수 a 1, a, p 1 그리고 q 11 를결정하기위하여 차 Talor 급수와같다고가정 Where dt d t f t t f t f ), ( ), ( ), ( 1! ), ( ), ( h t f h t f 식 3 개, 미지수 4 개 불능 a 의배정 ( 가정 ) 후나머지계수결정 차 RK 법은 a 의배정에따라무한히많은종류가있다
Runge-Kutta 법 ( 차 ) 반복이없는 Heun 법 a = 1/ a 1 = ½, p 1 = q 11 = 1 Where Heun 법에서수정자의반복이없는경우에해당 중점법 a = 1 a 1 =, p 1 = q 11 = 1/ Where 1 ( a1k1 ak ) h (k 1, k : 구간의시작 / 끝에서의기울기 ) 11 Chap. 미분방정식 Runge-Kutta 법 (4 차 ) 전형적인 4 차 Runge-Kutta 법 가장보편적으로사용되는 4 차 RK 법 Ex) P. 68 예제.3 Smpson 1/3 공식과유사함 간격에대한평균기울기를개선하기위해여러기울기값을추정 (Heun 법과유사 ) Where 1 Chap. 미분방정식 6
연립미분방정식의필요성 고차미분방정식의풀이 예 ) 1 자유도진동방정식 : 차미분방정식 d x dx m c kx dt dt (1) 종속변수의 1 차도함수를새로운변수로정의 () 식 () 를식 (1) 에대입하면 dv m cv kx dt : 연립미분방정식 13 Chap. 미분방정식 연립미분방정식의풀이 n차미분방정식 n개의연립미분방정식으로변환하여풀이시작값 t에서 n개의초기조건필요예 ) 번지점프를하는사람의속도와위치를모두결정해야하는경우 dx dv c d v g v dt dt m 초기조건 : x() = v() = 연립미분방정식의수치해법연립 Euler법연립 Runge-Kutta법 14 Chap. 미분방정식 7
Euler 법을사용한연립상미분방정식풀이 Ex) P. 611 예제.4 ( 번지점프문제 ) Analtc soluton: v( t) dx dv c d v g v x() = v() =, t = 1s, h = s dt dt m gm gc m d tanh t gcd x( t) ln cosh t c d m c d m x v 1 1 x h v h 15 Chap. 미분방정식 Euler 법을사용한연립상미분방정식풀이 Ex) P. 611 예제.4 ( 번지점프문제 ) t x true v true x Euler v Euler t (x) t (v) 4 6 8 1 19.1663 71.934 147.946 37.514 334.178 18.79 33.1118 4.76 46.9575 49.414 39.4 11.674 4.664 35.44 19.6 36.4137 46.983 5.18 51.313 1. 45.45 4.5 13.83 8.7 4.76 9.97 1.3 6.86 3.83 간격크기가커서결과가정확하지않음두번째반복을수행하기까지 x Euler 는 으로계산됨. 간격크기를줄이면결과를개선시킬수있음. 고차방법을사용하면상대적으로큰간격에대해서도좋은결과를얻을수있음 16 Chap. 미분방정식 8
연립 Runge-Kutta 법 R-K 법을사용한연립상미분방정식풀이 고차 RK 법은연립방정식의해를구하는데적용가능. 기울기를구하는데주의해야함. 간격의시작점에서모든변수에대해기울기 ( k 1 's ) 결정 k 1 's 를이용하여간격의중점에서의기울기 ( k 's ) 예측. 중점에서의새로운기울기 (k 3 's) 를예측. 간격의끝점에서의기울기 (k 4 's) 를예측. 모든 k 가증분함수에서합성되어간격끝에서의함수값이결정. Where k1 f ( t, ) 1 1 k f t h, k1h 1 1 k3 f t h, k h k 4 f t h, k 3h 17 Chap. 미분방정식 연립 Runge-Kutta 법 Ex) P. 61 예제.5 ( 번지점프문제 ) Analtc soluton: v( t) dx dv c d v g v x() = v() =, t = 1s, h = s dt dt m gm gc m d tanh t gcd x( t) ln cosh t c d m c d m 18 Chap. 미분방정식 9
연립 Runge-Kutta 법 Ex) P. 61 예제.5 ( 번지점프문제 ) k f ( t, ) 1 k 1 f t h, 1 k h 1 k 3 1 f t h, 1 k h t k 4 f h, k 3h 19 Chap. 미분방정식 연립 Runge-Kutta 법 Ex) P. 61 예제.5 ( 번지점프문제 ) t x true v true x RK4 v RK4 t (x) t (v) 4 6 8 1 19.1663 71.934 147.964 37.514 334.178 18.79 33.1118 4.76 46.9575 49.414 19.1656 71.9311 147.951 37.514 334.166 18.756 33.995 4.547 46.9345 49.47.4.1.4..5.19.37.51.49.38 Chap. 미분방정식 1
연립 Runge-Kutta 법 M Fle functon [tp,p] = rk4ss(ddt,tspan,,h,varargn) % rk4ss: 4th-order R-K for a sstem of ODEs f nargn<4,error('at least 4 nput arguments requred'),end f an(dff(tspan)<=),error( check tspan order'), end n = length(tspan); t = tspan(1);tf = tspan(n); f n == t = (t:h:tf)'; n = length(t); f t(n)<tf t(n+1) = tf; n = n+1; end else t = tspan; end tt = t; (1,:) = ; np = 1; tp(np) = tt; p(np,:) = (1,:); =1; whle(1) tend = t(np+1); hh = t(np+1) - t(np); f hh>h,hh = h;end whle(1) f tt+hh>tend,hh = tend-tt;end k1 = ddt(tt,(,:),varargn{:})'; md = (,:) + k1.*hh./; k = ddt(tt+hh/,md,varargn{:})'; md = (,:) + k*hh/; k3 = ddt(tt+hh/,md,varargn{:})'; end = (,:) + k3*hh; k4 = ddt(tt+hh,end,varargn{:})'; ph = (k1+*(k+k3)+k4)/6; (+1,:) = (,:) + ph*hh; tt = tt+hh; =+1; f tt>=tend,break,end end np = np+1; tp(np) = tt; p(np,:) = (,:); f tt>=tf,break,end end 1 Chap. 미분방정식 연립 Runge-Kutta 법 M Fle 실행 Step 1. 상미분방정식정의를위한 m-fle 작성 functon d = ddtss(t,) d = [(); 9.81-.5/68.1*()^]; Step. rk4ss.m-fle 실행 ( 구간 ~1, 간격 ) Step 3. 그래프를통한결과확인 ( 간격세분화 ) >> [t ] = rk4ss(@ddtss, [ 1], [ ],.1); >> plot(t,(:,1),t,(:,), ) >>legend( x(t), v(t) ) >> [t ] = rk4ss(@ddtss, [ 1], [ ], ); >> dsp ([t (:, 1) (:, )]). 19.1656 18.756 4. 71.9311 33.995 6. 147.951 4.547 8. 37.514 46.9345 1. 334.166 49.47 Chap. 미분방정식 11
Matlab 내장함수 ( 적응식 R-K 법 ) Ode3 Bogack와 Shampne, 1989 차와 3차 RK법사용 Ode45 Dormand와 Prnce, 199 4차와 5차 RK법사용가장널리사용 사용예 (ode45) >> [t, ] = ode45(odefun, tspan, ) tspan 정의방법 >> tspan = [t tf]; --- 구간정의 >> tspan = [t t1 tn]; --- 특정시간 Ode113 변동차수를갖는 Adams-Bashforth-Moulton 해법 엄격한오차공차나계산집약적상미분방정식을다루기에적합 3 Chap. 미분방정식 Matlab 내장함수 ( 적응식 R-K 법 ) Ex) P. 616 사례연구.6 ( 포식자 - 피식자모델 ) d1 d 1.1.61,.8.31, 1(), () 1 dt dt Step 1. 상미분방정식정의를위한 m-fle 작성 functon p = predpre(t,) p = [1.*(1).6*(1)*(); -.8*() +.3*(1)*()]; Step. 내장함수실행및그래프확인 >> tspan = [ ]; >> = [, 1]; >> [t,] = ode45(@predpre, tspan, ); >> plot(t,) >> plot((:,1),(:,)) % 상태공간그림 4 Chap. 미분방정식 1