본문바로가기
[코딩수학] 코딩으로 '미분방정식' 정복하기
수학동아 2019.06.06 13:48

 

코딩 수학 6번째 시간입니다. 다들 코딩을 이용해 수학 문제를 효율적으로 풀고 계신가요? 간단한 수학 문제는 손으로도 풀고 코딩으로도 풀 수 있지만, 실생활 문제나 어렵고 복잡한 함수가 포함된 수학 문제는 손으로 풀기 어려운 경우가 많아요. 코딩을 활용하면 쉽게 해결할 수 있으니 열심히 익히고 활용해 보자고요! 

코딩 명령어를 익힌 뒤엔 해당 코딩 명령어로 풀 수 있는 수학 문제를 찾아 댓글로 달고 친구들과 토론해 보세요.

 

 

 

 

침몰한 배를 찾아라! 미분방정식과 벡터장

 

 

이번에 코딩으로 정복해 볼 수학 개념은 '미분방정식'입니다. 미분방정식은 비행기기의 속도 파도의 움직임처럼 시시각각 변하는 대상과 그 대상의 변화율의 관계를 방정식으로 나타낸 거예요.

 

예를 들어 단위 시간 동안 붕괴하는 방사성 원자핵의 수, 즉 시간에 대한 원자핵 수의 변화율이 특정 시간 \large t에서의 원자핵 수를 나타내는 함수 \large N(t)에 비례하면, 미분방정식은 \large \frac{dN(t)}{dt}=kN(t)가 됩니다. 여기서 \large \frac{dN(t)}{dt}가 시간에 대한 원자핵 수의 변화율이고, \large k는 상수예요. 미분방정식을 푼다는 건 함수 \large N(t)를 알아내는 겁니다.

 

이렇게 미분방정식의 해는 대부분 수가 아닌 함수 같은 '수식'으로 나온다는 게 특징인데요. 이렇게 나온 수식에 특정 값을 대입할 때마다 원하는 값을 얻을 수 있습니다. 위 미분방정식에서 얻은 함수 \large N(t)에 \large t 대신 특정 시간을 대입하면 그 시간에 남아있는 원자핵의 수를 알 수 있는 것처럼요.

 

 

 

아이작 뉴턴(1642~1727)과 고트프리트 라이프니츠(1646~1716).

 

미분방정식을 만들고 해를 찾는 방법을 처음 제시한 사람은 영국 수학자 아이작 뉴턴과 독일 수학자 라이프니츠다. 특히 라이프니츠는 ‘미분방정식’이라는 용어를 처음 사용하고 \large dx, dy, \frac{dx}{dt}(변화율을 나타내는 기호)와 같이 현재 미분적분학에서 사용하는 기호들을 소개했다.

 

 

미분방정식은 현재 시점에서의 변화를 분석하면 가까운 미래를 정확하게 예측할 수 있어 날씨, 주가 예측 등 활용도가 무궁무진한데요, 문제는 미분방정식의 해를 구하는 게 무척 어렵다는 겁니다. 미국 클레이 수학 연구소가 뽑은 7대 난제 중 하나인 ‘나비에-스토크스 방정식’은 공기와 물 같은 유체의 움직임을 나타낸 미분방정식인데 해를 구하는 건 둘째 치고 유일한 해가 있는지 조차 아직 알 수 없지요.

 

복잡한 함수를 그래프로 그리면 함수의 성질을 쉽게 파악할 수 있죠? 비슷하게 미분방정식에서는 그래프 대신 ‘벡터장’, 다른 이름으로는 ‘방향장’을 그립니다. 벡터장은 유클리드 공간의 각 점에 벡터_힘, 속도처럼 크기와 방향을 가진 물리량를 대응시킨 건데요, 미분방정식을 알고 있으면 각 점에서 물체가 움직이는 방향과 그 크기를 벡터로 표현할 수 있습니다.

그러면 코딩을 이용해 미분방정식을 벡터장으로 나타내고 일반해를 구해볼까요?

 

 

 

 

.

.

.

 

 

 

 

코 딩 타 임

 

 

아래 상황을 읽고 미분방정식을 벡터장으로 나타내는 코딩 명령어를 이용해 침몰한 배의 위치를 추측해보자.

 

 

좌표평면 위의 한 지점 (50, 80)에서 침몰하는 배로부터 구조 신호가 왔다. 구조선이 신호가 온 지점 근처까지 가려면 약 1시간이 걸린다. 그런데 통신이 두절돼 1시간 뒤의 위치는 알 수 없다. 다행히 구조대는 근처 해역의 해류의 움직임을 나타내는 미분방정식을 알고 있었다!

 

\LARGE \frac{dx(t)}{dt}=-2x(t)+2y(t)         ................        \LARGE \frac{dy(t)}{dt}=2x(t)-5y(t)

 

 

 

 

0단계 기초 익히기

 

벡터장 그리기

var(‘변수’)는 변수를 정의하는 명령어로, ‘’ 사이에 변수를 띄어쓰기로 구분해 지정한다. 예를 들어 좌표평면은 var(‘x y’)이다.  

plot_vector_field((미분방정식), (x축 범위), (y축 범위), aspect_ratio=1)은 벡터장을 그리는 명령어로, 첫 번째 빈 괄호에 미분방정식을 쉼표로 구분해 입력한다.

두 번째와 세 번째 괄호에는 좌표평면의 표시 범위를 쉼표로 구분해 입력한다. 예를 들어 x, y축 모두 0부터 90까지 그리고 싶으면 (x, 0, 90), (y, 0, 90)이라고 쓴다.  

❸ aspect_ratio는 좌표평면의 가로와 세로의 비율로, 가로 길이를 세로 길이로 나눈 값이다.

 

일반해 구하기

❶ var(‘t’)를 입력해 t를 시간에 관한 변수로 지정한다.

 

x=function(‘x’)(t)는 x를 변수 t에 관한 함수로 지정하는 명령어다. 다음 줄에 x를 y로 바꿔 y를 t에 관한 함수로 지정하자.

diff(함수)==미분방정식은 미분방정식을 정의하는 명령어다. 괄호 안에는 함수의 분자, 분모를 쉼표로 구분해 입력하고, ==뒤에 식을 넣는다. 첫 번째 식의 이름을 de1이라고 정했다면 de1=diff(x, t)==-2*x+2*y라고 입력한다.

❹ desolve_system([방정식], [변수])은 미분방정식을 푸는 명령어로 일반해를 구하자. 첫 번째 중괄호에는 ❸에서 입력한 미분방정식을 쉼표로 구분해 입력하고, 두 번째 중괄호에는 변수를 넣는다.   

 

 

 

1단계 따라하기 _벡터장 그리기

 

①벡터장을 그리는 명령어를 이용해 해류의 움직임을 그림으로 나타내자.

 

②벡터장에서 배의 처음 위치 (50, 80)를 표시한 뒤 해류의 방향을 고려해 배가 어떻게 떠밀려갈지 빨간색 펜으로 표시해보자. 미분방정식의 일반해를 구하면 특정 시간에 배의 위치를 정확하게 구할 수 있다.

 

 

 

 

 

2단계 나아가기 _일반해 구하기

 

❶ var('t')를 입력해 시간 t를 변수로 지정한 뒤 function('x')(t)를 입력해 xt에 관한 함수로 지정하고 x라고 이름 붙이자. 마찬가지로 yt에 관한 함수로 지정하자.

❷ diff(x, t)== 뒤에 우리가 알고 있는 시간에 따른 x의 변화율을 입력하고 de1이라고 이름붙이자. 시간에 따른 y의 변화율은 x대신 y를 적으면 된다.

❸ desolve_system([de1, de2], [x, y])를 입력하면 일반해 x(t), y(t)를 알 수 있다. 일반해를 항상 구할 수 있는 건 아님을 유의하자. 일반해의 t에 1을 대입하면 1시간 후의 배의 위치를 알 수 있다.

 

 

 

 

3단계 마무리하기 _동시에 해보기

 

 

 

도전 문제 해역에서 해류의 x축의 시간에 따른 변화율이 y-50, y축 변화율이 –x+50으로 변했을 때 벡터장의 모습을 관찰하고, 어떤 자연현상과 비슷한지 생각해보자.

 

 

 

-끝-

 

 

 

 

 

  •  
    je Lv.1 2019.06.06 16:13

    이렇게 미분방적식을 풀어보니까 더 쉽게 이해가 되는 것 같아요ㅎ

    댓글 작성하기 좋아요1 댓글수0
  •  
    이용현 Lv.1 2019.06.06 16:18

    미분방정식 설명이 쉽게 풀이되어있어서 보기 좋았고 
    또한 코딩과 접목했기 때문에 직접 지식을 습득하는데 있어서 좋은 방법인 것 같습니다.

    다른 개념들도 잘 알려주셨으면 감사하겠습니다!!

    댓글 작성하기 좋아요0 댓글수0
  •  
    MMCCM Lv.1 2019.06.06 17:36

    정말 쉽게 미분방정식을 이렇게 쉽게 이해할 수 있었습니다.

    댓글 작성하기 좋아요0 댓글수0
  •  
    Undefined Lv.3 2019.06.21 02:17

    \frac{d^2x}{dt^2}이 들어가는 이차 이상의 미분방정식이나

    편미분방정식도 풀 수 있는지 궁금합니다.

    댓글 작성하기 좋아요1 댓글수3
    •  
      코딩수학_이재윤 Lv.1 2019.06.24 11:49

      Undefined 친구!  전세계의 많은 수학자분들이 SAGE를 이용해 미분방정식을 연구하고 있어요.

      아래와 같은 간단한 2차 미분방정식은 단 몇줄의 코딩으로 쉽게 풀 수 있습니다.

       x = var('x')
       y = function('y')(x)
       de = diff(y,x,2) - y == x
       desolve(de, y)
      
      
      아래와 같은 고차 미분방정식도 풀 수 있어요.  그리고 (변수를 정의한 후) 코딩을 단 한줄로 축약할 수도 있습니다. 
      

      x = var('x')
      y = function('y')(x)
      desolve(diff(y,x,2)+y*(diff(y,x,1))^3==0,y,[0,1,1,3]).expand()

      (매우 복잡한 미분방정식의 연산은 변환과정 및 여러가지 근사법 등 수학적인 테크닉이 더 필요해요)

       

      편미분방정식(PDE)의 경우에는 SAGE에서도 해결가능하지만, 연산 알고리즘을 직접 코딩해야 합니다.

      Python을 다룰 줄 안다면, FiPy 패키지를 설치하여, 이를 이용해 PDE를 푸는 것이 가장 쉬운 방법입니다.

      좋아요1
    •  
      코딩수학_이재윤 Lv.1 2019.06.25 11:21

      편미분방정식(PDE) 중에서 가장 간단한 예제인 열방정식(Heat Equation)에 대한, 코딩을 공유드립니다.

      코딩창에 붙여넣어서 직접 해보고,  원하는 편미분방정식의 형태로 수정해서 사용하시면 됩니다.

      var('x t') 
      f = function('f')(x)
      g = function('g')(t)
      z = f*g
      eq(x,t) = diff(z,x,2) == diff(z,t)
      show(eq(x,t))
      eqn = eq/z
      show(eqn(x,t))
      k = var('k')
      eq1(x,t) = eqn(x,t).lhs() == k
      eq2(x,t) = eqn(x,t).rhs() == k
      g(t) = desolve(eq2(x,t),[g,t])
      show(g(t))
      assume(k>0)
      show(desolve(eq1,[f,x]))

      좋아요0
    •  
      코딩수학_이재윤 Lv.1 2019.06.28 10:34

      추가로,  다른 친구들이 3차원 벡터장에 대해서 궁금해 할 거 같아  관련 코딩명령어를 공유합니다.

      var('x y z')
      plot_vector_field3d((2*x*y,-3*z,z),(x,0,2),(y,0,2),(z,0,2))

      좋아요0
  • 폴리매스 문제는 2019년도 정부의 재원으로 한국과학창의재단의 지원을 받아 수행된 성과물입니다.

  • ☎문의 02-6749-3911