極限 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月)]