첫 번째 수치 문제 컴퓨터시뮬레이션학과 담당교수 : 이형원 E304호, hwlee@inje.ac.kr 운동시뮬레이션 제2주 첫 번째 수치 문제 컴퓨터시뮬레이션학과 담당교수 : 이형원 E304호, hwlee@inje.ac.kr

제1장 첫 번 째 수치문제 방사성 붕괴 수치해석 방법 Modelica 표현 Modelica 클래스 작성 결과 해석

소개 물리학의 많은 문제는 미분방정식을 포함하고 있다. 일차 미분 방정식의 풀이 포사체운동, 단진동 운동, 행성운동 미분 방정식을 푸는 기술 습득 프로그램을 작성하는 과정 습득 Modelica를 활용한 문제 해결방법 습득

방사성 붕괴 많은 핵은 불안정하다는 것이 알려져 있다. 핵붕괴는 무작위(無作爲)과정이다. 동위원소 235 𝑈 (143중성자, 92양성자)는 확률이 작기는 하지만 두 개의 대충 반 크기의 핵과 양성자, 중성자, 전자, 알파입자로 붕괴한다. 핵붕괴는 무작위(無作爲)과정이다. 한 우라늄핵이 정확히 언제 붕괴할지 모른다. 붕괴할 확률을 줄 수 있다. 붕괴하는 데 걸리는 평균 시간을 줄 수 있다. 우라늄의 경우 수명이 1.0× 10 9 년이다.

방사성 붕괴 알파선 : 핵이 헬륨이온을 내면서 다른 핵자로 바뀔 때 나오는 방사선 알파선 : 핵이 헬륨이온을 내면서 다른 핵자로 바뀔 때 나오는 방사선 베타선 : 핵이 전자를 내면서 다른 핵자로 바뀔 때 나오는 방사 선 감마선 : 여기된 핵자가 안정되면서 나오는 방사선으로 전자기 파이다. 측정 장치 : 가이거 카운터(Geiger Counter) 소리1 소리2

방사성 붕괴

방사성 붕괴 수많은 핵자가 모여있는 것을 가정하면 핵자의 수는 다음 미분 방정식을 만족한다. 𝑑 𝑁 𝑈 𝑑𝑡 =− 𝑁 𝑈 𝜏 식(1.1) 직접 대입을 통하여 다음이 해인 것을 알 수 있다. 𝑁 𝑈 = 𝑁 𝑈 (0) 𝑒 −𝑡/𝜏 식(1.2), 𝜏는 평균 수명

수치해석 방법 위 미분 방정식은 해석적으로 해를 구할 수 있음. 수치적으로 해를 구하는 방법을 적용할 것임. 풀이의 목적은 𝑁 𝑈 를 임의시간에서의 값을 구하는 것임. 특별한 시간(보통은 0)에서 주어진 값을 안 후 이후의 시간에서의 값 을 구함. 이렇게 해를 구하는 것을 초기값문제(initial value problem) 이라고 함.

수치해석 방법 수치적인 해를 위하여 테일러(Taylor) 전개를 사용함. 𝑁 𝑈 ∆𝑡 = 𝑁 𝑈 0 + 𝑑 𝑁 𝑈 𝑑𝑡 ∆𝑡+ 1 2 𝑑 2 𝑁 𝑈 𝑑 𝑡 2 ∆𝑡 2 +⋯식(1.3) 위 식은 ∆𝑡에서의 핵자의 수는 0에서의 핵자의 수와 핵자 수 의 미분을 알면 구할 수 있음을 보여준다.

수치해석 방법 ∆𝑡가 매우 작은 경우 2차항 이상을 무시할 수 있다. 𝑁 𝑈 ∆𝑡 ≈ 𝑁 𝑈 0 + 𝑑 𝑁 𝑈 𝑑𝑡 ∆𝑡 식(1.4) 위 결과는 미분의 정의 𝑑 𝑁 𝑈 𝑑𝑡 ≡ lim ∆𝑡→0 𝑁 𝑈 𝑡+∆𝑡 − 𝑁 𝑈 (𝑡) ∆𝑡 ≈ 𝑁 𝑈 𝑡+∆𝑡 − 𝑁 𝑈 (𝑡) ∆𝑡 식(1.5) 식을 정리하면 𝑁 𝑈 𝑡+∆𝑡 ≈ 𝑁 𝑈 𝑡 + 𝑑 𝑁 𝑈 𝑑𝑡 ∆𝑡 식(1.6)

수치해석 방법 기하학적 의미 𝑁 𝑈 (∆𝑡) 𝑁 𝑈 (∆𝑡) 𝑒𝑥𝑎𝑐𝑡 ∆𝑡

수치해석 방법 미분의 정의 이 식에 의하여 (Euler method) 𝑑 𝑁 𝑈 𝑑𝑡 = lim ∆𝑡→0 𝑁 𝑈 𝑡+∆𝑡 − 𝑁 𝑈 (𝑡) ∆𝑡 ≈ 𝑁 𝑈 𝑡+∆𝑡 − 𝑁 𝑈 (𝑡) ∆𝑡 𝑁 𝑈 𝑡+∆𝑡 ≈ 𝑁 𝑈 𝑡 + 𝑑 𝑁 𝑈 𝑑𝑡 𝑡=𝑡 ∆𝑡

수치해석 방법 𝑁 𝑈 𝑡+∆𝑡 ≈ 𝑁 𝑈 𝑡 − 𝑁 𝑈 (𝑡) 𝜏 ∆𝑡 방사성 붕괴 방정식에 대입 하면 순차적으로 값 계산(initial value problem; 초기값 문제) 𝑑 𝑁 𝑈 𝑑𝑡 =− 𝑁 𝑈 𝜏 𝑁 𝑈 𝑡+∆𝑡 ≈ 𝑁 𝑈 𝑡 − 𝑁 𝑈 (𝑡) 𝜏 ∆𝑡 𝑁 𝑈 0 → 𝑁 𝑈 ∆𝑡 → 𝑁 𝑈 2∆𝑡 → 𝑁 𝑈 3∆𝑡 ⋯

수치해석 방법 위 미분방정식을 변형한 Euler 방법을 구현하는 알고리즘을 구 현한 프로그램 작성 수치해석 방법 Euler Method Runge-Kutta Method Differential Algebraic System Solver(DASSL)

수치해석 방법 위 과정으로 구한 해는 정확한 해가 아니고 근사해라는 것을 항 상명심해야 한다. 수치적인 해를 구하는 것에 있어서의 또 다른 목적은 수치적으 로 구한 근사해와 정확한 해 사이의 차이를 무시할 수 있도록 하는 방법을 찾는 것이다. 이 방법은 문제에 따라 달라진다.

결과 그래프

General guideline for test 1 결과가 그럴 듯 한가? 어떤 수치 계산을 하더라도 미리 결과가 대강 어떨 것이다 하는 것은 알고 있어야 한다. 결과가 직관과 직감에 맞는지 확인해야 한다. 이 과정을 통하여 문제에 대한 이해를 높일 수 있다. 남에게 결과를 설명할 때 상대방이 확신을 가질 수 있도록 할 수 있어 야 한다.

General guideline for test 2 결과가 아는 정확한 결과와 일치하는가? 핵붕괴 문제는 해석적인 정확한 해를 알고 있기 때문에 결과를 비교할 수 있다. 일반적으로는 특별한 경우에만 정확한 해를 알 수 있는 경우가 대부분 이다. 해를 아는 특별한 경우에 대하여 반드시 결과를 얻어서 정확한 것과 비교하는 것이 프로그램이 동작한다는 것을 보여주는 필요조건(충분조 건은 아님)이다.

General guideline for test 3 항상 크기가 다른 스텝으로 계산한 결과가 일치한다는 것을 점 검해야 한다. 핵붕괴 프로그램은 시간간격 dt 가 있다. 다른 프로그램도 유사하게 스텝 또는 그리드 크기가 있다. 올바른 프로그램은 스텝크기와 독립적으로 같은 결과를 주어야 한다.

Discretization(이산화) 핵붕괴 문제를 풀기 위하여 우리는 시간은 일정한 간격으로 나 누었다. (discretization) 즉 연속 미분 방정식을 식(1.7)과 같이 이산 차분 방정정식으로 바꾸었다. 이를 통하여 연속의 시간에서 핵자의 수를 구하는 것이 아니고 이산시간 점들(𝑛∆𝑡)에서만 핵자의 수를 구한다. 이러한 시간 또는 공간의 이산화는 일반적인 방법이다.

Discretization(이산화) 이렇게 이산화 때문에 발생하는 오차를 무시할 수 있는지를 어 떻게 알 수 있는가? 어떻게 이산화 한 변수의 스텝크기를 정할 수 있는가? 위 질문에 대한 정답은 없으며 단지 답을 구한 방향(guideline) 을 제시할 수 있다.

Guideline for discretization 1 같은 계산을 여러 스텝크기로 해보아야 한다. 다음 페이지의 그림은 핵붕괴 프로그램을 평균수명은 1초, 스텝크기를 0.5, 0.2, 0.05초로 해서 각각 계산한 결과이다. 스텝크기가 작아짐에 따라 결과가 정확한 값에 빠르게 수렴함을 알 수 있다. 식 (1.6)에서 무시한 오차가 ∆𝑡 2 이므로 예산된 결과이다. 즉 ∆𝑡를 줄 이면 결과가 좋아진다.

Guideline for discretization 1 시간이 𝑡만 큼 지난 점에서의 오차는 ∆𝑡스텝크기로 𝑡/∆𝑡 스텝만큼 갔 으므로 오차의 크기는 ~ 𝑡/∆𝑡 × ∆𝑡 2 ~∆𝑡 정도 이다. 이 오차를 전역 오차(global error)라고 한다. 따라서 오일러(Euler) 방법으로 계산한 결과와 정확한 값 사이의 오차 는 ∆𝑡를 반으로 줄이면 오차가 반으로 줄어들게 될 것이다.

Guideline for discretization 1 ∆𝑡=0.5𝑠 ∆𝑡=0.2𝑠 ∆𝑡=0.05𝑠 Solidline exact result

Guideline for discretization 1 핵붕괴 문제는 정확한 해가 있어서 정확도를 계산하는 것이 쉽 다. 정확한 해를 모를 경우에는 어떻게 하는가? 스텝크기를 줄였을 때 모든 점에서 한 값으로 수렴하는지를 판단해야 한다. 스텝크기는 문제의 특성시간보다 작게 잡아야 한다. 핵붕괴인 경우에 는 평균수명 특성 시간의 수 퍼센트가 적당하다.

Guideline for discretization 1 좋은 수치 알고리즘도 어떤 특정한 경우에는 매우 불안정한 경 우가 있음. 오일러 방법은 핵붕괴에 대해서는 잘동작하지만 3장에서 배우 는 진동문제에서는 답을 주지 못한다. 모든 문제에 다 잘 동작하는 수치 알고리즘은 없으며, 풀고자 하는 문제에 대한 이해와 여러분의 천재성을 발휘하여 알고리 즘을 잘 선택해야 한다.

Modelica 소개 Modelica는 주어진 미분 방정식에 대한 수치해를 구하는 프로 그램을 자동으로 생성해 준다. 사용자는 원하는 알고리즘만 선택하면 된다. dassl, euler, rungekutta, dasslwort, radau1,3,5, lobatto2,4 객체지향 방법을 기반으로 작성한다.

Modelica의 Solver DASSL : Differential Algebraic System Solver 𝐹 𝑡, 𝑦, 𝑦 ′ =0 의 해를 일반적으로 구해줌 EULER : Euler 방법으로 해를 구한다. 𝑦 ′ =𝑓(𝑡,𝑦) 의 해를 Euler 방법으로 구한다. Runge-Kutta : Runge-Kutta 방법으로 해를 구한다. 𝑦 ′ =𝑓(𝑡,𝑦) 의 해를 구한다.

Modelica의 Solver DASSLWORT : 내부에서 근 찾는 과정이 없는 DASSL과 동일 Radau1,3,5 : implicit RK 방법을 이용하여 변화가 심한 미분 방 정식의 해를 구한다. LOBATTO2,4 : RK를 개선한 Lobatto 방법에 의한 해를 구한다.

Modelica의 구성 요소 package : 용도별 클래스 모음 class : Modelica의 기본 구성 요소 model : class 와 동일함 block : model과 동일하지만 input/output 연결만 허용 function : 알고리즘을 구현하는데 사용 record : 자료저장을 위한 클래스 type : 자료형을 위한 클래스 operator : function 의 모음 connector : 방정식을 허용하지 않고 연결을 정의하는 클래스

Modelica의 구성 요소 Basic types Real : 실수 값 Integer : 정수 값 Boolean : 참, 거짓(true, false) String : 일반적인 문자열 Complex : 복소수

Modelica의 표현 미분의 표현 미분 방정식의 표현 방정식의 표현 Modelica의 식은 대입문이 아니고 수학의 방정식임. der(y) : y의 시간에 대한 미분 미분 방정식의 표현 der(y) = -y; 𝑑𝑦 𝑑𝑡 =−𝑦 의 표현식 방정식의 표현 A*x + b*y = c; 𝑎𝑥+𝑏𝑦=𝑐 의 표현식 Modelica의 식은 대입문이 아니고 수학의 방정식임.

Modelica 클래스 작성 𝑑𝑁 𝑑𝑡 =− 𝑁 𝜏 에 대한 클래스 작성

시뮬레이션 결과

결과 해석 핵자의 수가 이론적인 결과에서 알 수 있는 바와 같이 지수적으 로 감소함을 알 수 있다. 커서를 결과 그래프 위에 놓으면 다음 그림과 같은 수치를 볼 수 있다. 100× 𝑒 − 3.0061 3 =36.713 의 값과 거의 일치함을 알 수 있다.

Modelica 클래스 작성 과정 Motion.y2016.Week02 package 생성 Package 안에 NucDecay 클래스 작성 NucDecay 클래스 안에 방사성 붕괴를 기술하는 미분방정식 기 술 시뮬레이션 조건을 설정하고 시뮬레이션 실행

실습의 목표 Modelica로 방사성 붕괴 미분 방정식 풀기 결과 그래프 꾸미기 여러 결과 같이 그리기

불필요한 Package/Class 닫기 Motion 패키지가 사라졌음

Package 생성 패키지 명 Motion 패키지 선택 새로운 패키지 Motion이 root(global) 에 생성됨

Package 생성 패키지 명 y2016 패키지 선택 상위 패키지 새로운 패키지 y2016가 Motion 아래에 생성됨

Package 생성

Package 저장 D:\lec_hwlee\motion\y2015\week02 Motion.y2015.Week02.mo

클래스 생성 새로운 클래스 NucDecay가 Motion.y2016.Week02 아래에 생성됨

클래스 작성 시간에 따라 변하는 변수는 핵자의 수를 나타내는 N이다. 계의 조건에 따른 변수는 평균 수명 타우(𝜏)와 초기 핵자의 수 ( 𝑁 0 ) 이다.

클래스 작성 앞 클래스 기술에서 N은 변하는 임의의 변수로 구하고자 하는 시간의 함수이다. 타우(𝜏)와 𝑁 0 는 시뮬레이션 하는 동안 변하지 않는 값으로 파라 메터(parameter)로 선언한다. 6번 줄은 N이 만족하는 미분 방정식이다.

시뮬레이션 조건 설정

시뮬레이션 조건 설정 시뮬레이션 시작시간 종료시간 설정 미분방정식 Solver 선택 오차 한계

시뮬레이션 조건 설정 결과 출력 형식 mat : 이진 형식 출력, MATLAB, Octave 에서 사용 가능 plt : 일반 텍스트 출력 csv : 자료를 콤마로 구분하여 저장 empty : 출력하지 않음 적분 구간의 개수 ∆𝑡= 종료시간 −시작시간 인터벌 수 ∆𝑡= 30−0 500 =0.06초

시뮬레이션 실행 콤파일하고 시뮬레이션한 로그를 보여준다.

시뮬레이션 실행

결과 보이기 보고자 하는 변수 선택

결과 그래프 꾸미기 Setup을 통하여 그래프의 형식을 조정할 수 있다.

결과 그래프 꾸미기 그래프 타이틀 y축 타이틀 x축 타이틀 범례 선색 선굵기

결과 그래프 꾸미기

여러 결과 같이 그리기 여러 결과를 같이 그리기 위해서는 새로운 클래스를 생성하여 시뮬레이션을 해야 한다. 여러 결과를 같이 그리기 위해서는 새로운 클래스를 생성하여 시뮬레이션을 해야 한다. NucDecay1 클래스를 만들어 풀이 방법과 시간간격을 달리하여 그려 볼 것이다. NucDecay1 클래스는 NucDecay 클래스와 동일하다.

다른 시뮬레이션 조건 설정

다른 시뮬레이션 조건 설정 미분방정식 Solver를 Euler로 선택

다른 시뮬레이션 조건 설정 인터벌 수를 60으로 하여 적분 시간간격을 0.5로 설정

시뮬레이션 실행

시뮬레이션 결과

결과 표시(# of Intervals = 500)

Number of Intervals = 10 의 결과

Number of Intervals = 20 의 결과

결과 해석 적분 시간 간격을 어떻게 설정 하느냐에 따라 결과가 상이하게 나타나는 것을 알 수 있다. 적분 시간 간격을 어떻게 설정 하느냐에 따라 결과가 상이하게 나타나는 것을 알 수 있다. 따라서 요구하는 정확도에 따라 적분 방법과 적분 시간 간격을 잘 선택해야 한다.

도전해 보기 방사성 붕괴 프로그램에서 Euler 방법으로 0초부터 20초까지 시 뮬레이션 하면서 시간 간격을 0.5, 0.2, 0.05초로 계산한 결과와 0.5초 간격으로 rungekutta 방법으로 계산한 결과를 한 그래프 에 그리시오. 연습문제 1.1~1.6 해보기