演習1

Size: px

Start display at page:

Download "演習1"

ゆりかしばもと
4 years ago
Views:

1 神戸市立工業高等専門学校電気工学科 / 電子工学科専門科目数値解析演習 1 山浦剛 (tyamaura@riken.jp) 講義資料ページ

2 Fortran とは? Fortran(= FORmula TRANslation ) は 1950 年代に誕生した世界初の高級プログラミング言語数値計算プログラム作成に適している固定形式 : FORTRAN 66, 77 FORTRAN 自由形式 : Fortran 90/95, 2003, 2008 Fortran Fortran がなぜ数値計算プログラム作成に適しているか? 設計思想として数値計算用プログラミング言語に特化している組込関数や複素数配列操作など数値計算に便利な機能が予め組み込まれている他の高級言語 (C/C++/Java など ) と比較して数値計算を簡潔に記述することができ一般的にパフォーマンスに優れている ( 言語仕様上コンパイラにとって最適化しやすい ) LAPACK 等優れた数値計算ライブラリが豊富に利用可能移植性に優れ PC からスーパーコンピュータまで様々な環境で利用可能プログラミングにおいての誤りを犯しにくい

3 Fortran プログラムの基本構造大文字小文字の区別はないコメントは! 以降文の終了を宣言しない (; などをつけない ) 文を継続するには行末に & をつけるプログラムファイルの拡張子は f90 例 : sample01-1_hello.f90 Fortran コンパイラ (gfortran) でコンパイル例 : gfortran sample01-1_hello.f90 デフォルトの出力ファイル名は a.out 例 :./a.out -o オプションで出力ファイル名の変更も可基本形プログラム名型宣言処理を記述する program [ プログラム名 ] implicit none [ 型宣言 ] [ 処理 ] end program

4 Fortran の型宣言と演算よく使う変数の型単精度整数型 : integer 単精度実数型 : real または real(4) 倍精度実数型 : double precision または real(8) 単精度複素数型 : complex 定数 : ( 数の型 ), parameter 他にも文字列型や論理型などがある宣言の例 integer :: i, j, k real(8) :: x, y, z real(8), parameter :: pi = 演算代入 : xx=yy 加算 : xx+yy 減算 : xx yy 乗算 : xx yy 除算 : xx/yy べき乗 : xx yy 組み込み関数三角関数 : sin(xx), cos(xx), atan(xx), sinh(xx) 対数関数 : log(xx), ln(xx) 指数関数 : exp(xx) 絶対値関数 : abs(xx)

5 Fortran の論理構造条件分岐論理値を使って分岐を発生させる例 : if/else 文による絶対値の求め方 a = -5 if( a < 0 ) then b = -a else b = a 反復計算例 : x に 1~100 までの和を代入する do i = 1, 100 x = x + i 出力 (write 文 ) 引数 1: 出力先 ( 標準出力 =*) 引数 2: フォーマット ( 自動選択 =*) write(*,*) hello write(*,*) a=,a

6 Fortran の型変換型変換 Fortranでは castは自動で行われる異なる型の演算対象が含まれる演算では強い方の型 ( 種別も含む ) が使用される整数 : 弱実数 : 中複素数 : 強整数 + 実数実数複素数 + 実数複素数演算対象が実数もしくは複素数のみで構成される場合高い方の精度が採用される単精度実数 + 倍精度実数倍精度実数倍精度実数 + 単精度複素数倍精度複素数整数の除算では切り捨てが発生する 10 / / 2 1

7 復習 : 2 分法 2 分法のアルゴリズム 1. ff aa < 0, ff bb > 0 となる初期値 aa, bb を決める 2. cc aa+bb 2, dd aa bb 2 d < εε ならば 3 に移るを計算 d εε, ff cc < 0 ならば aa に cc を代入し (aa cc) 2 を繰り返す d εε, ff cc > 0 ならば bb に cc を代入し (bb cc) 2 を繰り返す 3. cc を数値解とする program bisection_method implicit none integer :: n real(8) :: a, b, c, fc a = 1.0 b = 2.0 do n = 1, 40 c = ( a + b ) / 2.0 fc = c**3-5.0 if( fc < 0.0 ) then a = c else b = c write(*,*) 'step number =',n,'; c=',c end program

8 復習 : ニュートン法ニュートン法のアルゴリズム 1. 初期値 xx 0 と許容する誤差 εε を決める 2. cc xx ff(xx) ff (xx) を計算 cc xx cc < εε ならば 3 に移りそうでなければ xx に cc を代入し 2 を繰り返す 3. cc を数値解とする program newton_method implicit none integer :: n real(8) :: x, fx, dfx x = 2.0 do n = 1, 10 fx = x**3-5.0 dfx = 3.0 * x**2 x = x - fx / dfx write(*,*) 'step number =',n,'; x=',x end program

9 課題 1. 2 分法のサンプルプログラムを参考に誤差 εε = 10 6 とし誤差未満まで計算できたときに終了するように修正せよ ( ヒント : do ループは exit 文で抜けることができる ) 2. ニュートン法のサンプルプログラムを参考に誤差 εε = 10 6 とし誤差未満まで計算できたときに終了するように修正せよ 3. 次の式の根についてニュートン法を用いて計算するプログラムを作成せよただし誤差 εε = 10 6 として誤差未満まで計算できたときに終了させること x の初期値を 2 とする ff xx = exp xx 2 sin xx 4. 次の式の根について 2 分法とニュートン法を用いて計算するプログラムを作成せよただし 2 分法での誤差はεε = 10 2, ニュートン法での誤差はεε = 10 6 として誤差未満まで計算できた dd(arctan xx ) ときに終了させること xの初期値を-4, 4とする ( ヒント : = 1 dddd 1+xx 2) ff xx = arctan(xx) 1 2

10 提出方法〆切 : 2019/05/16( 木 ) メールにプログラムを添付主題 : 演習 1 レポート ( 学籍番号 ) 宛先 : tyamaura@riken.jp 添付 : 学籍番号 _ 課題番号.f90 を 4 ファイル課題 1-1: r???????_kadai01-1.f90 課題 1-2: r???????_kadai01-2.f90 課題 1-3: r???????_kadai01-3.f90 課題 1-4: r???????_kadai01-4.f90

11 課題 1: 解答例ループの終了条件を意識誤差 ε 1.0E-6 (= ) program bisection_method implicit none! definition integer :: n real(8) :: a, b, c, d, fc real(8) :: eps! ! initial value a = 1.0 b = 2.0 eps = 1.0e-6! bisection method loop do n = 1, 40 c = ( a + b ) / 2.0 d = abs( a - b ) / 2.0 fc = c**3-5.0 if( fc < 0.0 ) then a = c else b = c! output write(*,*) 'step number =',n,'; c=',c if( d < eps ) then exit end program step number = 1 ; c= step number = 2 ; c= step number = 3 ; c= step number = 4 ; c= step number = 5 ; c= step number = 6 ; c= step number = 7 ; c= step number = 8 ; c= step number = 9 ; c= step number = 10 ; c= step number = 11 ; c= step number = 12 ; c= step number = 13 ; c= step number = 14 ; c= step number = 15 ; c= step number = 16 ; c= step number = 17 ; c= step number = 18 ; c= step number = 19 ; c= step number = 20 ; c=

12 program newton_method implicit none! definition integer :: n real(8) :: x, nx, fx, dfx real(8) :: eps! step number = 1 ; x= step number = 2 ; x= step number = 3 ; x= step number = 4 ; x= 課題 2: 解答例ループの終了条件を意識誤差 ε 1.0E-6 (= ) 鍵となる変数 (x) の扱いに注意! initial value x = 2.0 eps = 1.0e-6! newton method loop do n = 1, 10 fx = x**3-5.0 dfx = 3.0 * x**2 nx = x - fx / dfx! output write(*,*) 'step number =',n,'; x=',nx if( abs( ( nx - x ) / nx ) < eps ) then exit else x = nx end program

13 program newton_method implicit none! definition integer :: n real(8) :: x, nx, fx, dfx real(8) :: eps! step number = 1 ; x= step number = 2 ; x= step number = 3 ; x= step number = 4 ; x= step number = 5 ; x= step number = 6 ; x= step number = 7 ; x= 課題 3: 解答例式が異なるだけで課題 2 同様! initial value x = 2.0 eps = 1.0e-6! newton method loop do n = 1, 10 fx = exp(-x**2) sin(x) dfx = * x * exp(-x**2) - cos(x) nx = x - fx / dfx! output write(*,*) 'step number =',n,'; x=',nx if( abs( ( nx - x ) / nx ) < eps ) then exit else x = nx end program

14 課題 4: 解答例式 yy = atan(xx) はニュートン法で単純に解くと初期値によっては発散してしまう関数二分法で力業で解くこともできるが始めの数回を二分法その値を初期値としてニュートン法を行う実践的な解き方を学ぶのが目的 program bisection_newton_composite implicit none integer :: n real(8) :: x1, x2, d, x real(8) :: nx, fx, dfx real(8) :: eps1, eps2 x1 = -4.0 x2 = 4.0 eps1 = 1.0e-2 eps2 = 1.0e-6! bisection method loop do n = 1, 100 x = ( x1 + x2 ) / 2.0 d = abs( x1 - x2 ) / 2.0 fx = atan(x) if( fx < 0.0 ) then x1 = x else x2 = x write(*,*) 'bisection: step number =',n,'; x=',x if( d < eps1 ) then exit! newton method loop do n = 1, 100 fx = atan(x) dfx = 1.0 / ( x**2 ) nx = x - fx / dfx write(*,*) 'newton : step number =',n,'; x=',nx if( abs( ( nx - x ) / nx ) < eps2 ) then exit else x = nx end program bisection: step number = 1 ; x= bisection: step number = 2 ; x= bisection: step number = 3 ; x= bisection: step number = 4 ; x= bisection: step number = 5 ; x= bisection: step number = 6 ; x= bisection: step number = 7 ; x= bisection: step number = 8 ; x= bisection: step number = 9 ; x= bisection: step number = 10 ; x= newton : step number = 1 ; x= newton : step number = 2 ; x= newton : step number = 3 ; x=

演習2

演習2 神戸市立工業高等専門学校電気工学科 / 電子工学科専門科目数値解析 2017.6.2 演習 2 山浦剛 (tyamaura@riken.jp) 講義資料ページ h t t p://clim ate.aic s. riken. jp/m embers/yamaura/num erical_analysis. html 曲線の推定 N 次多項式ラグランジュ補間 y = p N x = σ N x x