[Python]미분방정식 - 테일러 급수, 오일러의 방법, 룽게-쿠타 방법

미분방정식 - 초기값 문제

• 테일러 급수
• 오일러의 방법
• 룽게-쿠타 방법

초기값 문제 Initial-Value Problem(IVP)

• 미분 방정식의 일반해를 이용, 초기조건을 이용하여 특수해 구하는 과정

• 일반해 + 초기조건 = 특수해

미분방정식

• x’ = x + 1

• x’ = 6t - 1

일반해

• \(x(t) = ce^t - 1\)

• \(x(t) = 3t^2 - t + c\)

초기조건

• x(0) = 0

• x(1) = 6

특수해

• \(x = e^t - 1\)
• \(x = 3t^2 - t + 4\)

오일러의 방법 실습

• 테일러 급수의 공식에서 앞에 2개의 항을 이용한 방법이다.

• \(x(t+h) = x(t) + hx’(t)\)

예시로 오일러 방법 사용해서 초기값이 x(1) = -4인 미분방정식 x' = 1 + x^2 + t^3에 대한

x(2)의 근삿값을 100단계를 적용하여 계산하라. 라고 할때? h=(2-1)/100 로 구할 수 있죠.

그리고 공식을 적용해주면?

• x(1) = -4

• x(1.01) = -3.82

• x(1.02) = …

• x(2) = …

이런식으로 구할수 있는 것이다.

다른 예시로 실습

# 초기값 : x'(t) = x, x(0) = 1, x(t) = e^t

import matplotlib.pyplot as plt
import numpy as np

def f(t, x): # x' -> x'(t) 구하는 함수
    return x # 초기값 보면 x

def Euler(a,b,x,n): # x는 x(0)
    h = (b-a)/n
    t = a
    answerT = [t]
    answerX = [x]
    for k in range(0,n):
        x = x+h*f(t,x) # x update -> x(t) 구함
        t = t+h # t update
        answerT.append(t)
        answerX.append(x)
    return answerT, answerX

# 초기값 : x'(t) = x, x(0) = 1, x(t) = e^t
# n = 10단계로, [0,4] 구간

n = 10
start = 0.0
end = 4.0

bX = []
bY = []
for i in range(0,n+1): # 기존 e^t 함수 그래프를 보여주기 위함
    x = start + (end-start)/n*i
    y = np.exp(x)
    bX.append(x)
    bY.append(y)

aT, aX = Euler(start, end, 1, n)

plt.plot(aT,aX, 'ro-', label='Euler\'s method')
plt.plot(bX,bY,'bo-',label='y=exp(x)')
plt.legend()
plt.show()

4차 룽게-쿠타법 RK4 실습

• 오일러랑 동일한데 미분 안하고 미분방정식을 해결한다는게 큰 특징이다.

• \(x(t+h) = x(t) + \frac{1}{6}(K_1 + 2K_2 + 2K_3 + K_4)\)

• K공식은 생략.. 아래 함수에 구현되어 있으므로 참고

# 초기값 : x'(t) = x, x(0) = 1, x(t) = e^t

import matplotlib.pyplot as plt
import numpy as np

def f(t, x): # x' -> x'(t) 구하는 함수
    return x # 초기값 보면 x

def RK4(a,b,x,n): # x는 x(0)
    h = (b-a)/n
    t = a
    answerT = [t]
    answerX = [x] # 여기까진 오일러랑 동일
    for k in range(0,n):
        k1 = h*f(t,x)
        k2 = h*f(t+0.5*h, x+0.5*k1)
        k3 = h*f(t+0.5*h, x+0.5*k2)
        k4 = h*f(t+h, x+k3)
        x = x+(k1+2.0*k2+2.0*k3+k4)/6.0
        t = t+h
        answerT.append(t)
        answerX.append(x)
    return answerT, answerX

전부 동시에 비교

# 전부 동시에 그려보겠음.
# 초기값 : x'(t) = x, x(0) = 1, x(t) = e^t
# n = 10단계로, [0,4] 구간

n = 10
start = 0.0
end = 4.0

bX = []
bY = []
for i in range(0,n+1): # 기존 e^t 함수 그래프를 보여주기 위함
    x = start + (end-start)/n*i
    y = np.exp(x)
    bX.append(x)
    bY.append(y)

aT, aX = Euler(start, end, 1, n)

rT, rX = RK4(start, end, 1, n) # 이것만 추가

plt.plot(aT,aX, 'ro-', label='Euler\'s method')
plt.plot(rT,rX, 'yo-', label='RK4')
plt.plot(bX,bY,'bo',label='y=exp(x)') # 구분 위해 점만 찍음
plt.legend()
plt.show()