数値計算演習1

数値計算の練習
Fortran言語によるプログラミング
倉橋貴彦
Fortranによるプログラミング
＜プログラム・入力データの作成 (Windows上における作業)＞
・Fortran言語によるプログラムの作成 → Windowsのエディタ(TeraPad等）を使用
（ファイル名は「○○○.f」とする（○は任意の半角英数字））
・入力データを作成する。→ input.dat
・nutechへのプログラム（ ○○○.f ），計算条件（ input.dat ）の転送
→ ffftp（ファイル転送ソフト）を使用
＜nutechを使用＞
・nutechにログインをし，プログラムのコンパイルを行う．
（コンパイルは「f95 ○○○.f」と入力し，enterキーを押す．
エラーがある場合は，コンパイル時にエラーが表示される．）
・プログラムによる計算は，「./a.out」と入力することにより実行される．
・答えのファイルは今回の問題では「output.dat」というファイルに出力される．
例：
以下のプログラムをWindows上のエディタ(メモ帳, Word Pad, Tera Pad)を
使用して作成する．
write(*,*) “hello.”
stop
end
test.fというファイル名で保存をし，ffftpを使用して，nutechにファイル転送する．
ファイルをドラッグ＆ドロップし，nutechに
ファイル(test.f)を転送する．
ホストにnutechと入力
ユーザー名とパスワードを入力
（※初回時は，windowが現れるが，OKを押して，
nutechに入る．）
①f95 test.fと入力しenterを押すことで
コンパイルが行われる．
②./a.outと入力することで，プログラムが実行される．
③ここではhello!と出力される．
同様に各課題のプログラムを作成し，
計算を実行する．
Fortranの基本的ルール
1234567
・プログラムは7カラム目から打つ。
implicit double precision ( a-h, o-z )
・6カラム目はプログラムが長くなり、改行の際に&
を入れる。（1行は72カラムまでに入力する。）
open(10,file='input.dat')
open(11,file='output.dat')
・先頭にcを付けるとコメント文となり、プログラムか
らは除去される。
c
c
read(20,*) x0
read(20,*) y0
c
eps = 0.000001d0
iter = 0
c
222 continue
・・・・・・
go to 222
c
write(11,110) aaa
110 format(f10.5)
・1～5カラムにはread文，write文のフォーマットの
数字や，continue文の文番号を入れる。
・implicit double precision(a-h,o-z)は倍精度による
プログラムを示しており，この定義をしないと単精
度の計算となる。
・a-h,o-zが頭文字になる変数は実数(ex a=1.1d0)，
I,j,k,l,m,nが頭文字になる場合は整数(ex i=120)を
表す．
・整数を入力する場合はd0を後ろに付ける．
・例えばルートの計算にはsqrtという関数を使用す
るが，倍精度の計算の場合dsqrtとする。（組み込
み関数(sqrt以外にもsin,cos等）については，使用
する場合に個々に調べること。）
read文 do文 write文
read文
read(*,*) aa
do文（c言語におけるfor文）
do I = 1,nx
*****
end do
算術演算
和＋
差－
積 *
商 /
累乗 **
例1 : 1～nまでの足し算
1234567
implicit double precision ( a-h, o-z )
c
write文
read(*,*) n
c
ia=0
write(*,*) aa
c
do i = 1,n
ia = ia + i
end do
write文の書式例
write(*,100) i, aa, bb
100 format(i5,1x,2f10.5)
c
＜意味＞整数iを5カラム内に出
力する．1カラム空ける．aa,bbの2
つの変数を10カラム内に小数点
以下5桁で出力する。
write(*,100) ia
100 format(i5)
stop
end
組み込み関数
exp(x) ⇒ ex
sin(x)
⇒ sin(x) （三角関数のxの値はラジアンとする。）
cos(x)
⇒ cos(x)
tan(x)
⇒ tan(x)
sqrt ⇒ 単精度の場合の平方根の計算
dsqrt ⇒ 倍精度の場合の平方根の計算
例2 a=sin(θ)の計算， b=√aの計算
implicit double precision ( a-h, o-z )
c
read(*,*) theta
c
rad = theta * 180.d0 / 3.1415926535d0
c
aa = sin(rad)
bb = dsqrt(aa)
c
write(*,110) aa, bb
110 format(2f15.10)
stop
end
If分の中身の例
aa .lt. bb aa<bb
aa .gt. bb aa>bb
aa. eq. bb aa=bb
aa .le. bb aa<=bb
aa .ge. bb aa>=bb
aa .gt. bb .and. cc .lt. dd
aa .gt. bb .or. cc .lt. dd
if文
If ( aa .lt. bb ) then
aa = cc
else
bb = cc
end if
＜上の文の意味＞
aaがbbより小さければaaにccの値を入れる。そ
うでなければ、bbにccの値を入れる。
例3 aa<bbならaaを出力し、aa<bbでなければ、bbを出力する。
implicit double precision ( a-h, o-z )
c
read(*,*) aa
bb=10.d0
if ( aa .lt. bb ) then
write(*,*) aa
else
write(*,*) bb
end if
c
stop
end
aa>bb and cc<dd
aa>bb or cc<dd
open文：ファイルからのデータの入力や，ファイルへのデータの出力を行う際にopen文
を使用する。
配列：ベクトルや行列をプログラム中において使用する際，配列を使用する。
parameter(md1=100)
implicit double precision (a-h,o-z)
dimension aa(md1,md1), bb(md1), cc(md1)
例4
c
open(10,file='input.dat')
open(11,file='output.dat')
c
read(10,'(a80)')
read(10,*) nx
read(10,'(a80)')
read(10,*) ((aa(i,j), i=1,nx),j=1,nx)
read(10,'(a80)')
read(10,*) (bb(i),i=1,nx)
c
do i = 1,nx
ww = 0.d0
do j = 1,nx
ww = ww + aa(i,j)*bb(j)
end do
cc(i) = ww
end do
c
write(11,*) 'Computational result'
do i = 1,nx
write(11,222) cc(i)
end do
222 format(f10.5)
c
stop
end
行列とベクトルの積の計算
 a11
a
 21
a12  b1   c1 
  

a22  b2  c2 
例 1 2 1  5
2 3 2  8

   
Input.dat
##### number of variable
2
##### matrix aa
12
23
##### vector bb
1
2
計算後
output.dat
Computational result
5.00000
8.00000
mod文：割り算の余りを使う際に使用する。
副プログラム：プログラムの一部を関数化したプログラムであり，call文とsubroutine文
から成る。
例5 : 反復的に行う計算
 a11
a
 21
a12  b1 
 
a22  b2 
(l )
b1 
 
b2 
( l 1)
例
(1)
5
 
8
( 2)
21
 
34
1 2 1 
 2 3  2 

 
1 2 5
2 3 8

 
( 2)
( 3)
左のようなステップ計算をして，
いくつかのステップか置きに
計算結果（右辺ベクトルの値）を出力する。
例5
parameter(md1=100)
implicit double precision (a-h,o-z)
dimension aa(md1,md1), bb(md1), cc(md1)
c
open(10,file='input2.dat')
open(11,file='output2.dat')
c
read(10,'(a80)')
read(10,*) imax
c
read(10,'(a80)')
read(10,*) iout
c
read(10,'(a80)')
read(10,*) nx
c
read(10,'(a80)')
read(10,*) ((aa(i,j), i=1,nx),j=1,nx)
c
read(10,'(a80)')
read(10,*) (bb(i),i=1,nx)
c
c
c
c
c===== ourput of computational result
c
if ( mod(istep,iout) .eq. 0 ) then
c
write(11,*) 'istep=',istep
do i = 1,nx
write(11,222) cc(i)
end do
222 format(f10.5)
c
else
end if
c
c===== step change
c
do i = 1,nx
bb(i) = cc(i)
end do
c
end do
c
stop
end
c
c================================
c
subroutine calcu
& ( nx, aa, bb, cc, md1 )
c
c================================
c
implicit double precision (a-h,o-z)
dimension aa(md1,md1), bb(*), cc(*)
c
do i = 1,nx
ww = 0.d0
do j = 1,nx
ww = ww + aa(i,j)*bb(j)
end do
cc(i) = ww
end do
c
return
end
do istep = 1,imax
c
c===== calculation of matrix vector produxt
c
call calcu ( nx, aa, bb, cc, md1 )
注：222番のフォーマットは，
計算結果が出力できるように
調整すること。
Input2.dat
##### imax
4
##### iout
2
##### number of variable
2
##### matrix aa
12
23
##### vector bb
1
2
計算後
output2.dat
istep= 2
21.00000
34.00000
istep= 4
377.00000
610.00000
レポート検討課題
（提出日：翌週の授業開始時に回収）
（※以後、本授業のレポートは最終回を除き、翌週の授業開始時に回収する。）
駐車場の料金計算
(1) 30分ごとに150円増加
(2) 6時間を超えて24時間までは1800円で一定
(3) 24時間以降は(1)からの繰り返し
駐車時間x（単位はh）を入力して，
料金yを計算し出力するプログラムを作成せよ。
※計算手順（アルゴリズム）を考え、手順を流れ図（フローチャート）にし、
プログラミングを行いなさい。
研究室のPCにTeraTerm（フリーソフト）をインストールし，ホスト名にnutechと入力すれば
学内の計算機(nutech)にアクセスできるので，授業時間で課題が終了しない場合は，
時間を見つけて実習室あるいは研究室にてレポート課題を行ってください．

Download Report