極限 limn→∞ (1+t/n)n を用いて、 指数関数 exp(t) を計算する。 とりあえず、 n と t を入力して f(t) = (1+t/n)n を計算するプログラムを 作ってみよう。
プログラムの中枢部は、function
文を用いて書く。
C ----------------------------------------------------- real function f(n,t) implicit none integer k,n; real t,r; if( n .le. 0 ) then write(0,*) 'エラー: n = ', n ; stop endif f = 1.0 ; r = 1.0 + t/n ; do k=1,n f = f * r enddo return end C -----------------------------------------------------
これで なぜ (1+t/n)n が求められるのか、しくみを考えること。
メインルーチンは、たとえば次のようになる。(なお、dt
は
結果の出力の間隔を決めるのに用いているだけで、計算自体には関係ない。)
C 指数関数を、次の極限を用いて計算: C exp(t) = lim (1+t/n)^n C ----------------------------------------------------- implicit none real t, tmax, dt, f; integer n, j, jmax; call get1i('大きな整数 n = ', n ) call get2r('tmax, dt = ', tmax, dt ) write(*,*) '# ', '(n, tmax, dt) = ', n, tmax, dt jmax = int(tmax/dt+0.5); do j=0,jmax t = j*dt; write(*,*) t, f(n,t) enddo return end C -----------------------------------------------------
極限による指数関数の定義を利用して、exp(it) を計算してみる。
計算の主要部分は、次のようになる:
complex function f(n,t) implicit none integer k,n; real t; complex z; if( n .le. 0 ) then エラー発生; stop; endif f = 1.0 z = complex(1.0, t/n) do k=1,n f = f * z enddo
途中の complex(実部,虚部)
という書き方に注目しよう。
メインルーチンで結果を出力するときには、次のように、
RealPart()
および aImag()
を用いて
実部と虚部を分ける:
do j=0,jmax t = j*dt; w = f(n,t); write(*,*) t, RealPart(w), aImag(w) enddo
あとは自分で書いてみよう。
くれぐれも implicit none
を忘れないように。
[数値計算] [力学演習(2002年11月)]