ODE + MATLAB

1 오일러 : first-order

2 Heun's method

3 Midpoint method

4 Runge-Kutta methods(Seconds, Third, Fourth, Higher)

---

 

 

 

1. F = ma  뉴턴의 제 2법칙 (힘의 합력 0이 아닐때)

 

2. 푸리에 열 방정식 : 열이 운동 특성, 열류 밀도 q 구할때, T는 온도의 미분: 높은곳에서 낮은곳 기울기 -가 된다. 단위 K.

 

여기서 q¯는 열류 밀도(heat flow density)[단위: W/㎡], κ는 열 전도도(heat conductivity: 물질의 고유한 특성)[단위: W/K/m], T는 온도(temperature)[단위: K]이다. 식 (1)을 잘 이해하려면 구배 연산자(gradient operator) ∇¯를 봐야 한다. 구배는 사실 미분(differentiation)이므로, 낮은 값에서 높은 값으로 가는 기울기는 (+)이다. 온도 T의 경우는 항상 높은 온도에서 낮은 온도로 변화가 일어나므로, 기울기 관점에서는 (-)가 된다. 높은 온도에서 낮은 온도로 가는 방향[기울기 (-), 구배 연산자의 값도 (-)]을 기준 방향 (+)로 만들기 위해 식 (1)의 구배 연산자 앞에 (-)를 붙여 열류 밀도 q¯가 (+)가 되도록 만든다. 식 (1)은 제안자인 푸리에 이름을 따서 푸리에의 열 전도 법칙(Fourier's law of heat conduction)이라 부른다. - 나무위키

 

추가자료: ghebook.blogspot.com/2012/07/heat-becomes-communication-heat-equation.html

열이 통신이 되다: 열 방정식(Heat Becomes Communication: Heat Equation)

3.픽의 확산 제 1법칙: 열역학에서 확산 과정을 나타내는 두 개의 법칙

Fick's law는 Fourier의 열전도법칙과 유사

 

* 관련 용어인 플럭스(Flux)

 

우선 기계공학 및 관련 분야에서는 금속을 용해할 때 금속이 산화하는 걸 막는 목적으로 금속 표면에 용제를 만들기도 하는데, 이를 플럭스라 부르며 염화물이나 플루오르화물, 수지 등이 주로 쓰인다. 이는 전자공학과 같은 분야에서도 기판에 땜질을 한 납의 유착과 융해를 막기 위해 플럭스를 사용하므로 통용되는 개념이다.

 

본문에서 말하는 플럭스는 위에 언급한 것과는 다른 개념으로, 사전적인 의미로는 주어진 방향에 대하여 수직인 단위 면적을 통하여 특정 물리량이 수송되는 비율을 말한다. 이 단위 면적은 해당 물리량이 '통과하거나', '넘어서는' 개념으로, 예컨데 흐르는 강이 있을 때 물이 흐르는 어느 한 지점에 횡단면을 만들었을 때 이 단면을 통과하는 물의 유량이 플럭스가 될 수 있다.

 

이 플럭스의 개념을 좀 더 자세히 파고들면, 수학적 정의와 스칼라적 개념과 벡터 필드 하에서의 개념, 수송 플럭스(Transport Flux)의 대표적인 여덟 개의 형태(Momentum Flux, Heat Flux, Diffusion Flux, Volumetric Flux, Mass Flux, Radiative Flux, Energy Flux, Particle Flux)에 관한 내용을 알 수 있으나 이 부분은 뉴턴의 점성법칙에서부터 맥스웰 정의까지도 다뤄야 하므로 본문에서는 넘어가도록 하겠다.

MHDV MREP

 

픽의 확산 법칙에서 다루는 플럭스의 개념은 위의 플럭스의 여덟 형태를 나열한 것 중 세 번째인 Diffusion Flux를 말하는 것으로, 이는 농도 기울기에 따라 고농도의 영역(region of high concentration)에서 저농도의 영역(region of low concentration)으로 넘어가는 플럭스를 가정한다. 공간(spatial), 즉 1차원(one dimension)에서 이 유량은 아래와 같이 표현된다.

 

 

(식 1. 1차원 공간을 통과하는 유동)

 

위 식을 다시 x로 미분하게 되면 픽의 제2법칙을 얻게 된다.

 

J는 확산 플럭스(diffusion flux)를 말하며, 이는 단위 시간(per unit time) 당, 단위 면적을 통과하는 단위 물리량을 말한다. 이는 어떠한 성분의 확산속도를 뜻하며, (#/cm²/s)의 단위로 표현될 수 있다. 몰(mol, mole)로서는 m-² s-¹ 으로 표현된다.

 

D는 확산 계수(diffusion coefficient), 혹은 확산 상수(diffusivity, cm²/s)로서 비례상수다.

 

N은 성분의 밀도(concentration)로서 mol/m³으로 표현될 수 있다.

 

x는 확산 방향으로의 거리를 뜻한다. cm 혹은 m의 단위로 표현될 수 있다.

 

위의 1번 식에서 볼 수 있듯이, 유동은 농도 기울기와 직접적인 연관성을 갖는다. 성분의 확산은 고농도 영역에서 저농도 영역으로 행해지며, D는 비례 상수이므로 특정 외부 조건에 그 값이 달라진다. D는 확산하는 성분인 확산 입자(diffusing particle)의 속도의 제곱(squared velocity)에 비례하며, 이는 스톡스-아인슈타인 관계식(Stokes-Einstein relation)에 따라 유체의 점성과 온도, 그리고 각 입자의 크기에 영향을 받는다.

 

1차원을 넘어서 2차원 혹은 그 이상의 차원에서는 반드시 기울기 연산자인 ∇(del, gradient operator)을 써줘야 한다. 이 때 위 1번 식은 다음과 같이 표현된다.

 

 

(식 2. 2차원 이상의 공간을 통과하는 유동)

 

위 식에서의 J는 확산 플럭스 벡터(diffusion flux vector)를 말한다

---

ODE 해보기.

 

상미분방정식의 수치적 해법 (Numerical Solution of Ordinary Differential Equations)😊

출처: https://egg-money.tistory.com/51 [완숙의 블로그]

*LHS, RHS

연산자 함수를 만들 경우 함수의 매개변수를 적어줄때 변수 이름으로 lhs(Left Hand Side) 와 rhs(Right Hand Side) 를 많이 사용한다.

매개변수가 두개일 경우 첫번째 매개변수에 lhs, 두번째 매개변수는 rhs 로 사용하고

매개변수가 하나일 경우 rhs 로 대부분의 개발자들이 많이 쓰고있다.

여기서 궁금한 것 한가지

나라면 매개변수가 하나일땐 lhs 로 사용할것 같다는것이다.

왜냐하면 사람이 글을 쓸때 왼쪽부터 오른쪽으로 쓴다. 코드도 한줄만 보면 왼쪽에서 오른쪽으로 짠다. class CBase 이것만 봐도 시작은 맨 왼쪽의 c 이다.

그렇다면 첫번째 매개변수는 lhs 로 사용 하는게 더 맞는것이 아닐까 ? 왜 대부분의 사람들은 rhs 로 사용을 하고있을까 그 이유가 궁금해졌다.

간단하게 설명을 하면,

---

int iVal;

iVal = iVal + 10 ; // iVal += 10;

---

 

iVal + 10 은

iVal.operator + (10); 이런 형식이 된다.

여기서 lhs 는 iVal 자기 자신이고 rhs는 10이다.

그렇기 때문에 lhs를 따로 쓸필요가 없다는것.

lhs 를 사용하는경우

lhs 를 사용하는 경우는 클래스 외부의 함수로 10 + iVal; 과 같이 했을 경우 (iVal 은 class 로 봐야함)

10.operator+(iVal) 이런건 존재 하지 않기 때문에 10을 계산 가능한 객체로 바꿔줘야 하기때문에

operator +(lsh, rhs) 라는 일반함수를 써야된다.

-blog.daum.net/trts1004/12109179

LHS 변수 와 RHS 변수 란 무엇인가?

-뒷 내용은 무슨 말인지 모르겠다.

---

컴퓨터는 이산적인 값에서 작동하는 기계입니다.

우리의 정신세계에서 이상적인 선은 존재할 수 있지만 현실에서는 가산적인 점의 집합이 결국 선이겠죠. 

연속적인 세계에서 정의된 해석적 미분방정식(Analytical)의 풀이 방법은 컴퓨터에서 적용할 수 없습니다. 

그래서 우리는 이산적인 세계에서 미분방정식의 해를 검출할 수 있는 다른 방법이 필요합니다.

 

 

Analytical은 공식을 통해 , 풀이 방법은 익혔으므로 생략하겠다.

 

Numerical을 살펴보자.

-참조 : egg-money.tistory.com/51

상미분방정식의 수치적 해법 (Numerical Solution of Ordinary Differential Equations)

 

 

(c = 항력계수)
떨어지는 물체가 있을 때, 종단속도에 관한 미분 방정식입니다.
직관적으로 보면 LHS = 가속도(속도의 변화)는,
RHS = 떨어지는 중력의 가속도 - 항력계수 * 속도 로 이해할 수 있습니다.

이 때, LHS는 속도의 함수라고 볼 수 있겠다

.

.여기서 Slope는 다음과 같습니다.

(속도는 시간의 함수므로)

정리하면, 다음과 같습니다.

Numerical Method

ODE 근사함수가 기본적으로 업데이트 방식으로 값을 얻어내기 때문에 경계조건이 필요합니다.

이 부분은 해석적으로 풀더라도 해의 유일성을 판단하기 위해 필요한 부분입니다.

 

 

문득, (종횡비(가로세로비율): axis square로 맞추기)

 

 

같게 그려진다.

새로운 문제 2;

---

오일러 적분 :

 

---

오차 개념 : 

수치해석은

어느 방정식의 정확한 해(Exact Solution)을 구하는 것이 아닌

근사적인 해(Approximate Solution)을 구합니다.

정확한 해를 구하는 것은 앞서 말한 해석적 방법을 통해 구한 해입니다.

따라서 수치적인 방법을 이용해 얻은 해는 오차가 존재하는 근사적인 해가 됩니다.

 

오차에는 크게 두가지가 있습니다.

하나는 반올림 오차(Round-off Error)이고,

다른 하나는 절단 오차(Truncation Error)입니다.

 

반올림 오차(Round-off Error)는 말 그대로

어느 특정 자릿수 이하 버리거나 반올림 또는 올림함으로서 발생하는 오차입니다.

흔히 파이(pi,π)=3.14 로 반올림하거나 1/3=0.3333 으로 반올림하는 것이 대표적인 것입니다.

 

절단 오차(Truncation Error)는

수학적인 연산의 근사적인 표현으로 발생하는 오차인데... 말이 좀 어렵죠??

절단 오차를 이해하려면 테일러 급수(Taylor Series)를 이해해야 합니다.

테일러 급수는 나중에 다룰 예정입니다만, 간단히 설명하자면

복잡한 형태의 연속적(continuous)이고 부드러운(smooth) 미분가능한 함수 f(x)를

다항식(Polynomial)의 합으로 근사한 것입니다.

테일러 급수는 무한급수이기 때문에 어느 항(term)이상을 절단(truncation)함으로서

함수를 근사할 수 있는데, 이때 발생한 오차를 절단 오차라 부릅니다.

 

예)

맨 마지막 항 O(x)가 절단오차 입니다.

 

오차해석이 중요한 이유는 무엇보다

근사적으로 구한 해를 실제로 사용할 수 있는 해인가를 판단하기 위함 입니다.

예를 들어 어느 방정식의 정확한 해가 5가 나왔는데, 수치적으로 구한 근사 해가 3이라면,

40%의 오차가 발생하게 됩니다. 이 40%를 감안하고 근사 해를 인정할 것인지, 아니면 더 정확한 해를 구해야 하는지는

엔지니어의 임무

-참조 : blog.naver.com/mykepzzang/220061309981

[수치해석] 2. 오차해석, Analysis of Error

---

오차의 정량화(Quantification of Error )
오차를 계산하는 방법에 대해 알아보겠습니다.

 

1. 참오차(True Error)

 

2. 참 상대오차(Relative True Error)

 

위의 오차계산은 참값을 알아야하죠.

즉, 해석적 방밥의 정확한 해를 알아야 오차가 계산 가능합니다.

하지만, 수치해석의 문제에서는 참값은 모른다는게 문제죠...

그래서 근사오차개념을 도입합니다.

 

3. 근사오차(Approximate Error)

 

4. 근사 상대오차(Relative Approximate Error)

 

사실 수치해석법은 반복법(Iteration)을 이용합니다.

반복법은 간단히 얘기하면, 똑같은 계산을 계속 반복합니다.

예를 들어, 이전 근사값을 계산해서 현재 근사값을 얻으면,

그 현재 근사값으로 다음 근사값을 얻는 방법으로 계속 반복 계산합니다.

그리고 나서 우리가 어느 정도 허용한 오차 범위내에 우리가 기대한 데이터가 들어오면 계산을 마칩니다.

즉, 많은 반복 계산을 거쳐서 근사 상대오차가 어느 정도 0으로 수렴(허용 오차 범위내에 들어오면)하면,

계산을 마칩니다.

그럼 그 허용한 오차범위를

라 한다면,

를 만족해야 계산을 끝마치고 우리가 원하는 어느 값을 얻을 수 있겠죠.

---

훈의 법칙 과 midpoint 법칙:

< 사다리꼴 공식(trapezoidal rule of integration) > 

수치적분(numerical integration)을 사용하는 이유는 해석적으로 적분이 안되는 경우가 대부분입니다. 고등학교 교과과정 및 대학교 1학년때 배우는 미적분학에서 거의 대부분 해석적인 적분이 가능합니다만, 함수가 조금만 복잡해도 적분불가능한 경우가 대부분입니다. 그렇기 때문에 공학에서 수치적분을 많이 이용합니다. 수치해석 첫번째 포스팅에서도 설명했듯이, 수치해석의 목적은 어떤 복잡한 식의 정확한 값은 아니더라도 근사적으로 이 식의 값이 얼마인가를 알고 싶기 때문입니다.

 

수치적분의 방법도 다양합니다. 가장 단순한 방법부터 가장 많이 사용되는 방법까지 이번 수치적분 포스팅을 통해 차근차근 설명해 드리도록 하겠습니다.

 

이번에 소개해 드릴 수치적분 방법은 '사다리꼴 공식(trapezoidal rule of integration)'입니다.

사다리꼴 공식은 가장 단순하지만 오차가 큽니다.

 

 

대부분 아시다시피 적분기호 바로 옆 a를 '아랫끝(lower limit of integration)', b를 '윗끝(upper limit of integration)'이라고 하고 f(x)를 '피적분함수(integrand)', 그리고 dx를 '미분소(differential of the variable x)'라고 부릅니다.

고등학교 교과과정에서 적분의 의미는 어떤 함수 f(x)의 넓이라고 배웠습니다. 아래 그림을 보시면,

 

 

사다리꼴 공식은 f(x)의 적분을 윗 그림처럼 초록색영역의 사다리꼴 넓이로 근사시키는 방법입니다. 실제 적분과의 오차는 노란색 영역만큼이 되겠습니다.

사다리꼴 넓이를 구하는 공식은 초등학교때 배우죠?

초록색영역의 넓이를 구하려면,

 

입니다.

 

따라서 구하고자 하는 적분 값을 다음과 같이 표현할 수 있습니다.

 

여기서 우변의 두번째 항을 보시면, 함수값의 평균값이라는 것을 알 수 있습니다.

 

미적분학에서 1변수함수 f의 구간 [a, b]에서의 평균값은 다음과 같이 나타낼 수 있습니다.

 

우변의 (b-a)를 양변에 곱해주면 다음과 같습니다.

참조 : blog.naver.com/PostView.nhn?blogId=mykepzzang&logNo=220583458401

 

**

---

합성 사다리꼴 공식 : 

사다리꼴 공식은 다양한 수치적분 방법 중 가장 단순하기 때문에 직관적으로 이해할 수 있지만, 오차가 크다는 것이 가장 큰 단점입니다. 아래의 예를 한번 볼까요?

 

실제 적분값은 부분적분을 이용해서 구해보면(설마...부분적분 모르시는 분은 없겠죠?),

입니다. 그래프를 이용해서 나타내면 다음과 같습니다.

 

 

바로 이전 포스팅에서 소개한 사다리꼴 공식을 이용해서 한번 구해봅시다.

입니다. 그래프를 이용해서 나타내면 다음과 같습니다.

그래프를 통해 보면 알겠지만, 사다리꼴 공식을 이용한 수치적분 방법은 적분구간에 따라 오차의 변동이 매우 심하고 어느 경우에는 오차의 크기(그래프의 노란색 영역)가 수치적분값(그래프의 초록색 영역)보다 더 큰 경우도 있습니다. 이처럼 사다리꼴 공식은 적분구간에 따라 실제 적분값을 잘 반영하지 못하는 매우 치명적인 단점을 가지고 있는 것이죠.

 

그래서 등장하게 된 것이 '합성 사다리꼴 공식(multiple segment trapezoidal integral)'입니다.

구하고자 하는 적분영역을 여러개로 나눈 뒤 각각의 쪼개진 영역에서 사다리꼴 방법을 이용해 값을 구하고 나중에 전부 합치는 것이죠. 마치 '구분구적법(mensuration by parts)'과 매우 비슷한 개념입니다. 아래 그림을 보시죠.

 

구하고자 하는 적분구간을 여섯부분으로 쪼개서 각 부분의 사다리꼴을 이용해 값을 구하고 최종적으로 여섯개의 사다리꼴을 전부 더하는 것입니다. 오차는 그래프만 보더라도 확실히 줄어든 것을 볼 수 있습니다. 물론 계산횟수가 늘어나는 것이 좀 귀찮을 뿐이죠.

 

합성 사다리꼴 공식을 한번 유도해 보겠습니다.

 

 

윗 식의 우변을 h/2 로 묶고 정리하면 다음과 같습니다.

 

h=(b-a)/n 이므로, 정리하면

 

입니다. 이렇게 공식유도가 되었습니다. 여기서 적분구간을 나눈 갯수, 즉 n의 값이 커질수록 적분값은 더 정확해집니다.

 

여기서 눈 여겨볼 것이 우변에서

 

이라는 것입니다. 즉, 그래프에서 구하고자하는 적분구간의 평균 높이(=함수값) 이라는 것이죠. 이게 왜 함수값의 평균이 되느냐하면,

분모에서 f(a)값 1개, f(b)값 1개, 그리고 시그마 부분에서 2(n-1)개 이므로, 분모의 항의 개수는 전부 2n개 입니다. 이 함수값의 합을 2n으로 나눠주니까 결국 함수값의 평균이 되는 것이죠. 따라서 다음과 같이 나타낼 수 있습니다.

 

 

한 가지 알아두셔야 할 사항은 사다리꼴 공식은 위에 그래프를 통해 알 수 있듯이, 주어진 함수 f(x)를 1차 다항식(linear polynomial)으로 근사시켜서 적분값을 구하는 방법

---

훈법칙 , mid point 정리:

좀 더 정확한 해를 얻기 위한 논리는 미분항 (dy/dx)을 개선하는 것입니다. i번째의 x와 (i+1)번째의 x사이의 기울기의 평균(혹은 대표값)을 사용하면 오차가 줄어든다는 것이죠.

 

 

오일러 방법과 마찬가지로 (i+1)번째의 y값을 추정하기 위해서

 

오일러 방법에서는 윗 식을 이용해서 바로 (i+1)번째의 y를 구하면 끝이지만, 호인 방법에서는 여기서 끝나지 않습니다. 호인 방법에서는 이 식을 다음 y값의 기울기를 예상하는 용도로만 사용합니다. 그래서 (i+1)번째 y에 상첨자 0을 붙였습니다. 호인 방법에서 이 식을 '예측 방정식(predictor equation)'이라고 합니다.

 

이 예측 변수 방정식을 이용해서 우선 다음을 구합니다.

 

그럼 이제 우리는 두가지의 기울기를 가지고 있습니다. i번째에서의 기울기와 (i+1)번째에서의 기울기를 가진 것이죠.

이 두 가지를 이용해서 평균기울기를 구합니다.

 

이 평균 기울기를 이용해 최종적으로 (i+1)번째 y값을 구할 수 있습니다.

 

호인 방법에서는 윗 식을 '수정 방정식(corrector equation)'이라고 합니다.

 

윗 과정을 그림으로 나타내면 다음과 같습니다.

 

자료출처 : Numerical Methods for Engineers 6th edition / Chapra, Canale / McGRAW-HILL

 

호인 방법을 이용해서 간단한 문제 하나 풀어보겠습니다.

 

x=2일 때 y의 값을 구하는 것이 문제입니다.

 

우선 예측 방정식을 이용하여,

이제 수정 방정식을 이용하면

 

아래 그림은 실제 해와 오일러 방법, 호인 방법을 비교한 그림입니다.

 

자료출처 : Numerical Methods for Engineers 6th edition / Chapra, Canale / McGRAW-HILL

 

중앙점 방법: midpoint method

 

오차가 큰 오일러 방법을 개선시킨 또 다른 방법은 '중앙점 방법(midpoint method)'입니다.

앞의 오일러 방법과 호인 방법을 잘 이해하셨다면, 중앙정법도 별 무리 없이 이해하실거라 생각됩니다.

아래 미분방정식을 생각해볼까요.

 

 

위 미분방정식에서 다음 (i+1)번째의 y값을 추정해야 하는데, 우선 오일러 방법을 이용해 주어진 h의 절반인 h/2에서의 기울기를 구해야 합니다. 즉 i+(1/2)번째의 y값을 먼저 구하는 것이죠.

 

여기서 구한 i+(1/2)번째 기울기를 통해서 i+1번째의 y를 추정할 수 있습니다.

이게 끝입니다. 개인적으로 매우 간단한 방법이라고 생각됩니다.

 

이 과정을 그림으로 나타내면 다음과 같습니다.

 

자료출처 : Numerical Methods for Engineers 6th edition / Chapra, Canale / McGRAW-HILL

 

위 그림 (a)에서 i+(1/2)의 기울기를 구하고, 그 기울기를 다시 i번째에 적용해서 (i+1)번째 y를 추정하는데 사용되고 있습니다.

 

문제 하나 풀어보겠습니다.

우선 초기값이 x=1에서 y=2의 값이 주어졌고, h의 크기는 1이므로, 초기값에서 h=0.5만큼 증가한 x=1.5에서의 기울기를 구해야 합니다.

 

이제 y(2)값을 구하기 위해 중앙점에서 구한 기울기 3을 대입합니다.

---

Runge-kutta methods: 

차례로 2차, 3차, 4차 방법 논의

 

2차 :

상미분방정식 풀이법은 '2차 룽게-쿠타(Runge-Kutta, 룬게-쿠타 또는 런지-쿠타(영미식 발음)) 법'입니다.

룽게-쿠타 법은 1900년경 독일의 수학자 룽게(C. Runge)와 쿠타(M. W. Kutta)의 이름을 따서 명명된 상미분방정식 풀이법입니다. 룽게-쿠타 법의 장점은 고차미분항 없이 상당히 정확한 해를 얻는다는 것 입니다.

 

우선 다음과 같은 상미분방정식이 주어졌다고 해봅시다.

 

 

(i+1)번째 y를 추정하기 위해 i번째 y의 값과 테일러 급수를 이용합니다. 1차 테일러 급수를 이용한 것이 바로 오일러 방법이죠.(오차 생략)

 

이제 1차 테일러 급수를 2차까지 확장해 보겠습니다.(오차 생략)

 

여기서 f'(x, y)는 연쇄법칙에 의해 다음과 같이 표현할 수 있습니다.

 

 

이제 이것들을 조합하면 다음과 같이 되겠습니다.

 

 

이 포스팅 초반에 룽게쿠타법은 고차미분항이 없다고 언급했습니다.

이제부터 전개될 내용이 룽게-쿠타 법의 핵심입니다. 미분항을 다음과 같이 가정합니다.

 

그러면 우리가 알고 있던 오일러 방법의 식이 다음과 같이 됩니다.

 

(a₁k₁+ a₂k₂)이 부분을 수치해석에서는 증분함수(increment function)라고 부릅니다.

오늘 다루는 주제가 '2차(2nd-order)'인 이유도 바로 두 개의 항으로 이루어진 a와 k 이기 때문이고, 또 테일러 급수 2차 확장까지 다루기 때문입니다. 여기서 a는 상수이고, k는 다음과 같습니다. 

 

k1은 dy/dx와 같지만(dy/dx=f(x, y)이므로), k2는 k1의 함수라는 것을 알 수 있습니다.

 

그럼 이 k2를 어떻게 처리하느냐가 관건인데, k2는 이변수 테일러급수를 이용합니다. 이변수에 대한 테일러급수는 미적분학 교과서를 참고하세요. 이변수에 대한 테일러 급수의 형태는 다음과 같습니다.

 

이제 k2를 이변수 테일러 급수를 이용해 1차까지 전개하면,(오차 생략)

 

윗 식을 빨간색으로 표시한 식에 대입합니다.

이 식을 정리하면,

 

이제 1식과 2식의 항과 항끼리 비교를 통해 다음을 얻을 수 있습니다.(이 과정은 직접 해보시기 바랍니다.)

 

 

여기서 a2를 1/2 라고 가정하고, p1, q11을 1이라고 두면, 이것은 호인 방법(Heun method)가 됩니다.

또 a2를 1로 가정하고, a1을 0 그리고 p1, q11을 1/2이라고 두면, 이것은 중앙점 방법(Midpoint method)가 됩니다.

이것을 표로 정리했습니다.

 

 

---

3차 룽게쿠타:

 

---

4차 룽게쿠타 :

 

룽게-쿠타 방법 중 가장 많이 사용하고 또 수치해석에서의 상미분방정식 초기값 문제를 통틀어서 가장 많이 사용되는 4차 룽게-쿠타 법을 알아보겠습니다. 일반적으로 룽게-쿠타 법이라고하면 4차 룽게-쿠타 법을 얘기합니다.

 

4차 룽게-쿠타 법을 직접 유도하지는 않고 식을 소개해드리는 정도로 이야기하려고 합니다.

다음과 같은 초기값형태의 미분방정식이 있다고 해봅시다.

 

지금까지 해왔던 것 처럼 초기값문제의 해를 추정하기 위한 iteration update 공식은 테일러 급수 확장을 이용합니다.(오차 생략)

 

여기서 dy/dx=f(x,y)를 이용하면,

 

이제 이 식을 증분함수(increment function)으로 나타내면, 다음과 같이 표현할 수 있습니다.

 

여기서 2차 룽게-쿠타 법과 마찬가지로 각 a와 k를 어떻게 구성하느냐에 따라 달라집니다.

 

룽게가 제안한 가장 흔한 4차 룽게-쿠타 법(classical fourth-order RK method)은 다음과 같습니다.

윗 식의 개념을 그림으로 묘사하면 다음과 같습니다.

 

자료출처 : Numerical Methods for Engineers 6th edition / Chapra, Canale / McGRAW-HILL

 

쿠타가 제안한 다른 종류의 4차 룽게-쿠타 법은 다음과 같습니다.

이 식은 마치 수치적분의 3/8 심프슨 공식과 비슷하게 4차 룽게-쿠타 3/8 공식이라고 부르기도 합니다.

---

**고차 룽게쿠타 법: 

---

진자의 운동/ 자유도 / 선형성 개념:

 

우선 운동방정식을 세우기 전에 '자유도(Degree of Freedom, DOF)'의 개념을 알아보겠습니다.

역학(力學)에서 말하는 자유도란 물체의 움직임을 표현하기 위해 도입한 개념으로 물체가 이동하는 방향의 개수와 회전할 수 있는 방향의 개수의 합을 나타냅니다. 동역학에서 자유도는 매우 중요하죠. 만약 어떤 물체가 특정한 한 방향으로만 쉽게 움직인다면 그 물체는 특정한 방향에 대해 '1자유도'를 가진다고 합니다. 예를 들어 엘리베이터의 운동은 1자유도로 가정할 수 있는것이죠. 부피를 가진 강체(rigid body)는 3차원 공간에서 병진운동(translation)과 회전운동(rotation)을 표현해야하므로 6자유도를 갖습니다.

 

 

 

3차원 공간에서 질점은 부피를 가지지 않기 때문에 회전운동이 없으므로 병진운동에 대한 3자유도를 갖습니다. 만약 2차원의 평면운동이라고 하면 강체의 경우에는 x, y축에 대한 병진운동과 z축에 대한 회전운동을 하기 때문에 3자유도를 갖고, 질점의 경우에는 x, y축에 대한 병진운동만 하므로 2자유도를 갖습니다. 1차원 직선운동에 대해서는 강체와 질점 둘 다 병진운동만 하므로 1자유도만 갖습니다.

 

이제 운동방정식에 대한 이야기를 해볼까요. 결론부터 얘기하자면 운동방정식은 결국 뉴튼의 제 2법칙인 "F=ma"입니다. 굳이 이 식을 해석하면 '질량에 작용하는 외력(force)의 합은 질량과 가속도의 곱과 같고, 이때 가속도의 방향은 외력이 작용하는 방향과 같다.'는 것이죠. 운동방정식은 어떤 시스템의 거동을 수학적으로 묘사하기 위한 식이라고 보시면 될 것 같습니다. 이처럼 어떤 자연현상 또는 물리적현상을 수학적으로 표현하는 것을 가리켜 '모델링'이라고 합니다.

 

우선 아래 그림을 통해 단순진자를 한 번 살펴보겠습니다.

 

 

(출처 : Wikipedia)

 

제가 진동이라는 학문을 처음 소개할 때, 진동의 물리적 의미는 위치에너지와 운동에너지의 상호관계라고 표현했습니다. 위 그림에서 진자가 움직임을 보면 진동의 물리적 의미가 딱 맞아떨어지는 것을 쉽게 알 수 있죠. θ가 최대가 될 때 질량 m을 가지는 질점의 위치에너지도 최대가 되고 운동에너지는 최소가 됩니다. 반면에 θ가 0일 때 위치에너지는 최소가 되지만 운동에너지는 최대가 되죠. 이렇게 위치에너지와 운동에너지의 상호작용이 바로 진동입니다. 윗 그림에서 단순진자는 질점이 θ에 대해서만 움직이므로 1자유도를 갖습니다.

 

학부 수준에서 진동은 '선형(linear)'라는 조건이 붙는다고 했습니다. 선형은 다음과 같은 성질을 가지고 있습니다.

 

사실 진자가 선형이라는 조건을 만족하려면, 위의 그림에서 θ는 매우 작아야합니다. 진자가 왕복운동하는게 보이지 않을 정도로요. 따라서 실제 상황에서 선형이라는 조건을 적용하기에는 무리가 있지만, 해석을 굉장히 단순화하기 때문에 '선형'이라는 가정을 많이 사용하는 것이죠. 우리가 살고있는 이 세상에 선형시스템을 만족하는게 얼마나 있을까요? 저는 지구의 90%이상은 '비선형(nonlinear)'시스템이라고 생각합니다. 그만큼 선형이라는 가정은 실제와 조금 동떨어져 있지만, 우리는 진동을 처음 배우는 입장이기 때문에 해석을 단순화하고자 선형이라는 가정을 받아들이도록 합시다. 선형이라는 가정이 가져다 주는 장점은 바로 '중첩의 원리(superposition principle)'입니다. 미분적분학이나 미분방정식에서도 나오지만 중첩의 원리는 나중에 따로 설명하도록 하겠습니다. 참고로 비선형 진동해석은 '푸리에 급수(fourier series)'를 이용해서 표현합니다.

 

그럼 진자가 '선형'이라는 가정을 전제로 운동방정식을 한 번 유도해봅시다.

 

 

여기에서 m은 부피는 없고 질량m인 질점을 나타내고, 원점O와 질점m을 연결하는 줄의 길이가 l이고 줄의 질량은 고려하지 않는다고 가정합시다. 그리고 g는 중력가속도를 나타냅니다.

 

그럼 진자에 대한 운동방정식을 구해보겠습니다. 핵심은 원점O에 대한 모멘트 식을 세우는 것입니다.(모멘트에 대한 내용은 정역학(Statics) 교과서를 참고하시기 바랍니다.) 원점O에서 반시계방향을 +로 두면,

 

 

 

위에 빨간색으로 표시한 식이 바로 단순진자의 운동방정식입니다. 즉 진자의 움직임을 표현하는 식이죠. 식을 유도하는 도중에 매우 작은 θ(very small θ)라고 했죠. 바로 '선형'이라는 가정때문에 가능한 것입니다. 그럼 왜 sinθ가 매우 작으면 θ가 되는지는 테일러급수(정확히는 매클로린 급수)를 통해 알 수 있습니다.

여기에서 θ는 매우 작기 때문에, θ을 거듭제곱한 항들은 거의 0에 가까워서 고려를 하지 않습니다.

그리고 단순진자의 운동방정식을 보면 질점의 질량 m과는 무관하다는 것을 알 수 있습니다.

참조: blog.naver.com/PostView.nhn?blogId=mykepzzang&logNo=220627110813&parentCategoryNo=&categoryNo=&viewDate=&isShowPopularPosts=false&from=postView

---

Adaptive RK method:

---

요약하자면.

 

 

---

절반까지 왔다!!

 

---