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

Numerical Solution of Ordinary Differential Equations

Bello-
암호를 이해하기 위해서는 ODE(나중에 포스팅 하겠ㅅ..)가 필요합니다.

Intuition Concept

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

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

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

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

ODE의 풀이는 생략하고 미분방정식의 꼴에서 부터 생각해보아요.

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

이 때, LHS는 속도의 함수라고 볼 수 있겠네요.

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

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

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

---

Function

적분의 근사방법이 굉장히 많았듯이 ODE도 근본적으로 근사하는 방법이기에 다양한 방법이 있습니다.

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

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

Solve the ODE
[t y] = solver_name(ODEfunc,구할t의범위,y0)
(where y0 at t0)

---

Example1

1 < t < 3 의 범위내에서

를 풀어라. (B.C y = 4.2 at t = 1)

위의 원리를 만족하는 함수를 만들어보자.

function[t, y] = ode_self(n, initial, final, f, y0)
dt = (final-initial)/n;
t(1) = initial;
y(1) = y0;
    for i = 1: n
    t(i+1) = t(i) + dt;
    y(i+1) = y(i) + dt * f(t(i),y(i));
    end
end

그리고 실행해보자.

실행코드
dydt = @(t,y)(t^3 -2*y)/t
[t,y] = ode_self(100, 1, 3, dydt, 4.2)

plot(t,y)

ODE45로도 풀어보자.

dydt = @(t,y)(t^3 -2*y)/t
[t,y] = ode45(dydt, [1:0.01:3], 4.2)

plot(t,y)

같게 그려지는 것을 알 수 있다.

Example2

(MATLAB_5th_Edition_An_Introduction)

x와 v와의 관계를 원하므로 약간의 조작을 가했다.

K = 30
m = 1500
v0 = 90*1000/3600
% 단위 변환
dvdx = @(x,v)(-1)*(K/m)*(v^2)*(x+1)^3
range_x = [0:0.001:3]
[x, v] = ode45(dvdx,range_x,v0)

plot(x,v)
xlabel('x')
ylabel('v')

그럼
Poopaye-