第1回レポート●以下のプログラムを作成しメールで提出「級数展開によりsin(x π)を求め、表示する」 ● sin(x π)の級数展開:以下で近似する.
●   xとnは実行時にキーボードから入力. ●実数の計算は全て倍精度で行う. ● 提出先: z5nt02in＠s.kyushu-u.ac.jp ●締め切り: 6月1日(水)必着
●実行例: ● nが1以上なら計算結果を表示% ./rep1 x : 0.3333 n: 10 sin(x * pi) is about 0.8660202090600015 ● nが 0以下の場合,エラーを表示して終了% ./rep1 x : 0.25 n: 0 Error!  n must be > 0エラーメッセージ
●字下げ,改行を使い,見やすいプログラムを書く. ●必ず自分で実行し,正しい結果が得られることを確認●  xやnに様々な値を入力して確認する. ● sin関数を使った結果と比較● わざとnに負の値などを入力してエラーの表示も確認する. ● ソースファイルをメールに添付して提出。プログラム名を間違えないよう注意●末尾に ~が付いたファイルは最新版ではない●末尾が. f90でないファイルはソースファイルではない●本文には苦労した点などを書く.
●必ず返信を確認する. ● 翌日の午前中までに返信予定間違い個所の指摘,ヒントなど● ●締め切り前であれば何度でも再提出可
● Webページも確認.
●提出状況. ● 翌日の午前中までに反映予定● レポートに関する追加情報.
●日時 5月31日(火) 17:00〜19:00 ●場所第1講義室●レポートや講義内容に関する質問受け付け. ●それ以外の時間の質問は,とりあえずメールで.
program report1 implicit none integer :: n, i, j real(8) :: x, total, pi, rad, f, g intrinsic dble, atan pi = 4.0D0 * atan(1.0D0) write(*, *) 'Enter x (times pi): ' read(*, *) x rad = pi * x何度も利用するpi * xを予め計算しておくwrite(*, *) 'Enter n: ' read(*, *) n if (n <= 0) then write(*, *) 'Error: n<=0' stop endif
total = 0.0D0 f = 1.0D0 do i = 1, n g = 1.0D0 ! calculate (2i-1)! do j = 1, 2 * i - 1 g = g * dble(j) end do total = total + f * rad**(2*i-1) / g f = -1.0D0 * f ! flip flag end do (-1.0D0)**(i+1)だとi+1回の掛け算が必要write(*, '(a, f7.3, a, f15.10)') & 'sin (', x,'*pi) is about ', total stop end program細かい注意点●整数と実数●字下げ● ● ● ● ● ● ● ���� ) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● o do i=1,n a = 1.0D0 do j=1,2*i-1 a = a * (x*pi/j) end do total = total + (-1)**(i+1) * a end do do i=1,n a = 1.0D0 do j=1,2*i-1 a = a * (x*pi/dble(j)) end do total = total + (-1.0D0)**(i+1) * a end do -実数式の中で整数を使う場合は実数に変換する. -ただし,べき乗する数は整数でも可
-字下げの幅は一定にしたほうが見やすい. -文の意味が変化する箇所に改行を入れると見やすい.
●実数の表示形式●桁あふれ対策1 :整数より倍精度実数●桁あふれ対策2:計算の順序
●書式を指定せずに実数を出力した場合●プログラムwrite(*, *)  'sin (', x,'*pi) is about ', total ●出力結果sin ( 1.00000000000 *pi) is about  3.458669144327514E-16 ●  3.458669144327514E-16   =3.458669144327514   x 10-16
● 実数の表示形式をＦにする
●プログラム全体で15桁,小数点以下10桁write(*, '(a, f7.3, a, f15.10)') & 'sin (', x,'*pi) is about ', total o出力結果sin (  1.000*pi) is about 0.000000000小数点以下10桁では表示できないくらい小さい
● nやxが大きい場合(2i-1)!が非常に大きくなり,計算機で普通に扱える範囲を超える. →桁あふれ。データの型(整数,単精度実数,倍精度実数)によって扱える範囲が違う
●整数より実数●計算機で扱える値の範囲整数-231  ～ 231 - 1 (= -2,147,483,648 ～ 2,147,483,647) ●倍精度実数    約±1.8 x 10³08まで●ちなみに絶対値の小さい値は約±4.9 x 10 -3²²まで▪それより小さくなると0となる. (underflow)
●計算の順番を変える● 足し算と引き算, 掛け算と割り算を交互に行う. ●レポートの例
program report1 implicit none integer :: n, i, j real(8) :: x, total, pi, rad, f, g intrinsic sin, atan pi = 4.0D0 * atan(1.0D0) write(*, *) 'Enter x (times pi): ' read(*, *) x rad = x * pi write(*, *) 'Enter n: ' read(*, *) n if (n <= 0) then write(*, *) 'Error: n<=0' stop endif total = 0.0D0 £ = 1.0D0 do i = 1, n g = 1.0D0 do j = 1, 2 * i - 1 g = g * rad / dble(j) end do total = total + f * g f = -1.0D0 * f       ! flip flag end do write(*, '(a, £7.3, a, £15.10)') & 'sin (', x,'*pi) is about ', total stop end program
●計算の節約● nが整数か否か判別する●正しい値が入力されるまで繰り返す●ガンマ関数を使ったべき乗計算
●同じ計算を何度も行っている
計算回数と所要時間。このプログラムでは掛け算と割り算をそれぞれ n²回( ∑ n 2i-1回)サーバah.s.kyushu-u.ac.jpの計算スピードは, ●どんなに頑張っても約2.2x10 9回/秒●実際は割り算に時間がかかるので約3.0x107回/秒程度. ●所要時間:約3秒n = 1000 の時: n = 10000の時:  約300秒●
●直前の計算結果を利用total, = x | Y 1 = X X X | Yi = (−1) • Yi-1 • (2i−2) (2i−1) ( i = 2 … n)
program report1 implicit none integer :: n, i, h real(8) :: x, total, pi, rad, y intrinsic sin, atan pi = 4.0D0 * atan(1.0D0) write(*, *) 'Enter x (times pi): ' read(*, *) x rad = x * pi write(*, *) 'Enter n: ' read(*, *) n if (n <= 0) then write(*, *) 'Error: n<=0' stop endif total = rad y = rad do i = 2, n ! calculate Y_i = -1 * Y_(i-1) * x * x / ((2i-1)*(2i-2)) h = (2 * i - 2) * (2 * i -1) y = -1.0D0 * y * rad * rad / dble(h) ! calculate TOTAL_i = TOTAL_(i-1) + Y_i total = total + y end do write(*, '(a, f7.3, a, f15.5)') & 'sin (', x,'*pi) is about ', total stop end program
real(8) :: n • • • write(*,*) '級数のnは？' read(*,*) n c = ceiling(n) if (n <= 0 .or. n /=c) then write(*,*) "nは自然数にしてね♪" stop end if
write(*,*) "n=?" read(*,*) n do while (n <= 0.0D0) write(*,*)"Error! n must be > 0. Please try again!" write(*,*) "n=?" read(*,*) n end do補足:   do while(条件)   〜  end do条件が満たされている間, end doまでの処理を繰り返す.
● -1.0～1.0になるまで1.0を加減y = x do while (y >= 1.0D0) y = y -1.0D0 end do do while (y <= -1.0D0) y = y +1.0D0 end do ●他の方法: 。切り捨てた値を利用y = x - int(x)
degtorad = pi / 180.0D0 write(*, *) '入力値xの種別' write(*, *) '1:x radian/2: πx radian/other:degree ? ' read(*, *) type select case(type) case('1') xrad = x xpi = x / pi xdeg = x / degtorad case('2') xrad = x * pi xpi = x xdeg = x * 180.0D0 case default xrad = x * degtorad xpi = x / 180.0D0 xdeg = x end select計算方法の選択も同様
total = 0.0D0 f = 1.0D0 do i = 1, n, 1 total = total + f * rad**(2*i-1)/gamma(2.0D0*dble(i)) f = -1.0D0 * f end do ●ガンマ関数: l e-t t x-l dt ( Γ( x ) = ●ガンマ関数の性質|Г(х+1) = хГ(х)より|Г( n +1) = n!