Julia에서 Euler 방법을 시도했습니다.

1. Euler법이란?



Wikipedia에 따르면,

오일러법(오일러 쪽, 영국: Euler method)이란, 상미분 방정식의 수치 해법의 하나이다. 이 방법은 수학적으로 이해하기 쉽고 프로그래밍 방식으로 간단하기 때문에 수치 해석의 초보적 학습 문제로 잘 다루어집니다.
그러나, 1계단수 상미분 방정식의 수치 해법으로서는 오차가 축적되기 때문에, 정밀도가 나쁘고, 원래의 미분 방정식에 따라서는 어떠한 h 를 매우 원래 방정식의 해에 수렴하지 않는 것도 있는 방법이므로, 학습 목적 이외에 별로 사용되지 않는다.

얽히게 Euler법을 쓰면 미분의 정의를 이용해,
$$
y' =\frac{y(t +\Delta t) - y(t)}{\Delta t}
$$
이것을 식 변형하여 정리하면,
$$
y(t +\Delta t) = y(t) + y'\Delta t
$$
된다. 이것을 프로그램을 작성하기 위해 배열로 작성하면
y[n+1] = y[n] + y'*dt 
#y'は関数により、Δt = dt と表記する。

오차는 생기지만, 수치 계산의 초보적인 연습에는 적합하다고 한다. 최근 저는 Julia에서의 수치 계산을 연습하기도 하고 있으므로 사용해 보자. (Pkg로 ODEProblem 라든지 있지만, 이번에는 사용하지 않아)

보다 정확하게 수치 계산 한다면 개량 Euler법이라든지 Runge-Kutta법으로 해결하는 것이 좋을까 생각합니다.

2. 코드를 작성하기 위한 준비



이하의 순서에 따라 준비를 진행해 나가자.

2.1 이번에 사용하는 1층 선형 미분 방정식



여기서는
$$
y' = 3y + 2~~ , ~~y(0) = y_0
$$
를 예로 한다. 이것을 Euler법으로 생각해 간다.

2.2 Euler법에의 떨어뜨림


  • 에서 도입 한 식에 대입하면,

    $$
    y(t +\Delta t) = y(t) + (3y + 2)\Delta t
    $$

    자명하게 $y ​​= y(t)$ 이므로, $y(t)$ 로 오면,

    $$
    y(t +\Delta t) = (1 + 3\Delta t)y(t) + 2\Delta t
    $$

    이다. 이 식을 배열로 다시 작성하면
    y[n+1] = (1 + 3*dt)*y[n] + 2*dt   
    #Δt = dt
    

    이다. 여기까지 준비할 수 있다면 거의 완성된 것일 것이다.

    2.3 엄밀해도 준비해 두자



    수치해와의 비교도 하고 싶으므로, 엄밀해도 동시에 플롯하자. 주어진 미분 방정식을 미분 적분학의 지식을 사용하여 풀면,

    $$
    y(t) =\left(y_0 +\frac{2}{3}\right)e^{3t} -\frac{2}{3}
    $$
    된다.

    3. 코드 작성



    그럼, 지금까지 준비해 온 것을 사용해 코드를 기술해 가자.

    euler.jl
    using Plots;gr()
    
    # 点数
    N = 1000
    
    # 刻み
    dt = 0.001
    
    # 初期化
    y = zeros(N)
    
    # 初期値(y(0) = y_0 = 1 とした)
    y[1] = 1
    
    # 数値計算
    for n = 1:N-1
        y[n+1] = (1 + 3*dt)y[n] + 2*dt
    end
    
    # x = 開始値:間隔:終了値
    x = 0:dt:(N-1)dt
    
    # 厳密解
    function f(x)
        5//3 * exp.(3*x) - 2//3
    end
    
    #色々と指定してplot.
    #plot!はグラフを追加する機能をもつ。
    plot(x, f, color=:black, linewidth=3, label="Exact",xlabel="t",ylabel="y" , fmt = :png)
    plot!(x, y, color=:red, linewidth=2, linestyle=:dash, label="Numerical" , fmt = :png)
    

    이것을 실행하면, 제대로 그래프가 출력된다.


    단지 $t\in [0,1]$라고 오차는 없을 것 같다. 출력을 확인하면 pdf로 저장하자(png라면 화질이 거칠었다).

    euler.jl
    savefig("myplot.pdf")
    

    에서 실행중인 디렉토리에 저장됩니다. (이런 느낌)


    4. 마지막으로



    그런데, 여러분의 환경하에서도 실행할 수 있었습니까? 역시 코드가 움직이면 즐겁네요!Julia 에는 미분 방정식을 풀기 위한 Pkg도 있으므로, 그쪽에서 풀어준 ver도 앞으로 써 가려고 합니다.

    (추기)
    썼다 → 여기
  • 좋은 웹페이지 즐겨찾기